Milan Korda

Didier Henrion

Colin N Jones

Inner approximations of the region of attraction for polynomial dynamical systems

In a previous work we developed a convex infinite dimensional linear programming (LP) approach to approximating the region of attraction (ROA) of polynomial dynamical systems subject to compact basic semialgebraic state constraints. Finite dimensional relaxations to the infinite-dimensional LP lead to a truncated moment problem in the primal and a polynomial sum-of-squares problem in the dual. This primal-dual linear matrix inequality (LMI) problem can be solved numerically with standard semidefinite programming solvers, producing a hierarchy of outer (i.e. exterior) approximations of the ROA by polynomial sublevel sets, with a guarantee of almost uniform and set-wise convergence. In this companion paper, we show that our approach is flexible enough to be modified so as to generate a hierarchy of polynomial inner (i.e. interior) approximations of the ROA with similar convergence guarantees.

Introduction

Given an autonomous nonlinear system and a target set, the region of attraction (ROA) is the set of all states that end in the target set at a given time without leaving the state constraint set 5 . The ROA is one of the principal sets associated to any dynamical system and goes by many other names in the literature (e.g., backward reachable set or capture basin [START_REF] Aubin | Viability theory: new directions[END_REF]).

In [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF] we showed (in a controlled setting) that there is a genuinely primal convex characterization of the ROA. Optimization over system trajectories is formulated as optimization over occupation measures, leading to an infinite dimensional linear programming (LP) problem in the cone of nonnegative measures. Finite dimensional relaxations of the dual of this problem then provide a converging sequence of outer approximations to the ROA. For a description of alternative techniques for numerical approximations of the ROA, please consult [START_REF] Chesi | Domain of attraction; analysis and control via SOS programming[END_REF] or [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF] and the many references therein.

In this paper we show, within the same measure-theoretic framework, that there exists an infinite dimensional LP whose finite-dimensional relaxations provide a converging sequence of inner approximations to the ROA. This paper can therefore be seen as a complement to [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. To simplify our developments and to emphasize our contribution, we focus only on the uncontrolled setting. The main idea is to construct a converging sequence of outer approximations to the complement of the ROA. There are certain difficulties, topological in nature, associated with this approach. A careful distinction had to be made between trajectories leaving the constraint set and trajectories hitting its boundary. This then translates to a (sometimes subtle, but necessary) distinction between open and closed semialgebraic sets. Fortunately, the LP formulation proposed in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF] was flexible enough to allow for these modifications.

Generally speaking, and consistently with our previous work [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF], we believe that the main virtues of our approach are overall convexity, conceptual simplicity and compactness. Both primal and dual finite-dimensional relaxations turn out to be linear matrix inequalities (LMI), also called semidefinite programming (SDP) problems, with no tuning parameters besides the relaxation order and no initialization data besides the defining ingredients of the problem. In addition, the inner approximations obtained are particularly simplethey are given by a sublevel set of a single polynomial of a predefined degree. Therefore, an ROA approximation in analytic form can be readily obtained by solving a single LMI using freely available software (e.g., SeDuMi [START_REF] Pólik | SeDuMi: a package for conic optimization[END_REF]).

Problem statement

Consider the autonomous system

ẋ(t) = f (t, x(t)), x(t) ∈ X ⊂ R n , t ∈ [0, T ] (1)
with a given vector field f with polynomial entries

f i ∈ R[t, x], i = 1, . . . , n, final time T > 0. The state trajectory x(•) is constrained to a nonempty open 6 basic semialgebraic set 7 X := {x ∈ R n : g X (x) > 0}, (2)
where the polynomial

g X ∈ R[x] is such that the set X := {x ∈ R n : g X (x) ≥ 0} ⊃ X is compact 8 .
The vector field f is polynomial and therefore Lipschitz on the compact set X. As a result, for any x 0 ∈ X there exists a unique maximal solution x(•) to ODE [START_REF] Ambrosio | Transport equation and Cauchy problem for non-smooth vector fields[END_REF]. The time interval on which this solution is defined contains the time interval on which x(t) ∈ X.

Region of attraction (ROA)

Given a final time T and an open bounded basic semialgebraic target set

X T := {x ∈ R n : g T (x) > 0} ⊂ X,
the region of attraction (ROA) is defined as

X 0 := x 0 ∈ X : ∃ x(•) s.t. ẋ(t) = f (t, x(t)), x(0) = x 0 , x(T ) ∈ X T , x(t) ∈ X, ∀t ∈ [0, T ] .
(3) In words, the ROA is the set of all initial states from X for which the unique solution to (1) stays in X for all t ∈ [0, T ] and ends in the target set at time T .

Complement ROA

The idea to get inner approximations of the ROA X 0 is to construct outer approximations of the complement ROA X c 0 := X \ X 0 . By continuity of solutions to (1), the set X c 0 is equal to

X c 0 = x 0 ∈ X : ∃ x(•) s.t. ẋ(t) = f (t, x(t)) and ∃ t ∈ [0, T ] s.t. x(t) ∈ X ∂ and/or x(T ) ∈ X c T ,
where

X c T := {x ∈ R n : g X (x) ≥ 0, g T (x) ≤ 0} is the complement of X T in X and X ∂ := {x ∈ R n : g X (x) = 0}.
In words, X c 0 is the set of initial states that give rise to trajectories which do not end up in X T at time T and/or violate the state constraint at some point between 0 and T .

Occupation measures

In this section we introduce the concept of occupation measures and show how the nonlinear system dynamics can be equivalently described by a linear equation on measures. Notation We will use the following notation. The vector space of all signed Borel measures with support contained in a Borel set K is denoted by M (K). The support (i.e., the smallest closed set whose complement has a zero measure) of a measure µ is denoted by spt µ. The space of continuous functions on K is denoted by C(K) and likewise the space of continuously differentiable functions is C 1 (K). The indicator function of a set K (i.e., the function equal to one on K and zero otherwise) is denoted by I K (•). The symbol λ denotes the n-dimensional Lebesgue measure (i.e., the standard n-dimensional volume). The integral of a function v w.r.t a measure µ over a set K is denoted by K v(x) dµ(x). Sometimes we for simplicity omit the integration variable and/or the set over which we integrate if they are obvious from the context. Now assume x 0 ∈ X and define the first hitting time of X ∂ as

τ (x 0 ) := min T, inf{t ≥ 0 : x(t | x 0 ) ∈ X ∂ } , (4)
where x(• | x 0 ) denotes the unique trajectory starting from x 0 (which is well defined on the time interval [0, τ (x 0 )]). Then we define the occupation measure associated to the trajectory starting from x 0 by

µ(A × B | x 0 ) := τ (x 0 ) 0 I A×B (t, x(t)) dt for all Borel 9 sets A × B ⊂ [0, T ] × X.
The interpretation is that the occupation measure measures the time spent by the trajectory x(• | x 0 ) in subsets of the state space.

The occupation measure enjoys the following important property: for any measurable function v(t, x) the equality

τ (x 0 ) 0 v(t, x(t)) dt = [0,T ]× X v(t, x) dµ(t, x | x 0 ) (5)
holds. In words, the time integral of a function evaluated along the trajectory x(• | x 0 ) is equal to the integral of the function w.r.t. the occupation measure associated to x 0 . Therefore, loosely speaking, all information about the trajectory x(• | x 0 ) is encoded by the occupation measure µ(• | x 0 ). Now suppose that the initial state is not a single point but that its spatial distribution is given by an initial measure µ 0 ∈ M ( X). Then we define the average occupation measure µ ∈ M ([0, T ] × X) as

µ(A × B) := X µ(A × B | x 0 ) dµ 0 (x 0 ).
Lastly, we define the final measure

µ T ∈ M ([0, T ] × X) by µ T (B) := X I B (x(T | x 0 )) dµ 0 (x 0 ).
To derive an equation linking the three principal measures, consider a test function v ∈ C 1 ([0, T ] × X) evaluated along a trajectory. Using the chain rule and equation [START_REF] Chesi | Domain of attraction; analysis and control via SOS programming[END_REF] we

obtain v τ (x 0 ), x(τ (x 0 ) | x 0 ) -v(0, x 0 ) = τ (x 0 ) 0 d dt v(t, x(t | x 0 )) dt = τ (x 0 ) 0 ∂v ∂t + grad v • f (t, x(t | x 0 )) dt = [0,T ]× X ∂v ∂t + grad v • f (t, x) dµ(t, x | x 0 ) = [0,T ]× X Lv(t, x) dµ(t, x | x 0 )
where the linear operator L :

C 1 ([0, T ] × X) → C([0, T ] × X) is defined by v → Lv := ∂v ∂t + grad v • f.
Integrating the above equation w.r.t. µ 0 leads to the equation

[0,T ]× X v(t, x) dµ T (t, x) - X v(0, x) dµ 0 (x) = [0,T ]× X Lv(t, x) dµ(t, x) ∀ v ∈ C 1 ([0, T ] × X), (6)
which is a linear equation linking the measures µ 0 , µ and µ T . Equation ( 6) is sometimes referred to as Liouville's equation.

Primal LP

In this section we follow the approach developed in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF] and derive an infinite-dimensional linear programming (LP) characterization of the complement ROA X c 0 . Certain sublevel sets of feasible solutions to the dual of this LP then yield inner approximations to the ROA X 0 .

The basic idea is to maximize the mass of the initial measure µ 0 under the constraint that it is dominated by the Lebesgue measure, i.e., µ 0 ≤ λ. System dynamics is captured by Liouville's equation ( 6) and state and terminal constraints are handled by suitable constraints on the support of the measures. The key idea is then to split the final measure in two measures such that each measure is supported on a suitable compact basic semialgebraic set. More explicitly, we let

µ T := µ 1 T + µ 2 T with µ 1 T ∈ M ([0, T ] × X ∂ ) and µ 2 T ∈ M ({T } × X c T ). That is, we require that spt µ 1 T ⊂ [0, T ] × X ∂ and spt µ 2 T ⊂ {T } × X c T .
The interpretation is that measure µ 1 T models the trajectories that leave X, whereas measure µ 2

T models the trajectories that end in X c T (i.e., not in X T ). These support constraints on the final measure(s) along with system dynamics enforce that the support of the initial measure µ 0 must be contained in X c 0 . Since there are no other constraints on µ 0 besides µ 0 ≤ λ, maximization of its mass should yield the restriction of the Lebesgue measure λ to X c 0 . The constraint µ 0 ≤ λ can be rewritten equivalently as µ 0 + μ0 = λ for some nonnegative slack measure μ0 ∈ M (X). This is equivalent to requiring that w dµ 0 + w dμ 0 = w dλ for all test functions w ∈ C( X). In addition, we can drop the time argument from the definition of µ 2 T since its time component is supported on a singleton. This leads to the following optimization problem:

p * = sup 1 dµ 0 s.t. v dµ 1 T + v(T, •) dµ 2 T -v(0, •) dµ 0 = Lv dµ ∀ v ∈ C 1 ([0, T ] × X) w dµ 0 + w dμ 0 = w dλ ∀ w ∈ C( X) µ 0 ≥ 0, µ ≥ 0, µ 1 T ≥ 0, µ 2 T ≥ 0, μ0 ≥ 0 spt µ ⊂ [0, T ] × X, spt µ 0 ⊂ X, spt μ0 ⊂ X spt µ 1 T ⊂ [0, T ] × X ∂ , spt µ 2 T ⊂ X c T , (7)
where the supremum is over the vector of nonnegative measures

(µ 0 , µ, µ 1 T , µ 2 T , μ0 ) ∈ M ( X) × M ([0, T ] × X) × M ([0, T ] × X ∂ ) × M (X c T ) × M ( X).
Problem ( 7) is an infinite-dimensional LP in the cone of nonnegative measures. Indeed, the objective is linear, the first two constraints are linear equality constraints and the remaining constraints are conic constraints (the set of nonnegative measures supported on a given set is a positive cone in the vector space of all measures supported on the same set).

The discussion leading to problem [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] is formalized in the following result.

Theorem 1 The optimal value of LP problem ( 7) is equal to the volume of the complement ROA X c 0 , that is, p * = λ(X c 0 ). Moreover, the supremum is attained by the restriction of the Lebesgue measure to the complement ROA X c 0 .

Proof: Closely follows arguments in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. By definition of the relaxed complement ROA, the unique trajectory x(•) associated to any initial condition x 0 ∈ X c 0 either hits X ∂ at some t ∈ [0, T ] or ends in X c T . Therefore for any initial measure µ 0 ≤ λ with spt µ 0 ⊂ X ⊂ X there exist an occupation measure µ, final measures µ 1 T , µ 2 T and a slack measure μ0 such that the constraints of problem (7) are satisfied. One such measure µ 0 is the restriction of the Lebesgue measure to X c 0 , and therefore p * ≥ λ(X c 0 ). Now we show that p * ≤ λ(X c 0 ). Take a vector of measures (µ 0 , µ, µ 1 T , µ 2 T , μ0 ) feasible in [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] and suppose that λ(spt µ 0 \ X c 0 ) > 0. Since any level set of a polynomial has a zero Lebesgue measure we have λ(X ∂ ) = 0 and

λ(spt µ 0 \ (X c 0 ∪ X ∂ )) = λ(spt µ 0 \ X c 0 ) > 0.
By a superposition principle [1, Theorem 3.2] using arguments of [6, Appendix A, Lemma 4], there exists a family of admissible trajectories of the ODE (1) starting from µ 0 generating the occupation measure µ and the final measure µ T = µ 1 T + µ 2 T . However, this is a contradiction since spt µ 0 \ (X c 0 ∪ X ∂ ) ⊂ X 0 , which means that all trajectories starting from spt µ 0 \ (X c 0 ∪ X ∂ ) neither hit X ∂ nor end in Xc T . Thus, λ(spt µ 0 \ X c 0 ) = 0 and so λ(spt µ 0 ) ≤ λ(X c 0 ). Combining this with the constraint µ 0 ≤ λ we get µ 0 (X) = µ 0 (spt µ 0 ) ≤ λ(spt µ 0 ) ≤ λ(X c 0 ) for any feasible µ 0 . Therefore p * ≤ λ(X c 0 ) and thus in fact p * = λ(X c 0 ).

Dual LP

In this section we derive a dual LP on continuous functions, prove the absence of a duality gap between the primal and dual LPs and relate feasible solutions to the dual to the indicator function of the complement ROA X c 0 . By standard infinite-dimensional LP theory (see, e.g., [START_REF] Anderson | Linear programming in infinite-dimensional spaces: theory and applications[END_REF]), the dual to LP [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] reads

d * = inf X w(x) dλ(x) s.t. Lv(t, x) ≤ 0, ∀ (t, x, u) ∈ [0, T ] × X w(x) ≥ v(0, x) + 1, ∀ x ∈ X v(T, x) ≥ 0, ∀ x ∈ X c T v(t, x) ≥ 0, ∀(t, x) ∈ [0, T ] × X ∂ w(x) ≥ 0, ∀ x ∈ X, (8)
where the infimum is over (v, w)

∈ C 1 ([0, T ] × X) × C( X).
The intuition is that given x 0 ∈ X c 0 the constraint Lv ≤ 0 forces v to decrease along trajectories as long as it does not hit X ∂ or end in

X c T . Because of the constraint v ≥ 0 on [0, T ] × X ∂ ∪ {T } × X c T we must have v(0, •) ≥ 0 on X c 0 . Consequently, w(x) ≥ 1 on X c
0 . This instrumental observation is formalized in the following Lemma.

Lemma 1 If Lv ≤ 0 on [0, T ]× X, v ≥ 0 on ([0, T ]×X ∂ )∪({T }×X c T ) and w ≥ v(0, •)+1 on X, then w ≥ 1 on X c 0 .
Proof: Take x 0 ∈ X c 0 and consider the first hitting time of X ∂ , τ := τ (x 0 ), defined by [START_REF] Aubin | Viability theory: new directions[END_REF]. By definition of X c 0 the trajectory starting from x 0 will either hit X ∂ or end in

X c T . Therefore x(τ ) ∈ ([0, T ] × X ∂ ) ∪ ({T } × X c T ) and x(t) ∈ X for t ∈ [0, τ ]. Therefore v(τ, x(τ )) ≥ 0, Lv(t, x(t)) ≤ 0, ∀ t ∈ [0, τ ] and so 0 ≤ v(τ, x(τ )) = v(0, x 0 ) + τ 0 Lv(t, x(t)) dt ≤ v(0, x 0 ) ≤ w(x 0 ) -1.
The following result is of key importance for subsequent developments.

Theorem 2 There is no duality gap between primal LP problems (7) on measures and dual LP problem (8) on functions, that is, p * = d * .

Proof: Here we only outline the basic steps; for a detailed argument in a similar setting see [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]Theorem 2]. Since the supports of all measures are compact, the initial measure is dominated by the Lebesgue measure and the final time is finite, we have

µ 0 ( X) ≤ λ( X) < ∞, µ T ([0, T ] × X) = µ 0 ( X) < ∞ and µ([0, T ] × X) ≤ T µ T ([0, T ] × X) < ∞,
where the last two inequalities follow by plugging in v(t, x) = 1 and v(t, x) = t in Liouvillel's equation [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. Therefore p * < ∞ and the feasible set of problem ( 7) is weakly-* bounded. Furthermore, the feasible set of ( 7) is nonempty since (µ 0 , µ, µ 1 T , µ 2 T , μ0 ) = (0, 0, 0, 0, λ) is a trivial feasible point; therefore 0 ≤ p * < ∞. The absence of a duality gap then follows from [2, Theorem 3.10] using Alaoglu's theorem (see, e.g., [START_REF] Luenberger | Optimization by vector space methods[END_REF]Chapter 5]) and the weak-* continuity of the adjoint of the operator L.

Next, we establish our first convergence10 result.

Theorem 3 There is a sequence of feasible solutions to problem [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] such that its w-component converges from above to I X 0 in L 1 norm and almost uniformly.

Proof: Follows by the same arguments as Theorem 3 in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF].

LMI relaxations

In this section we derive finite dimensional semidefinite programming (SDP) or linear matrix inequality (LMI) relaxations to the infinite dimensional LPs ( 7) and ( 8) and establish several convergence results relating these relaxations to the infinite dimensional LPs and to the ROA.

In what follows, R k [•] denotes the vector space of real multivariate polynomials of total degree less than or equal to k.

Derivation of the finite dimensional relaxations is standard and the reader is referred to [6, Section 5] or to the comprehensive reference [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]; therefore we only highlight the main ideas. First of all, since the supports of all measures feasible in ( 7) are compact, these measures are determined by their moments, i.e., by integrals of all monomials (which is a sequence of real numbers when indexed in, e.g., the canonical monomial basis). Therefore, it suffices to restrict the test functions w(x) and v(t, x) in ( 7) to all monomials, reducing the linear equality constraints of ( 7) to linear equality constraints on the moments. Next, by the celebrated Putinar Positivstellensatz (see [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF][START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF]), the constraint that the support of a measure is included in a given compact basic semialgebraic set is equivalent to the feasibility of an infinite sequence of LMIs involving the so-called moment and localizing matrices, which are linear in the coefficients of the moment sequence. By truncating the moment sequence and taking only the moments corresponding to monomials of total degree less than or equal to 2k we obtain a necessary condition for this truncated moment sequence to be the first part of a moment sequence of a measure with the desired support. This procedure leads to the primal SDP relaxation of order k

p * k = max (y 0 ) 0 s.t. A k (y, y 0 , y 1 T , y 2 T , ŷ0 ) = b k M k (y) 0, M k-d X i (g X , y) 0 M k (y 0 ) 0, M k-d X (g X , y 0 ) 0 M k (y 1 T ) 0, M k-d T (g X , y 1 T ) 0 M k-d T (-g X , y 1 T ) 0 M k (y 2 T ) 0, M k-d T (g X , y 2 T ) 0 M k-d T (-g T , y 2 T ) 0 M k (ŷ 0 ) 0, M k-d X (g X , ŷ0 ) 0 M k-1 (t(T -t), y) 0, M k-1 (t(T -t), y 1
T ) 0 (9) where the notation 0 stands for positive semidefinite and the minimum is over moment sequences (y, y 0 , y 1 T , y 2 T , ŷ0 ) truncated to degree 2k corresponding to measures µ, µ 0 , µ 1 T , µ 2

T and μ0 . The linear equality constraint captures the two linear equality constraints of ( 7) with v(t, x) ∈ R 2k [t, x] and w(x) ∈ R 2k [x] being monomials of total degree less than or equal to 2k. The matrices M k (•) are the moment and localizing matrices, following the notations of [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] or [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. In problem [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF], a linear objective is minimized subject to linear equality constraints and LMI constraints; therefore problem ( 9) is an SDP problem.

The SDP problem dual to problem (9) turns out to be the sum-of-squares problem

d * k = inf w • l s.t. -Lv(t, x) = p(t, x) + q 1 (t, x)t(T -t) + q 2 (t, x)g X (x) w(x) -v(0, x) -1 = p 0 (x) + q 0 1 (x)g X (x) v(t, x) = p T 1 (x) + q T 1 (t, x)t(T -t) + r(x)g X (x) v(T, x) = p T 2 (x) + q T 2 (x)g X (x) -q T 3 (x)g T (x) w(x) = s 0 (x) + s 1 (x)g X (x), ( 10
)
where l is the vector of Lebesgue moments over X indexed in the same basis in which the polynomial w(x) with coefficients w is expressed. The minimum is over polynomials

v(t, x) ∈ R 2k [t, x], w(x) ∈ R 2k [x]
, over the polynomial r(x) and polynomial sum-of-squares p(t, x), q 1 (t, x), q 2 (t, x), q 0 1 (x), p T 1 (x), p T 2 (x), q T 1 (x), q T 2 (x), q T 3 (x), s 0 (x), s 1 (x) of appropriate degrees. The constraints that polynomials are sum-of-squares can be written explicitly as LMI constraints (see, e.g., [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]), and the objective is linear in the coefficients of the polynomial w(x); therefore problem [START_REF] Luenberger | Optimization by vector space methods[END_REF] can be formulated as an SDP problem.

Theorem 4 There is no duality gap between primal LMI problem ( 9) and dual LMI problem [START_REF] Luenberger | Optimization by vector space methods[END_REF], i.e. p * k = d * k .

Proof: Follows by the same arguments based on standard SDP duality theory as Theorem 4 in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF].

Now we prove our main convergence results.

Theorem 5 Let w k ∈ R 2k [x] denote the w-component of a solution to the dual LMI problem [START_REF] Luenberger | Optimization by vector space methods[END_REF] and let wk (x) = min i≤k w i (x). Then 1 -w k converges from below to I X 0 in L 1 norm and 1 -wk converges from below to I X 0 in L 1 norm and almost uniformly.

Proof: It follows from Theorem 3 and from the density of polynomials in the space of continuous functions on compact sets (for a detailed argument in a similar setting see [6, Theorem 5]) that w k and wk converge from above to I X c 0 in L 1 and almost uniformly on X, respectively. Therefore 1 -w k and 1 -wk converge from below to I X 0 = 1 -I X c 0 on X in the same manner.

The next Corollary follows immediately from Theorem 5.

Corollary 1

The sequence of infima of LMI problems [START_REF] Luenberger | Optimization by vector space methods[END_REF] Proof: Follows the proof of Corollary 1 in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. Monotone convergence of the dual optima d * k follows immediately from Theorem 5 and from the fact that the higher the relaxation order k, the looser the constraint set of the minimization problem [START_REF] Luenberger | Optimization by vector space methods[END_REF]. To prove convergence of the primal maxima observe that from weak SDP duality we have d * k ≥ p * k and from Theorems 5 and 2 it follows that d * k → d * = p * . In addition, clearly p * k ≥ p * and p * k+1 ≤ p * k since the higher the relaxation order k, the tighter the constraint set of the maximization problem [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. Therefore p * k → p * monotonically from above. Our last results establishes set-wise convergence of inner approximations to the ROA.

Theorem 6 Let w k ∈ R 2k [x]
denote the w-component of a solution to the dual LMI problem [START_REF] Luenberger | Optimization by vector space methods[END_REF] and let X 0k := {x ∈ X :

w k (x) < 1}. Then X 0k ⊂ X 0 , lim k→∞ λ(X 0 \ X 0k ) = 0 and λ(X 0 \ ∪ ∞ k=1 X 0k ) = 0.
Proof: Follows the proof of Theorem 6 in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. From Lemma 1 we have w k (x) ≤ 1 ⇒ x ∈ X 0 for all x ∈ X, and therefore I X 0k ≤ I X 0 . Since w k ≥ 0 on X we also have 1 -w k ≤ I X 0k on X and therefore 1 -

w k ≤ I X 0k ≤ I X 0 on X. From Theorem 5, we have 1 -w k → I X 0 in L 1 norm on X. Consequently, λ(X 0 ) = X I X 0 dλ = lim k→∞ X 1 -w k dλ ≤ lim k→∞ X I X 0k dλ = lim k→∞ λ(X 0k ) ≤ lim k→∞ λ(∪ k i=1 X 0i ) = λ(∪ ∞ k=1 X 0k ).
But since X 0k ⊂ X 0 for all k, we must have

lim k→∞ λ(X 0k ) = λ(X 0 ) and λ(∪ ∞ k=1 X 0k ) = λ(X 0 ),
which proves the theorem.

Numerical examples

In this section we present two numerical examples. The primal problems on measures were modeled using Gloptipoly 3 [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] interfaced with the SDP solver SeDuMi [START_REF] Pólik | SeDuMi: a package for conic optimization[END_REF]; this solver also returns the solution to the dual SDP relaxation. In Section 7.3 we then investigate how tight low order approximations can be obtained.

Univariate cubic dynamics

Consider the system given by ẋ = x(x -0.5)(x + 0.5), the constraint set X = [-1, 1], the final time T = 10 and the target set

X T = [-0.3, 0.3].
The ROA in this setup is X 0 = [-0.5, 0.5]. Polynomial approximations to the complement ROA for degrees d ∈ {16, 20, 24, 28} are shown in Figure 1. As expected, the functional convergence of the polynomial to the discontinuous indicator function is rather slow. A slightly better convergence is observed in the volume error of the sublevel set approximation to the ROA documented in Table 1. The relatively slow convergence could be significantly improved if a tighter constraint set X was employed; see Section 7.3 below. Alternative polynomial bases (e.g. Chebyshev polynomials) would also allow tighter higher order approximations; see [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] for more details.

-

Van der Pol oscillator

As a second example consider a scaled11 version of the uncontrolled reversed-time Van der Pol oscillator given by ẋ1 = -2.0x 2 , ẋ2 = 0.80x 1 + 10(x 2 1 -0.21)x 2 .

We take T = 1 and X T = {x ∈ R n : x 2 ≤ 0.50} and X := {x ∈ R n : x ≤ 1.1}. The ROA is bounded, having the characteristic Van der Pol shape. Plots of polynomial sublevel set approximations of degrees d ∈ {9, 12, 15, 18} are shown in Figure 3. We observe a relatively fast convergence to the ROA, which is also documented by the relative volume errors reported in Table 2. Figure 2 then shows a degree 18 polynomial approximation to the indicator function of the complement ROA.

Low order approximations

In the examples above, relatively high order polynomials had to be used to obtain tight approximations, which can limit subsequent applicability of the approximations. There are several ways to obtain low order approximations of similar quality. First of all, since the integral of a polynomial w is minimized over the constraint set X, it is desirable that X be a good outer approximation of the ROA. Of course, selecting X is possible only if it is an artificially specified outer approximation of the ROA, not a constraint set coming from physical requirements on the system. More importantly, notice that in problem (10) the system dynamics enters the constraints on the polynomial v(t, x), whereas the polynomial w(x) is only upper-bounding v(t, x) + 1 for t = 0. Since the inner approximations are given by sublevel sets of w, it is possible and plausible to choose different degrees of w and v -low for w and higher for v. Both techniques are illustrated in Figure 4; in Figure 4 (a) we consider the univariate cubic dynamics and we both shrink the constraint set X and choose low order w while keeping v of higher order. In Figure 4 (b) we consider the Van Der Pol oscillator, keeping the constraint set X unchanged and only selecting low order w. The inner approximations obtained are indeed significantly tighter for the given degrees (compare with Figures 1 and3). 1). Right: low order w only (compare with Figure 3).

Conclusion

This paper presented an infinite dimensional convex characterization of the region of attraction (ROA) for uncontrolled polynomial systems, following the approach initiated in our previous work [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. Finite dimensional dual relaxations yield a converging sequence of inner approximations to the ROA, thereby complementing the outer approximations of [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. One of the virtues of the approach is its conceptual simplicity -the resulting approximation is the outcome of a single SDP or LMI problem with no free parameters except for the relaxation order. The approximations itself are also simple, given by sublevel sets of polynomials of predefined degrees.

Nevertheless, this approach does not escape the curse of dimensionality -indeed, whereas the number of variables of the LMI relaxations grows polynomially with the relaxation order, this number grows exponentially with the state dimension. Tailored structureexploiting SDP solvers could enable this approach to reach higher dimensions. In addition, a different choice of basis functions (e.g., Chebyshev polynomials rather than monomials) would improve numerical conditioning of the LMIs, allowing higher oder relaxations to be computed.

Future research directions include inner approximations in a controlled setting and the related problem of robust region of attraction / reachable set computation with either unknown but constant uncertainty or a time-varying disturbance. The cases of asymptotic region of attraction and maximum (controlled) positively invariant set computation are amenable to similar tools.

  converges monotonically from above to the supremum of the LP problem (8), i.e., d * ≤ d * k+1 ≤ d * k and lim k→∞ d * k = d * . Similarly, the sequence of maxima of LMI problems (9) converges monotonically from above to the maximum of the LP problem (7), i.e., p * ≤ p * k+1 ≤ p * k and lim k→∞ p * k = p * .

Figure 1 :

 1 Figure 1: Univariate cubic dynamics -polynomial approximations (solid line) to the complement ROA indicator function I X c 0 = I [-1,-0.5] + I [0.5,1] (dashed line) for degrees d ∈ {16, 20, 24, 28}.

Figure 2 :

 2 Figure 2: Van der Pol oscillator -degree 18 polynomial approximation to the indicator function of the complement ROA.

Table 2 :

 2 Van der Pol oscillator -relative error of the inner approximation to the ROA X 0 as a function of the approximating polynomial degree.

Figure 3 :

 3 Figure 3: Van der Pol oscillator -polynomial inner approximations (light gray) to the ROA (dark gray) for degrees d ∈ {9, 12, 15, 18}.

  Van der Pol oscillator -deg w = 8 (deg v = 18). Volume approximation error 5.46 %.

Figure 4 :

 4 Figure4: Low order approximations to the ROA. Left: tighter constraint set X and low order w (compare with Figure1). Right: low order w only (compare with Figure3).

Table 1 :

 1 Univariate cubic dynamics -relative error of the inner approximations to the ROA X 0 = [-0.5, 0.5] as a function of the approximating polynomial degree.

	degree	16	20	24	28
	error 11.4 % 6.4 % 4.84 % 4.54 %

Laboratoire d'Automatique, École Polytechnique Fédérale de Lausanne, Station 9, CH-1015, Lausanne, Switzerland. {milan.korda,colin.jones}@epfl.ch

2 CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse; France. henrion@laas.

fr 3 Université de Toulouse, LAAS, F-31400 Toulouse;

France 4 Faculty of Electrical Engineering, Czech Technical University in Prague, Technická 2, CZ-16626 Prague, Czech

Republic 5 There are various modifications on this setup (e.g., one may consider asymptotic convergence instead of finite-time reachability, with or without constraints, in the presence of disturbances and/or uncertainty, or in a controlled setting); most of these modifications are amenable to the methods presented in this paper, sometimes with different qualitative results.

The requirement of the constraint set being open is merely technical, for this considerably simplifies the developments and the proofs.

For clarity of exposition we consider the constraint set given by a single superlevel set of a polynomial. The approach can, however, be straightforwardly extended to constraint sets defined by the intersection of finitely many polynomial superlevel sets.

Note that the closed semialgebraic set X = {x : g X (x) ≥ 0} can be strictly larger than the closure of the open semialgebraic set X = {x : g X (x) > 0}, consider in R e.g. g X (x) = (1 -x 2 )(2 + x)[START_REF] Anderson | Linear programming in infinite-dimensional spaces: theory and applications[END_REF] . For a similar reason, note also that X bounded does not imply X bounded. Indeed, in R 2 with g X (x) = (1 -x 2 1 -x 2 2 )(2 + x 1 ) 2 we have X = {x : ||x|| < 1} and X = {x : ||x|| ≤ 1} ∪ {x : x 1 = -2}.

For brevity we drop the adjective "Borel" in the sequel.

Please refer to[START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF] or, e.g.,[START_REF] Ash | Real analysis and probability[END_REF] for definitions of the various types of convergence relevant in this context.

The coefficients were chosen so that the ROA fits within the box [-1, 1] 2 .


Pascal Berto
email: pascal.berto@fresnel.fr

Esteban Bermudez Ureña

Pierre Bon

Romain Quidant

Hervé Rigneault

Guillaume Baffou
email: guillaume.baffou@fresnel.fr

Esteban Bermúdez Ureña

Quantitative absorption spectroscopy of nano-objects

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

The likelihood of interaction between light and a nanoparticle (NP) can be quantified by different optical cross sections. The absorption cross section (ACS) σ abs quantifies the ability of the NP to absorb an incident photon, while the scattering cross section σ sca refers to the ability of the NP to scatter a photon. These cross sections are defined such that the power P abs absorbed by the NP and the power P sca scattered by the NP simply read P abs = σ abs I and P sca = σ sca I, where I is the irradiance (power per unit area) of the incoming light. The total (or extinction) cross section is σ ext = σ sca + σ abs . Depending on the nature of the NP, its morphology, and the excitation wavelength, the relative efficiencies of absorption and scattering can be very different from one NP to another.

Light absorption by a nano-object is often considered as a source of energy loss since it usually corresponds to a conversion of light into heat. However, since recent developments in thermo-plasmonics, 1 this transduction pathway has been promoted since it makes it possible to use a metal NP as a nano-source of heat, one of the most fundamental tool in nature. 2 Many emerging research areas and applications are now based on photothermal processes of noble metal NPs. In this field, the knowledge of the maximum absorption wavelength and the corresponding ACS are important to optimize heat generation. However, measuring the absorption cross section of a nano-object remains a difficult task, mainly because the information of interest is no longer radiated and detectable in the far field. Several techniques have been developed recently to estimate ACSs of metal NPs. [3][4][5] However, they present major limitations: i) They rely on the assumption that ACS can be assimilated to extinction cross section. This assumption is valid only for NPs much smaller than the illumination wavelength. ii) They do not lead to the absolute value of the ACS. iii) They use a complex instrumentation. iv) They are relatively slow.

We report here on an optical microscopy technique able to quantitatively and rapidly measure the ACS of absorbing nano-structures, whatever their nature, their size and their morphology. After introducing and describing the principle of the technique, we illustrate its capabilities first with an ensemble measurement on a regular pattern of lithographic gold structures, and then with isolated nanostructures. Then, we discuss the sensitivity of the measurements and explain how to optimize it. A final part is dedicated to make a comparison with previous related techniques reported in the literature.

II. EXPERIMENTAL TECHNIQUE

We recently introduced an optical microscopy technique able to map the distributions of temperature and heat power density around metal nanostructures. 6 This technique is based on quadriwave lateral shearing interferometry (QWLSI). Our present work follows this development and utilizes a similar experimental setup, which is presented in Fig. 1, to achieve quantitative absorption spectroscopy. The principle consists in probing the thermally induced distortion of a planar optical wavefront passing through the medium of interest. The wavefront distortion is due to light absorption by the NP that generates a temperature increase and a subsequent variation of the refractive index of the surrounding medium. The samples we used consist of lithographic gold nanostructures lying on a planar glass substrate. A drop of glycerol was deposited on this substrate and an additional glass cover slip was pressed on top of the sample to end up with a sandwich-like structure (inset of Fig. 1). The use of a liquid as the surrounding medium ensures a substantial thermally induced variation of the refractive index even for temperature increases as small as 1 K. The sample is illuminated from the top by a collimated broad-band light emitting diode (λ 0 = 600 ± 50 nm). The distorted wavefront profile of this probe beam is then imaged using the QWLSI wavefront sensor. More precisely, a modified Hartmann diffraction grating 7 produces an interferogram on a CCD camera that can be processed to quantitatively retrieve the optical wavefront profile of the probe beam. 6 In parallel, a CW Ti:Sapphire laser is used to heat the light absorbing metal structures. The wavelength of this infrared (IR) laser can be tuned from λ 0 =720 up to 950 nm. We designed the gold nanostructures to adjust their plasmonic resonances in this wavelength range. To perform an image, a first reference interferogram is recorded without heating (IR laser turned off) from which a reference wavefront profile is computed. Then, the IR laser is turned on to heat the structures and a second interferogram is recorded which yields a second wavefront profile. From the subtraction of these two images, the distribution of the optical path difference (OPD) through the sample experienced by the probe beam is retrieved.

In the case of a point-like NP, the thermally induced OPD can be easily modeled theoretically. In such an ideal case, the temperature profile in the surrounding liquid is simply given by:

T (r) = T ∞ + σ abs I 4πκ r ( 1
)
where r is the radial coordinate, T ∞ the ambient temperature, I the irradiance of the illumination and κ = (κ 1 + κ 2 )/2 the average of the thermal conductivities of the surrounding liquid and the substrate. Considering the average thermal conductivity κ is not an approximation. It is the exact solution of the problem, which can be derived using the image method with a source of heat localized right at the interface. 8 An analytical expression of this associated distorted wavefront can be easily derived if one considers that the refractive index of glycerol varies linearly with the temperature, δn = β(T -T ∞ ). By integrating δn over the liquid height we find that the OPD δℓ experienced by the beam reads:

δℓ(r) = P 4πκ β ln h/r + 1 + (h/r) 2 (2) δℓ(r) = P G ℓ (r) (3)
where P = σ abs I. Eq. 3 defines the Green's function G ℓ discussed further in this work. Eq. 2 can be directly used to quantitatively estimate the delivered power σ abs of a point-like NP since the liquid thickness h and the irradiance I can be easily estimated, and the liquid parameters κ and β are presumably known. One just has to fit the radial profile of the OPD with σ abs as the unique fitting parameter. This technique is not limited to retrieve the ACS of point-like NPs. In the case of an extended absorbing structure of any shape and size, the formalism can be refined to retrieve as well the overall ACS. The OPD profile in the surrounding medium appears in this case as a convolution between the heat power density p(x, y) of the sample and the Green's function G ℓ (Eq. 3):

δℓ(r) = [p ⊗ G ℓ ](r) (4)
Hence, the distribution of p(x, y) can be retrieved by a deconvolution procedure. Once p(x, y) is known, the overall ACS of the extended pattern reads

σ abs = 1 I p(x, y)dxdy (5)

III. RESULTS AND DISCUSSION

We shall now illustrate the technique by performing quantitative absorption spectroscopy in the two above mentioned cases: First an extended distribution of nanostructures and then a single isolated NP.

Figure 2 presents a series of measurements performed on an array of identical lithographic gold dimer structures. Each dimer consists of two rods of length L and width 50 nm separated by a gap of 30 nm (Fig. 2a). Several arrays characterized by different rod lengths L have been fabricated and investigated. The central part of the array was illuminated by a uniform circular beam of diameter D = 14 µm to generate optical absorption and heat. The liquid height was set to h = 1.0 mm. Measurements were performed at various wavelengths close to the resonance of the nanostructures in order to retrieve the absorption spectra. Figure 2a presents a SEM image of an array of gold dimers with L = 320 nm. The size of the illuminated area is represented by a circle. At each wavelength, the distribution of heat source density p(x, y) (Fig. 2c) is retrieved from the distribution of the thermal-induced OPD (Fig. 2b) using Eq. 4 and a deconvolution algorithm. The total absorbed power P is estimated quantitatively by spatially integrating the heat power density (Eq. 5). Since the extended structure is composed of identical NPs, we can get a first estimation of the single-NP ACS by considering the quantity σabs = P/(I N ) where N is the estimated number of NPs in the illuminated area. This procedure was repeated at various wavelengths to get a quantitative NP absorption spectra. Different nanostructures with different sizes have been investigated as represented in Fig. 2d. In this series of measurements, we wish to illustrate that different absorption modes can be identified. The resonance observed in Fig. 2d corresponds to the dipolar resonance of the structure. When increasing the length L of the structure, this resonance red-shifts and disappears far in the IR while a second resonance arrives from the visible range of the spectrum. This resonance corresponds to a 3λ 0 /2 plasmonic mode of the structure, 9 as represented in Fig. 2e for different L. Finally, while increasing even more the length of the structure L, the weak 5λ 0 /2 resonance can also be observed and quantified (Fig. 2f). In all the measured absorption spectra, we also plot the associated extinction spectra measured using conventional means. In the case of the dipolar resonance, a slight shift of 20 nm between the absorption and the extinction spectra is observed, which is the expected behavior according to the literature. 10 We shall now describe absorption spectroscopy performed on single NPs. Figure 3 presents a series of measurements performed on single dimer nanostructures (Fig. 3a), as the ones investigated in the previous paragraph. For this measurement, the liquid thickness was set to h = 7 µm and the illumination diameter D = 6.2 µm. This time, the OPD distribution is much more confined as observed in (Fig. 3b). As mentioned before, the estimation of the ACS is much simpler for single NPs since the computation of the heat power density distribution (using a deconvolution procedure) is not required. The heat power can be directly retrieved from the OPD image. The OPD radial profile is simply fitted using Eq. 2 with the delivered power P as the only fitting parameter. A perfect matching between Eq. 2 and the OPD radial profile is usually observed as shown in Fig. 3c. Then, the ACS is simply given by σ abs = P/I. The procedure was repeated for various wavelengths, which gave the absorption spectra represented in Fig. 3d. In this figure, the absorption spectra obtained in the previous section from an array of NPs is also represented and is, as expected, identical. The measured extinction spectra is also plotted, showing the same resonance wavelength. At the current stage of its development, ACSs as small as ca. 10 3 nm 2 can be measured. Using a liquid with a higher refractive index dependence on temperature β and a CCD camera with a faster readout (in order to collect more photons) could improve the sensitivity of the method in the future.

We wish to underline a limitation that can occur in specific conditions while performing absorption measurements. While investigating single NPs, it may happen that the energy absorbed by the illuminated surrounding liquid becomes more important than the energy absorbed by the NP itself. In this case, the temperature increase is also due to liquid heating and the wavefront distortion is no longer a specific signature of the NP. Absorption spectra of different liquids were measured using a standard photospectrometer. Results are represented in Fig. 4a and show that strong absorption usually occur after 850 nm for usual solvents such as water, glycerol and ethanol.

The occurrence of such a side effect can be clearly identified on the OPD image as illustrated in Fig. 4b where two OPD profiles are represented: A first OPD profile taken at λ 0 = 720 nm and a second at λ 0 = 950 nm where liquid absorption dominates. In the latter case, the heated liquid layer creates a linear slope on the OPD profile at high radial coordinate r that cannot be fitted using Eq. 2. This abnormal slope is a sign of liquid heating. Hence, absorption spectroscopy in the IR would require the use of a non-absorbing solvent. The importance of such an effect can be predicted by estimating the order of magnitudes of the powers absorbed by the NP itself P NP and by the liquid layer P liq . It yields:

ξ = P liq P NP = SI(1 -e αh ) σ abs I ≃ αSh σ abs (6)
where S = πD 2 /4, D is the diameter of the beam (assumed to be cylindrical over the liquid height) and α the extinction coefficient of the liquid. According to this expression, and as expected, the influence of the liquid heating can be minimized (ξ ≪ 1) by reducing the height h of the liquid or the diameter D of the incoming beam. In the experiments reported in Fig. 2, α ≃ 0.018 cm -1 and σabs = 10 3 nm 2 at λ 0 = 800 nm. Given that 152 NPs where illuminated it yields ξ ≃ 0.02, which is weak enough. On an isolated NP, like in the experiments reported in Fig. 3, such values of D and h would have led to ξ ≃ 2.7 > 1. This is why we chose to reduce the height of the liquid layer down to h = 7 µm and the beam diameter to D = 6.2 µm. In such conditions, ξ = 0.004 ≪ 1 at λ 0 = 800 nm and ξ remained small up to λ 0 = 850 nm.

We wish finally to compare our technique to the other reported techniques aimed at measuring ACSs, in order to clarify the advantages and the drawbacks of each approaches.

Let us first depict an important class of photothermal imaging (PTI) techniques, introduced in a pioneer article by Orrit, Lounis and coworkers in 2002. 11 In this article, the authors report on a method able to detect single nanoparticles smaller than 3 nm by optically heating them. Such an achievement combined high-frequency modulation and a polarization interference contrast (DIC) measurement. In 2004, Berciaud, Lounis and coworkers made this approach even simpler and more sensitive by implementing a heterodyne detection and measuring a modulated scattered intensity. 12 This improved version of the technique was then used by the same group to qualitatively measure absorption spectra of nano-objects such as metal NPs, carbon nanotubes 13 or semiconductor nano-crystals. 14 In 2010, Gaiduk, Orrit and coworkers optimized the sensitivity of this heterodyne technique and managed to detect single dye molecules, with an ACS as weak as 1 nm 2 5 This series of reported experiments is closely related to our QWLSI technique since it is based on the same physics: a photothermal induced variation of the refractive index surrounding the absorbing structure. The fundamental differences is that we do not modulate the heating and we do not perform the same kind of measurements: While the PTI techniques measure either a modulated scattered intensity, 5,13 or a combination of both the modulated intensity and phase gradient by DIC, 11 the QWLSI technique measures a quantitative 2D phase image. 6 Minding the phase and not the scattered inten-sity leads us to the possibility to quantitatively estimate the heat power delivered by the NP. Indeed, in the case of the PTI technique, a direct estimation of the ACS seems not possible as all the ACS measurements reported by the Berciaud, Lounis and coworkers were in arbitrary units. 13,14 In ref. 5, Gaiduk et al. overcame this problem by previously measuring the photothermal signal of a reference gold spherical nanoparticle of a given size, in the same environment, and computing the theoretical ACS using Mie theory. 5 This gave them a calibration coefficient linking ACS and photothermal signal, assuming a linearity of the signal. Then, they performed measurements on single molecules, measured the photothermal signal, and used the calibration coefficient to indirectly retrieve an absolute value of the ACS of the molecule. The advantage of our technique is to perform a direct measurement of the ACS, in the sense that we do not need to perform a preliminary measurement on a model system (like a gold spherical NP) before investigating the absorbers of interest.

Other experimental approaches are based on the assumption that ACS can be assimilated to extinction cross section, like the approach reported by Arbouet et al.. 3 An extinction cross section is indeed simpler to measure. This assumption is valid only for very small NPs when scattering is negligible, for example for gold NPs smaller than 30 nm. Hence, such techniques cannot be applied, for instance, to larger structures or to extended distributions such as the ones reported in the present work.

Moreover, the QWLSI technique is straightforward to implement on any microscope, since one just has to use a CCD-based device. The use of an acousto-optic modulator, a single photon counting module, a heterodyne detection or to scan the sample to get an image, like the approaches listed above, is not required.

Finally, the QWLSI technique is faster than the above mentioned techniques. A entire image is obtained in less than 1 second while other PTI techniques require to scan an image with an integration time of approx. 10 ms to 200 ms per pixel. For this reason, our approach can be conveniently used to search and find out the resonance wavelength of an absorbing structure, in a reasonable time. The direct consequence is that our technique remains approximatively 10 3 times less sensitive than other techniques that can detect single molecules. Our technique stands for this reason as a complementary technique compared to the existing field of techniques, a technique valuable when the ACS of interest are not too weak, typically for lithographic metal nanostructures larger than 100 nm.

IV. CONCLUSION

In summary, we recently introduced an optical microscopy technique able to map the distributions of temperature and heat power density around metal nanostructures. 6 This technique is based on quadriwave lateral shearing interferometry (QWLSI). Here, we follow this development and detail how to utilize a similar experimental procedure to achieve quantitative absorption spectroscopy on single absorbers. At the current stage of its development, our approach can measure absorption cross sections as small as 10 3 2 . We decided to illustrate the technique on gold nanoparticles featuring plasmonics resonances. Different plasmonics modes were evidenced and quantified on a gold dimer. This sensitivity could be improved by modifying the liquid environment or using a faster CCD camera to collect more photons. We also evidenced that single NP absorption measurements can be affected by the absorption properties of the environment, especially when working under infrared illumination. The environment contribution can be minimized by adjusting geometrical and illumination parameters. This technique is straightforward to implement on any conventional microscope. The emerging field of thermo-plasmonics, which uses metal nanoparticles as nano-sources of heat, should benefit from this technique since the conception of highly absorbing nanostructures is still a subject of active investigations.

FIG. 1 .

 1 FIG. 1. Experimental configuration used to perform quantitative absorption spectroscopy measurements. The inset illustrates the wavefront distortion thermally induced by the heated NP.

FIG. 2

 2 FIG.2. a) Scanning Electron Microscopy (SEM) image of an array of lithographic gold dimers (rod length L = 320 nm) made by e-beam lithography. b) Measured distribution of the thermally-induced optical path difference (OPD) when heating the array with a uniform circular beam at λ0 = 750 nm. c) Associated heat power density computed from image b. d) Experimental absorption spectrum (solid line) of dimers, along with the associated extinction spectra (dash line) measured using standard extinction spectroscopy. The resonance corresponds to the plasmonic dipolar mode of the structure. e) Same as image d for 4 sizes of dimers corresponding to the 3λ/2 plasmonic mode. f) Same as image d for 2 sizes of dimers corresponding to the 5λ/2 mode.

FIG. 3

 3 FIG. 3. a) SEM image of a gold dimer (L = 320 nm) made by e-beam lithography. b) Measured distribution of the thermally-induced OPD when heating a single dimer with a uniform circular beam at λ0 = 740 nm. c) Radial profile of the OPD distribution in image b in good agreement with the fit function (Eq. 2). d) NP absorption spectra (solid lines) measured on a single dimer and averaged on an array of dimers along with the measured extinction spectrum (dash line).

  FIG. 4. a) Absorption spectra of various liquids. b) Radial profile of the OPD distribution around a single gold NP at λ0 = 720nm and λ0 = 960nm (D = 6.2 µm et h = 7 µm). In the second case, the heating of the liquid is responsible of an abnormal linear slope at large distances.

S. Berciaud, L. Cognet, and B. Lounis, Nano Lett. 5, 2160 (2005).