Raphaël Charrondière
email: raphael.charrondiere@inria.fr

Florence Bertails-Descoubes

Sébastien Neukirch

Victor Romero

Numerical modeling of inextensible elastic ribbons with curvature-based elements

Keywords: Slender elastic structure, developable surface, curvature-based elements, constraint-based nonlinear optimization MSC: B, M, N, K

come

Numerical modeling of inextensible elastic ribbons with curvature-based elements

. Introduction

From cable plies and fashion accessories to hair ringlets, flexible band-shaped structures, namely ribbons, are widely present in our daily environment. While the mechanical study of ribbons was initiated in the s [ ], there is recently a renewed interest in understanding and modeling such fascinating structures [ ]. Formally, a ribbon-like structure lies between rod and plate models: its thickness remains negligible in comparison to its width, the latter being itself negligible in comparison to its length. Under clamped-free boundary conditions, it is reasonable to consider the ribbon as a developable surface [ , Sec. . . ], that is, a surface that can be flattened onto a plane without any local distortion. What is more, this surface remains inextensible, meaning that it cannot be stretched nor sheared, but it deforms isometrically, by pure bending (in contrast to situations where wrinkles appear, see e.g. [ ]). We assume that this surface is always smooth, that is, at least twice differentiable. In this paper, we further limit ourselves to the case of a rectangular ribbon geometry with a fixed thickness, composed of an isotropic and linear elastic material, see Fig. .  Developability implies that whatever the state of the ribbon, its surface is ruled, i.e. it is the union of infinitely thin segments called rulings. Note that these rulings do not correspond to material lines of the ribbon, but only to geometric lines which continuously evolve during deformation. The presence of rulings is an important feature of the ribbon surface and was discussed in connection with the shape of the Möbius strip soon after its introduction, see e.g. [ , , ]. Computing the equilibrium shape of the Möbius strip is partly responsible for the recent interest in elastic ribbon equations, see for example [ , ] where equilibrium equations are written in their strong formulation and the boundary value problem is solved using the Auto package [ ]. Other recent works include [ ] who show that ribbon equations can be written in a form similar to the Kirchhoff equations for elastic rods, albeit yielding different behaviors under the same loading [ , , ]. While the number of numerical models for elastic rods has exploded from the s to recent years, see e.g. [ , , , ,

, to our knowledge very few computer codes have been designed specifically for ribbons. An exception is the dynamic ribbon simulator proposed in [ , ], which relies on a nodal discretization of the ribbon surface. This model requires a number of elements that increases as the surface is getting curved, making its use inefficient in the case of strongly curled ribbons. In addition, inextensibility of the ribbon is taken into account by means of additional constraints, which are non-linear. Our approach relies on a radically different discretization scheme, where curvatures (instead of positions) are chosen as degrees of freedom. Compared to inextensible rod models, only two (instead of three) curvatures are sufficient to describe the actual degrees of freedom of a ribbon, due to the developability constraint. Our final choice for a reduced set of curvatures allows us to capture precisely the geometry of the ribbon for arbitrary natural curvatures, while guaranteeing the developability and inextensibility of the ribbon in an exact and automated manner.

Contribution. Our approach is based on the theoretical Wunderlich ribbon model presented in [ ], which we briefly recall in section . We numerically compute stable equilibrium solutions of ribbons, drawing inspiration from the curvature-based models developed for inextensible elastic rods [ , ].

In particular, we choose to use material curvatures of the ribbon as main variables of our system, making it easier to recover a smooth and inextensible ribbon centerline, regardless of the chosen resolution -which contrasts with classical nodal models. The developability constraint is simply accounted for by reducing the number of variables to two scalar fields: one material curvature, and one variable accounting for the direction of the rulings. A difficulty, however, is to reconstruct the geometry of the ribbon accurately from this reduced set of degrees of freedom. Inspired by the super-clothoid rod model [ ], we use power series for the integration of the ribbon centerline from curvatures, while keeping a fine control on the numerical error (section ). Thanks to the ruling direction as a primary variable, the developable surface of the ribbon can then be constructed from the centerline with excellent accuracy. A last difficulty consists in ensuring the non-crossing of the rulings during the deformation of the ribbon, so as to guarantee that the ribbon surface remains smooth. We treat this condition as an inequality constraint in our energy minimization problem. We analytically compute the first and second derivatives of our discrete ribbon energy and use an interior point minimization method relying on the energy Hessian [ ]. Numerical solutions are then found in a few seconds on a standard machine (see inset figure), and this for a wide range of parameters (section ). Finally, we show that our code compares favorably to standard shooting and collocation methods, and to experiments. Furthermore, our results highlight significant differences compared to equilibria generated by thin elastic rod models, as well as a substantial speed-up compared to a more general shell numerical model, which confirms the benefit of designing a numerical model dedicated to ribbons (section ).

. Wunderlich ribbon parameterized by curvatures

The ribbon model of Sadowsky [ , ] defines a one dimensional elastic energy derived for the case h w L, with h the thickness, w the width, and L the length of the ribbon. It enforces a developable mid-surface, as opposed to elastic rods models. Our approach relies on the Wunderlich ribbon model [ ], as described in [ ], which is a generalization of the Sadowsky model when w grows towards L. In the following we show how to exploit this D ribbon energy, without using existing auto-differentiation tools, to reach a desired and controlled precision at a minimal cost.

Notation. Scalars will be denoted by plain small letters, vectors of R 3 by bold small letters, and matrices of the rotation group SO 3 by slant capital letters. For instance, we shall write v ∈ R, a ∈ R 3 , and R ∈ SO 3 . Scalar product between two vectors u and v of R 3 will be denoted by u • v. We further define [•] × : R 3 → A 3 as the linear operator transforming the vector product into a matrix multiplication u × v = [u] × v, with A 3 the set of real asymmetric matrices of size by . Finally, x (s) = dx ds will denote the spatial derivative of the function x(s).

. . Base geometric parameterization

We consider a rectangular ribbon with constant thickness h. Let r : [0, L] → R 3 be the centerline of the ribbon, parameterized by the arc length s ∈ [0, L], with L the total length of the ribbon. The material (Cosserat) frame R : [0, L] → SO 3 describes the surface locally around the centerline. For each point r(s), we denote by d 1 (s), d 2 (s), and d 3 (s) the column vectors of the rotation matrix R(s). Geometrically, the component d 3 (s) stands for the unitary vector tangent to the centerline at s, d 2 (s) for the unitary surface normal, while Note that q(s) is not normed in general, its norm being 1 + η 2 (s). From the generatrix, the ruling then reads L(s) = {r(s) + v q(s); v ∈ S(s)}, where the interval S(s) ∈ R is specified in the following. Finally, the surface of the ribbon, representing the set of all rulings, is given by the function

d 1 (s) = d 2 (s) × d 3 (s) is imposed by the SO 3 structure.
φ : [0, L] × S → R 3 (s, v) → r(s) + v [d 1 (s) + η(s)d 3 (s)] q(s) . ( )
In the present case of a rectangular ribbon, the set S(s) boils down to a constant interval. Indeed, consider a deformed ribbon with its ruling

L(s) = {r(s) + v [d 1 (s) + η(s)d 3 (s)]; v ∈ S(s)}.
The interval S(s) is actually easier to evaluate on a flat configuration, so we deform the ribbon isometrically until it becomes flat, and we denote variables in the flat configuration with bars. Here consider the deformation of the initial ruling L(s), which reads L(s) = {r(s) + v [ d1 (s) + η(s) d3 (s)]; v ∈ S(s)}. Note that in this expression, only the material frame has changed compared to the deformed configuration. All distances and angles are preserved, so the ruling orientation η(s) and v are kept the same, and as a result the interval S(s) is invariant through the isometric flattening operation. We show in figure -left, the ruling on the flat configuration. At the edge point r+ (s) = φ(s, w 2 ) we have v = w 2 , so we deduce S = [-w 2 , w 2 ], which is independent of s.

The reader can remark that rulings are not material lines: when the ribbon deforms, the set of material points forming the rulings changes. Moreover, L(s) denotes a valid set of rules for the flat configuration, and they might not be parallel. Wunderlich (or smoothness) constraint. As illustrated in figure -right, the rulings of a smooth developable surface cannot intersect inside the surface. To be convinced, the reader can make a simple experiment: take a paper sheet, bend the edges to form a cone, and try to get the top of the cone (singularity point) inside the sheet; unless you crease or rip the sheet, this is impossible. This non-intersection constraint can be expressed globally on the flat configuration as

d3 d1 w 2 w 2 r+ (s) r-(s) 1 η(s) q ( s ) v vη(s) φ(s, v) r(s) s r(s a ) r(s b ) I w 2 w 2 d3 d1
∀s a , s b ∈ [0, L], s a = s b , ∀v a , v b ∈ S, φ(s a , v a ) = φ(s b , v b ),
or, say otherwise, the point I lying at the intersection of two rulings passing through s a and s b , if it exists (non-parallel rulings), should lie outside the ribbon surface (see figure -right). By noticing that in the flat configuration, the material frame { d1 , d2 , d3 } is uniform, and using expression ( ), one realizes that intersection can only happen if v a = v b . Furthermore, the global constraint expressed above is equivalent to a local constraint where s a and s b are infinitely close, that is, for s b = s a + ds. Expanding φ(s a + ds, v a ) at the first order, the intersection condition reads ∂ φ ∂s (s a , v a ) = 0, which yields 1 + v a η (s) = 0. The condition for point I to lie outside the ribbon, |v a | > w/2, then reads

∀s ∈ [0, L] |η (s)| < 2 w , ( )
which also encompasses the case when rulings are parallel (η (s) = 0).

. . Curvature-based parameterization

Like for curvature-based rod models [ , ], our idea is to use the material curvatures of the ribbon as main variables of our ribbon model, instead of taking the base geometric variables (r, R, η). This choice allows us to better stick to the actual degrees of freedom of the ribbon, and in particular, to capture the inextensibility of the centerline intrinsically. Like for rods, we first need to relate material curvatures to the geometry (r, R, η) of the ribbon. In contrast to rods however, we also need to reduce the set of curvatures in order to account for the inextensibility and developability of the ribbon surface.

From curvatures to the centerline. Inextensibility of the centerline reads r (s) = d 3 (s), which integrates to r(s) = r(0)+ s t=0 d 3 (t)dt. Besides, the SO 3 structure implies the existence of a vector Ω(s) ∈ R 3 such that R (s) = [Ω(s)] × R. The vector Ω(s) is called the Darboux vector, and reads Ω(s) = κ 1 (s)d 1 (s) + κ 2 (s)d 2 (s)+κ 3 (s)d 3 (s) where κ 1 (s), κ 2 (s), and κ 3 (s) are respectively the normal curvature, the geodesic curvature and the geodesic torsion of the ribbon at arc length s -here we simply refer to these quantities as material curvatures, or just curvatures. It is convenient to manipulate curvatures directly by transforming the relation into R (s) = R(s) [κ(s)] × , where we have introduced the curvature vector κ(s) = (κ 1 (s), κ 2 (s), κ 3 (s))

T , such that Ω(s) = R(s)κ(s). Once the curvature vector field κ(s) and the initial conditions (r(0), R(0)) = (r 0 , R 0 ) are given, we obtain a Cauchy integration problem for (R, r), called the Darboux problem,

P =    (r(0), R(0)) = (r 0 , R 0 ) R (s) = R(s) [κ(s)] × r (s) = R(s) e z .
( ) Here e z = (0, 0, 1) T stands for the third vector of the canonical basis, allowing us to select the third vector of the material frame, as d 3 (s) = R(s) e z . The Darboux problem ( ) possesses a unique solution. However it is a linear nonautonomous differential system and apart from specific cases, the solution has no closed-form expression in the general case , i.e. it cannot be represented with a finite number of elementary operations. Nevertheless, when the curvature is a polynomial function of s, it can be shown that the solution is analytic , that is, it can be expressed as a convergent power series [ ]. In section . , we show how to adapt the power series integration method developed in the case of rods [ ] to the present case of ribbons.

Isometry of the surface. Geometrically, a major difference between rods and ribbons is that ribbons are characterized by additional constraints: their surface should deform isometrically from a planar configuration, and thus remain

In particular, due to the non-commutativity of SO(3), taking the exponential of the primitive of the curvature matrix does not work here as it works for autonomous linear systems of ODEs, see e.g. [ , , ]. The term analytic is more inclusive compared to a closed-form expression as it encompasses infinite sets of elementary operations. In this paper, we further misuse the term analytic to characterize a computational method which relies on a truncated power series where the sum of truncated terms is numerically negligible: machine precision is reached while omitting these terms.

inextensible and developable. In the following we simply provide a synthesis of these conditions, and refer the reader to [ ] for mathematical justifications.

On the one hand, inextensibility of the ribbon surface imposes the invariance of the geodesic curvature. Since κ 2 = d 3 • d 1 , this implies κ 2 = d 3 • d1 , which in the case of a rectangular ribbon boils down to

κ 2 = 0. ( )
On the other hand, developability implies a zero Gaussian curvature for the surface, which can be written as

κ 3 = ηκ 1 . ( )
Reduced degrees of freedom. Given the additional constraints above, we propose in our model to use η and κ 1 as degrees of freedom. This choice allows us, via equations ( ) and ( ), to compute κ(s) = (κ 1 , 0, ηκ 1 ) T for all values of η and κ 1 . Note that choosing κ 1 and κ 3 instead would cause numerical difficulties for computing η as κ 1 vanishes (likewise for the choice of η and κ 3 ).

. . Energies

We consider a naturally curved elastic ribbon clamped under gravity. We do not take into account contact of the ribbon with another object, nor self-contact.

Elastic energy. According to [ , Sec. . ], the Wunderlich elastic energy of a rectangular ribbon with natural normal curvature κ N reads

E W = Dw 2 L 0 κ 2 1 (1 + η 2 ) 2 1 wη ln 1 + η w 2 1 -η w 2 -2κ N κ 1 (1 + νη 2 ) ds, ( )
where ν is the Poisson ratio, Y the Young modulus, and D = Y h 3 12(1-ν 2 ) the bending modulus of the surface. In the limit η w → 0, the Wunderlich energy E W boils down to the so-called Sadowsky energy,

E S = Dw 2 L 0 κ 2 1 (1 + η 2 ) 2 -2κ N κ 1 (1 + νη 2 )ds. ( )
Gravitational energy. We define the potential gravitational energy in a standard way as

E G = - L 0 w 2 -w 2 ρh φ(s, v) • g da(s, v), ( )
where ρ is the material density, h the thickness of the ribbon, g the gravitational acceleration, and da(s, v) = (1 + v η (s))dvds the area of the infinitesimal (trapezoidal) integration element. Integrating ( ) along v yields

E G = - L 0 ρhw r(s) + η (s)w 2 12 q(s) • g ds. ( )
The length of the ribbon being much greater than its width, we consider in the following that the term η (s)w 2 12 q(s) is negligible, and omit it in our computations for the sake of simplicity (see section . ). In practice, we have not observed any discrepancy related to this approximation.

Remarks. Our goal in the following is to compute stable configurations of the ribbon, that is, local minima of its total mechanical energy. Before discretizing our set of equations as done in next section, we can already anticipate some difficulties. First, the non-convexity of the energy appears as a major pitfall for standard optimization techniques (such as Newton). Another challenge arises from the logarithmic term in the Wunderlich energy, which diverges when the constraint ( ) becomes active. On the one hand, we shall take this constraint into account while minimizing energy, and on the other hand we shall build a limit case for the Wunderlich energy when the constraint is active. In sections . and . , we present a solving method able to handle all these issues in a robust and efficient way.

. Discretization

To tackle our energy minimization problem numerically, we adopt a Galerkinlike discretization strategy. The ribbon is cut in p (trapezoidal and curved) elements, each element being described by a linear normal curvature and a quadratic geodesic torsion. We first explain how this discrete scheme, combined with a robust power series summation algorithm, leads to a consistent and accurate approximation of the ribbon kinematics. Then we show how it can be extended to compute accurately the mechanical energy as well as its first and second derivatives.

. . Discrete kinematics

We cut the centerline in p segments S i = {r(s), s ∈ [s i , s i+1 ]} of length s i+1 -s i = i , with i ∈ {0, ..., p -1}. In the transverse direction, the two boundaries of each ribbon element E i are determined by the rulings passing through r(s i ) and r(s i+1 ), respectively. Each element E i thus takes the form of a trapezoid on the flat configuration (see figure -left).

Linear and quadratic form functions. Our goal is to approximate our primary variables η and κ 1 by basis functions parameterized by a finite number of degrees of freedom. The simplest choice, consisting in taking piecewise uniform functions as often done in the case of rod models [ , , , ], is however inappropriate here as it leads to a discontinuous ribbon surface. Indeed, if the function η is discontinuous at element boundaries, then elements cannot even connect to each other properly (see figure -right). In order to guarantee a continuous and sufficiently smooth ribbon surface, we choose η and κ 1 as continuous and piecewise linear functions of s, i.e. ∀s ∈

[s i , s i+1 ], κ 1 (s) = a i s + b i and η(s) = c i s + d i , with a i , b i , c i , d i ∈ R. As a result, the geodesic torsion κ 3 = η κ 1 takes the form of a piecewise quadratic function of s, i.e. ∀s ∈ [s i , s i+1 ], κ 3 (s) = a i c i s 2 + a i d i + b i c i s + b i d i . Finally, the curvature vector κ is a polynomial of degree 2 in s on each element, ∀s ∈ [s i , s i+1 ], κ(s) = λ i 0 + λ i 1 s + λ i 2 s 2 ( ) with λ i 0 = (b i , 0, b i d i ) T , λ i 1 = (a i , 0, a i d i + b i c i ) T
, and λ i 2 = (0, 0, a i c i ) T . Taking a polynomial curvature vector of degree ≥ 1 per element however implies that the Darboux problem ( ) has no closed-form solution on the element (see section . ). The recourse to a numerical integration method therefore proves necessary, with the risk of turning this crucial kinematical step, which needs to be accurately resolved, into an overwhelming step. Fortunately, an analytic integration method, based on power series, has been deployed in [ ]. We show that this method can be adapted to solve our ribbon kinematics in a both accurate and efficient way.

Numerical integration with power series. For the sake of clarity, in this section we consider a ribbon with a single element, parameterized by its arc length s ∈ [0, L], and keep previous notation while omitting the subscript i. The extension to a chain of several elements is discussed in section . .

The curvature vector κ being a polynomial, it is a convergent power series. Then, it can be shown that r and R are convergent power series too [ , Theorem ]. Expanding r(s) and R(s) as power series r(s)

= +∞ n=0 r n s n , R(s) = +∞ n=0 R n s n , we rewrite the Darboux problem ( ) as    r(0) = r 0 , R(0) = R 0 +∞ n=1 nR n s n-1 = +∞ n=0 R n s n [κ(s)] × +∞ n=1 nr n s n-1 = +∞ n=0 R n s n e z .
( ) By uniqueness of the coefficients of a power series, and using the linearity of the operator [•] × , solving problem ( ) amounts to solving the following recurrences,

                 R 0 = R 0 R 1 = R 0 [λ 0 ] × R 2 = 1 2 R 1 [λ 0 ] × + R 0 [λ 1 ] × ∀n 3, R n = 1 n 2 k=0 R n-1-k [λ k ] × r 0 = r 0 ∀n > 0, r n = 1 n R n-1 e z .
( )

One main difference with the rod model derived in [ ] is the three-term (instead of a two-term) recurrence for the material frame R. This order increase is due to our quadratic geodesic torsion, which itself stems from the Wunderlich constraint ( ), inherent to ribbons. Such a difference justifies the need for revisiting the algorithm proposed in [ ]. But still, both problems share many common features. First, due to the structure of the Darboux problem, the ribbon series ( ), like the rod series, are characterized by a fast convergence of their remainder, and this, for a relatively low rank n [ ]. In practice, it is therefore sufficient to consider less than a hundred terms for evaluating the series accurately. Second, if we naively sum the different terms of these series using floating-point arithmetic, we are quickly faced with a drastic precision loss. This problem is known as catastrophic cancellation, and has been carefully analyzed in [ ]. For the sake of completeness, we recall the base lines of this issue in the following.

Catastrophic cancellation. Problems of precision loss appear when summing terms with very different orders of magnitude and opposite signs. In modern computers, a real number is represented by a mantissa M , an exponent n and a sign s = ±1. The value of the real is then s × 1, M × 2 n where the comma is the same as the one that is used to write the decimal part of a number (except that it is written in binary format here). In double precision, the mantissa M is of size m = 52 bits, allowing to represent numbers not larger than 2 52 ≈ 10 16 . Consequently, evaluating 1 + 10 16 numerically yields 10 16 , a phenomenon called absorption. The problem is that if we now subtract 10 16 from the initial sum, we find 0 instead of 1, which corresponds to a precision loss of 100 %! The series encountered in ( ) includes terms of alternating signs (due to skew matrices), whose absolute value quickly reaches a very large value as the rank n increases, before becoming very small beyond a certain rank. As a result, all the ingredients are present for catastrophic cancellation to arise when summing the terms of the series.

Error control through element subdivision.

In [ ], it is explained how to avoid catastrophic cancellation when integrating the Darboux problem using power series. A key observation is that numerical issues only occur when s is greater than a certain threshold. Hence the idea of subdividing the element into subelements, and to solve the Cauchy problem iteratively on each sub-element, starting from the end of the previous sub-element (shifting first point of the current sub-element to s = 0). If each sub-element is small enough, integration is guaranteed to be performed safely, i.e. without precision loss. The question is how to estimate a suitable size for each sub-element. One idea is to derive a theoretical upper-bound s max for the length of a given sub-element, such that the terms of the series to be summed do not exceed in norm the maximal range covered by floating-point arithmetic. Said otherwise, if one wants a precision of = 2 -m 2 (where m is the number of bits for expressing the mantissa), one should make sure that the largest value of R n s n is bounded by 2 m 2 . In his PhD thesis, Casati has derived an upper-bound for R n s n which is valid for any degree of the curvature polynomials [ ], 2 otherwise. Consequently, the more the ribbon is curled, the more sub-elements are needed to preserve accuracy.

R n s n ≤ e C(s) with C(s) = 2|s|

. . Discrete mechanical energy

We now have to compute the total potential energy corresponding to our discrete ribbon kinematics. As before, we focus here on a ribbon made of a single element, before generalizing to a chain of elements in section . . Recall that on each element of the ribbon, η and κ 1 are polynomials of degree 1 in s, κ 1 (s) = as + b and η(s) = cs + d.

Discrete Wunderlich energy. By deriving the expressions given in section , the elastic Wunderlich energy ( ) becomes

E W = Dw 2 [T (c)P (a, b, c, d) -Q(a, b, c, d)] ( a)
with

T (c) = 1 cw ln 1 + cw 2 1 -cw 2 , ( b) P (a, b, c, d) = 0 (as + b) 2 1 + (cs + d) 2 2 ds, ( c) Q(a, b, c, d) = 0 2κ N (as + b) 1 + (cs + d) 2 ds. ( d)
The only difficulty here is the computation of T which is not defined if c = ± 2 w . Note that in theory, the smoothness constraint ( ) should prevent us from reaching these bounds. In practice, we have however observed that when getting close to the bounds, numerical problems arise. In such cases, we instead use Taylor's expansion formula T (c) = 1 + 1 3 cw 2 2 + o(c 3 ). We switch to this formula as soon as cw 2 < T , where T is an empirical threshold value. In practice, we have found that choosing T = 10 -5 leads to satisfactory results in double precision.

Discrete gravitational energy.

As mentioned earlier, the length of the ribbon being much greater than its width, we approximate the potential energy of gravity ( ) by concentrating the mass on the centerline, E G -L 0 wρh r(s) • g ds. This integral can be computed from the power series

E G = - +∞ k=0 E k s k • g, with    E 0 = 0 E 1 = wr 0 ∀k > 0, E k+1 = 1 k+1 wr k . ( )

. . Composition of the pieces

We now consider the general case of a chain of p elements. Recall that for accuracy purpose, each element itself is subdivided in sub-elements. In previous section we have defined, for each sub-element, four ready-to-use functions a, b, c, d). The parameters R 0 , r 0 , b, d all depend on previous sub-elements. This dependency complexifies the derivation of the energy, as the functions are composed as many times as there are sub-elements.

R(R 0 , a, b, c, d), r(R 0 , r 0 , a, b, c, d), E G (R 0 , r 0 , a, b, c, d) and E W (
To connect the elements into a single ribbon, we either need to express each element as a function of the previous elements, or we can connect each element to its previous element using a nonlinear constraint. To avoid any accuracy loss due to the linearization of such a constraint, we choose the first option (meaning that ( ) is the only constraint in our problem).

To retrieve a smooth surface at connections between two elements, it is sufficient to enforce curvature continuity. For the sake of simplicity, we express κ 1 and η on each element E i as functions of a local arc length

s i ∈ [0, i ]. We thus have κ 1 (s i ) = a i s i +b i loc and η(s i ) = c i s i +d i loc , where b i loc = b i + i-1 n=0 a n n and d i loc = d i + i-1 n=0 c n n . Continuity conditions thus read ∀i > 0, b i loc = a i-1 i + b i-1 loc ( a) d i loc = c i-1 i + d i-1 loc , ( b) hence the full system is parameterized by (b 0 , d 0 , a 0 , c 0 , a 1 , c 1 , • • • , a p-1 , c p-1 ).
When the energy is differentiated, like for sub-elements, composition between elements should be taken into account, leading to some numerical complexity: typically, the computation of the Hessian of the system is cubic in the number of elements. Fortunately, the number of elements remains low (around 20). Because of the nature of the system, the numerical implementation is long and error-prone. During implementation, we have carefully checked analytic differentiation using comparisons with finite differences.

Despite these challenges, we show in the following that computing the energy, its gradient and Hessian analytically was worth it: this proved to be crucial in the success of our overall method for computing the statics of a ribbon.

. Computating stable equilibria of ribbons . . Robust energy minimization

Because of non-convexity and variable bounds, our local minimization problem is particularly challenging. We have tried many minimization algorithms before being able to solve it in a robust way. We eventually rely on the opensource library Ipopt [ ], an interior-point method which, in addition to regularization and line-searching, handles equality and inequality constraints. We feed this algorithm with both the gradient and the Hessian of the ribbon energy, calculated in an analytic way. This algorithm proves to work well in a large number of scenarios and outperforms other minimization strategies when considering the following criteria: it never diverges, and is the fastest (table ).

To better understand the role played by the algorithm with respect to our problem, we have conducted a number of comparison experiments with various minimization methods, as reported below. These comparisons allowed us to draw some conclusions on the key ingredients that are necessary for a successful computation of equilibria of our ribbon model,

• Use a second-order approach. We have always observed poor convergence for methods that do not exploit the information of the exact Hessian;

• Use regularization to deal with a non-positive Hessian robustly;

• Use line-searching to improve convergence;

• Finally, perform minimization on a dimensionless problem. For secondorder minimization strategies, we have observed a speed gain between ×3 and ×10, meaning that non-dimensionalization greatly improves the conditioning of the linear system to be solved.

In the following we present and compare various standard minimization methods to justify our final selection of Ipopt.

. . Evaluation of various local minimization methods

Many algorithms exist to find a local minimum of a function. We have compared our final choice, named Ipopt-Hess, against five standard methods that we remind to the reader below. All these methods correspond to an instance of the generic minimization algorithm provided in Algorithm . For a more comprehensive survey and analysis of minimization algorithms, we refer the reader to [ ]. Here, we assume that the energy is twice-differentiable. At each step a direction d is computed and the current point x is moved along this direction. The function stepSize may for example be a line-search that regulates the amplitude of this movement in order to respect variables bounds and/or energy decrease. The convergence criterion used is often based on the norm |∇E(x)| of the gradient. When this norm is below a given tolerance, a critical (i.e. minimum, saddle, or a maximum) point has probably been reached.

Input

First-order methods. First-order methods only use the value of the energy and its gradient to perform the minimization. The standard gradient descent algorithm, called GradDesc, uses a first-order approximation of the energy and follows the opposite of the gradient to decrease the energy: sysMat is the identity, and we take stepSize = 1. After a few iterations, first-order methods are however generally slow as the gradient becomes small when approaching a minimum.

Second-order methods. To increase convergence speed, a second-order approximation of the energy can be used. In this case, the function sysMat computes the Hessian HE(x) of the energy. If we still take stepSize = 1, we have a pure Newton method, called Newton. Near the solution, the error |∇E(x)| is small and the method benefits from a quadratic convergence. However, this method only converges to an extremum, not necessarily a minimum. It may also diverge for a number of reasons: (i) convergence is guaranteed only for convex functions, which is not the case here, (ii) the second-order approximation is local, if d is too large, the new point x can be unreasonable, i.e. with too large curvature variables, and (iii) the term ( b) in the energy can diverge if the constraint |n| ≤ 2 w , see ( ), is not taken into account. Since the pure Newton method does not necessarily provide a descent direction in the non-convex case, we introduce a variant called NewtonC. When d is not a descent direction, i.e. d T HE(x) d < 0, NewtonC switches to a gradient descent step. Moreover, a line-search is used and the bound check |n| ≤ 2 w is turned on. Convergence is improved in some cases, but the algorithm often remains stuck in a gradient descent region, where convergence speed is not satisfactory.

Finally, NewtonR is yet another variant introduced to cope with nonconvexity. There, sysMat computes the Hessian of the energy and its smallest eigenvalue µ. If µ < 10 -8 , the Hessian is regularized by adding 10|µ|I to it. A line-search complements this method which is virtually as robust and efficient as Ipopt-Hess.

As a side note, we have also tested NetwonCG, a truncated Newton method using a modified version of the Conjugate Gradient algorithm. While performing better than NewtonC, this method proved to be slower than NewtonR, presumably due to the fact that our Hessian is non-sparse.

Quasi-Newton methods. As explained in section . , the cost of calculating the gradient increases quadratically with the number of elements, while the cost for the Hessian increases cubically. Quasi-Newton methods avoid the computation of the Hessian of the energy. A well-known example is the BFGS algorithm [ ] which builds an approximation of the Hessian, and corrects it at each iteration. Under a convexity assumption, the approximation is known to converge to the real Hessian. This method nevertheless requires a larger number of steps which means a larger number of evaluations of the energy and its gradient, and we show in table . that the L-BFGS method from Ipopt is unsatisfactory.

Evaluation. We compare our different numerical strategies to find a target stable equilibrium configuration. For the sake of simplicity we use here the simpler Sadowsky model and do not take the smoothness constraint ( ) into account (the influence of this constraint is later discussed in section . ). We normalize the physical quantities with the length L of the ribbon and its bending rigidity Dw. Consequently, once the Poisson ratio is chosen (we use ν = 0.4) there are only two dimensionless parameters, Γ = M gL 2 /(Dw) and K N = κ N L, where M = ρhwL is the total mass of the ribbon and g = 9.81 m/s 2 is the acceleration of gravity. Here we only consider a uniform natural curvature dκ N /ds = 0, ∀s. As illustrated in figure , we study three cases of increasing complexity, from a gently curved, sagging ribbon (K N = -10, Γ = 100) to a highly curved, heavy ribbon (K N = -100, Γ = 7500). The structure is clamped at one end, and free at the other. We choose the orientation g = -e Table : To test the added value of having access to the Hessian of the energy, we compute the static equilibrium of a ribbon with Ipopt in two ways: (first line) feeding Ipopt with the Hessian HE, or (second line) not providing the Hessian but using the L-BFGS option of Ipopt. We use a ribbon with length L = 50 cm, width w = 0.5 cm, thickness h = 0.1 mm, area density ρ = 2000 kg/m 3 , Young modulus Y = 1.8 GPa, Poisson ratio ν = 0.5, and natural normal curvature κ N = 80 m -1 . The ribbon is clamped with a α = π 4 angle with respect to the vertical axis.

and d 3 (s = 0) = e z , with α = π/4. In the minimization routines, the two parameters K N and Γ are directly set to their targeted values and the routines are initiated with the natural shape of the ribbon. For the first case we use p = 10 elements, while for the two others we use p = 20.

Results are listed in Table and we see that methods based on pure Newton, gradient descent steps (or a mix of the two) are unable to converge to the equilibrium starting from the natural shape. By contrast, methods monitoring the Hessian matrix and enforcing its positiveness converge in all three cases. Line-search procedures, present in NewtonR and Ipopt-Hess for example, sometimes reduce the step size taken during the Newton iteration, thereby preventing the algorithm from sending the actual point too far in the energy landscape. Table : Reaching a stable equilibrium with various algorithms, on the three test cases depicted in figure , using the Sadowsky energy. Among minimization methods applied on our discrete model (top), only NewtonR and Ipopt-Hess perform robustly: they both combine exact second-order information and regularization. Among continuation methods using a strong formulation of ribbons (bottom), Auto proves to be very fast as it even outperforms best minimization methods applied to our model for a large number of elements. Current continuation methods are however limited to the Sadowsky energy, and cannot easily accommodate the Wunderlich energy nor the smoothness constraint ( ), which in contrast are well accounted for by our approach.

Finally, in order to test whether the numerical difficulties of the first three minimization methods arise from the non-quadratic elastic energy of the ribbon model, we have performed the same equilibrium search with the Kirchhoff model for elastic rods discretized with linear-curvature elements [ ]. Overall convergence is somewhat faster and the pure Newton method manages to converge in this case, but no substantial difference is observed.

. Results and validation

We have compared our numerical model against various alternative approaches. These comparisons show the accuracy, robustness, practicability and versatility of our model, which allows us to envision interesting applications such as the fast exploration of multiple ribbon static configurations (please refer to our accompanying video for a better visualization of some results). We additionally show recent experimental study conducted with the aim of confronting the model to real scenarios.

. . Comparison with the Kirchhoff model for thin elastic rods

We consider a naturally curved ribbon hanging under its own weight. The ribbon is clamped at one extremity and free at the other. We use the following physical parameters: length L in the range [1, 50] cm, width w = 1 cm, thickness h = 0.1 mm, density ρ = 1000 kg/m 3 , Young's modulus Y = 1.8 GPa, Poisson's ratio ν = 0.5, acceleration of gravity g = 9.81 m/s 2 . The natural curvature is chosen to be κ N = 60 m -1 , hence in its natural state the ribbon is winding on itself, making approximately 5 turns. We numerically compute the equilibrium state of the ribbon using our approach with the Wunderlich energy. The structure is divided into 20 elements, and we start the minimizing procedure with a flat and horizontal initial state. For a given length L, convergence is reached after 36 iterations in 2 seconds. Equilibrium states are shown in blue for various lengths in figure -right.

Next we use the inextensible and unshearable Kirchhoff model [ , ] for twisted elastic rods to compute the equilibrium of the system. The Kirchhoff model is accurate for rods with, for example, circular or rectangle cross-sections as long as the aspect ratio of the section is not too large. For rectangle sections with w/h 1, the ribbon and the Kirchhoff models generate static configurations that can be different from each other, with experiments validating the ribbon model [ , ]. The major difference between the two models comes from the twist-curvature coupling present in the case of the ribbon. In the Kirchhoff rod model, the elastic energy is a quadratic function of the curvatures

E K (κ 1 , κ 2 , κ 3 ) = L 0 1 2 K 11 (κ 1 -κ N ) 2 + K 22 κ 2 2 + K 33 κ 2 3 ds ( )
where the stiffness coefficients are, for a rectangular section [ ]

K 11 = Y I 1 , K 22 = Y I 2 , K 33 = Y 2(1 + ν) J ( a) I = h 3 w 12 , I 2 = hw 3 12 , J = c J h 3 w 3 ( b) c J = 1 - 192 π 5 h w ∞ k=0 tanh (2k+1)πw 2h (2k + 1) 5 . ( c)
In the present case where w = 100h, we have J h 3 w 3 . We choose to recalibrate the rod bending stiffness Y I 1 with its equivalent for ribbons Dw by taking

K 11 = Y I1
1-ν 2 = Dw. Both models then have the same rigidity compared to bending in the soft direction. Using physical parameters of the structure, we evaluate K 11 = 0.002, K 22 = 15, and K 33 = 0.006 mN m 2 . For the Kirchhoff rod model, we use the same spatial discretization (i.e. super-clothoids [ ]) as above and also perform the approximation of putting the mass on the center-line. The energy of gravity is then the same in both models, only the elastic energies differ. Using the same minimization procedure, we obtain the equilibrium solutions shown in orange in figure -right.

For short lengths, the elastic energy dominates and the structure adopts a shape close to its natural curvature. Then, from L = 5 cm, gravity is no longer negligible and the structure unfolds, but remains in the plane. Finally from L ≈ 29 cm the ribbon swings out of the plane, and torsion appears. We notice that the two models perfectly coïncide at the beginning and deviate from each other as soon as torsion comes into play.

One could still wonder whether a more general rod model, which allows for stretching and shearing, could better approximate the behavior of a ribbon. In Appendix B, we have performed some simulations based on the Antman rod model [ ], and we show that such a model almost exactly behaves like the Kirchhoff rod model: the shearing/stretching modes are of no help for modeling the correct behavior of a ribbon, and we observe a major discrepancy between our ribbon simulations and rod simulations, whatever the rod model used -Kirchhoff or Antman. Intuitively, this is due to the fact that none of these models guarantees the developability of the ribbon surface: as shown in [ ], a ribbon is actually more constrained than the Kirchhoff rod model: relaxing the Kirchhoff rod model in the hope of capturing ribbon deformations is thus vain. .

. Comparison with the Naghdi thin elastic shell model

Next we study how finite element computations of elastic shells compare with the elastic ribbon model and with our numerical scheme. We use the first case of section . (K N = 10, Γ = 100, and ν = 0.4, Sadowsky's formulation, see figure -left) and compare it to the equilibrium solution of the Naghdi shell model [ ] -a general thin shell theory relying upon a mid-surface combined with straight (but not necessarily normal) transverse directors (so-called Reissner-Mindlin kinematic assumptions) -with w/L = 0.03 and w/h = 10, 20, 50, 100, and 200. Elastic ribbon models, designed for structures with L w h, are derived from shell theories where the shell stretching and shearing energies are discarded and where the curvature energy is integrated over the width w. They have been shown to compare well with experiments when the ratio w/h is large enough, typically w/h > 75, but fail to capture the behavior of the structure when w/h < 30 [ ]. We use the Python-based FEniCS-Shell finite element library to compute the static solution of the Naghdi shell model. This implementation includes director parameterization, in the line of [ ]. The shell is discretized over a triangular 100 × 4 mesh, using the the partial reduced selective integration (PRSI) proposed in [ ] and adapted to nonlinear shells in [ ]. Continuation takes from two (w/h = 10) to four (w/h = 200) minutes to reach {K N = 10, Γ = 100}, to be compared with the 0.1 seconds needed for our ribbon simulator (Table ). In figure , we see that the shell equilibrium differs from the ribbon equilibrium at low w/h ratios, but eventually converges to it at large w/h ratios.

. . Comparison with continuation methods

An alternative way to find equilibrium configurations is to derive differential equations corresponding to the first variation of the energy of the system. In Appendix A we list such differential equations and boundary conditions for the Sadowsky energy ( ). To our knowledge, continuation has never been applied to the Wunderlich energy combined with the inequality constraint ( ) and natural curvature. In this sense, it is less general than our approach.

Continuation methods. The first continuation method we have used is a classical shooting method. Differential equations are integrated with a stiff solver in Mathematica [ ] and the unknown initial values (n(s = 0) and m(s = 0)) are found with the numerical root solver of Mathematica. This shooting method is easy and quick to set up but not robust. The second continuation method we have used is the Fortran-based Auto [ , , ] package. The spatial discretization is based on orthogonal collocation involving Lagrange polynomials of degree N COL and comprising N T ST segments. Once discretized, the boundary value problem takes the form of a nonlinear algebraic system of size N DIM × N COL × N T ST , where N DIM is the size of the differential system (N DIM = 20 in the present case). The algebraic system is solved with a Newton-Chord method using numerical differentiation for the computation of the Jacobian (a.k.a. Hessian in the minimization formalism) matrix.

However, a major difference between energy minimization and continuation methods is that the former searches for stable configurations, while the latter searches for any critical configuration. Hence, when comparing the two, we restrict ourselves to stable equilibria. Cross-validation. We compare the equilibrium shape obtained using our approach with the shape obtained with the shooting method. We use the following parameters: thickness h = 0.1 mm, width w = 1.5 cm, length L = 50 cm, natural curvature κ N = 30 m -1 , density ρ = 1000 kg/m 3 , Poisson's ratio ν = 0.5, and Young's modulus Y = 1.8 GPa. Both shapes, plotted in figure , are identical, which validates our approach on the Sadowsky model. Timing comparisons. We consider again the three test cases studied in section . and illustrated in figure . This time we seek to obtain the stable equilibrium configurations by continuation.

In the continuation routines K N is set directly to its targeted value and Γ is gradually increased from Γ = 0 to the desired target value. For each intermediate Γ value an equilibrium configuration is computed in our shooting routine in Mathematica, and in Auto with a polynomial degree N COL = 4 and N T ST = 10 elements for the first case, N T ST = 15 elements for the second case, and N T ST = 35 elements for the third case. As reported in table , the fastest routine is the continuation method Auto even if the number of steps (i.e. intermediate values of Γ visited) is larger than in the minimization routines. (Note that in addition, a Newton correction is used in Auto at each intermediate Γ value.)

On the Sadowsky problem, Auto thus proves more efficient than best minimization algorithms performed on our discrete energy, especially for a large number of elements. However, continuation approaches cannot easily handle inequality constraints in general. In particular the inequality ( ) of the Wunderlich model is challenging to account for. To the best of our knowledge, only [ ] so far has proposed to tackle this constraint within a continuation method, by adding a regularization term to the Wunderlich energy. However, the authors acknowledge that this ad-hoc solution perturbs the ribbon equilibrium solutions as the width w increases, which makes it only suitable for a moderate range of widths. For large ribbons (large w), the Wunderlich constraint becomes crucial, and we show in the following that removing this constraint from the model may lead to fatal numerical issues. In contrast, our more versatile approach allows us to account for this constraint properly.

. . Active Wunderlich constraint

We now use the Wunderlich model and show a case where the inequality constraint ( ) is active , i.e. ∃s * /|η (s * )| w = 2. We use the following physical parameters: length L = 32 cm, width w = 2 cm, thickness h = 0.1 mm, density ρ = 2000 kg/m 3 , Young's modulus Y = 2.016 GPa, Poisson's ratio ν = 0.4, uniform natural curvature κ N = -31.25 m -1 , and g = -g e x with the acceleration of gravity g = 9.81 m/s 2 . The clamp at s = 0 is oriented with As the Wunderlich elastic energy diverges when the constraint is active and since the solver only handles large inequalities, we numerically set the bound to 2 -with = 10 -8 The natural state is shown in purple, the initial state for the solving routine in green, and the equilibrium state in navy blue/turquoise. Please note that for visualization purposes, ribbons are displayed here with an increased thickness.

The nonlinear problem of finding equilibria of elastic ribbons may have several concurrent solutions. We illustrate this possibility by selecting parameter values for which at least five (plus two symmetric solutions) different solutions exist. We work with dimensionless quantities and set K N = -8, Γ = 85, ν = 0.4, and p = 20. Some of these solutions have proven difficult to find: shooting procedures based on the Sadowsky formulation failed to find them.

Consequently, we turn to the Wunderlich formulation along with the inequality constraint ( ) (that sometimes prevents the minimization procedure to go haywire) and therefore fix the dimensionless width, w/L = 0.03. We choose the orientation g = -e x and α = π/2. We use Ipopt to find the stable equilibrium solutions. Feeding the algorithm with different initial states, we isolate two such solutions, see (A) and (E) in figure . These two solutions are D and twisted (η(s) = 0) and each has a symmetrical twin, not shown here. Additionally, we isolate three planar (r(s) • e y = 0, ∀s) and untwisted (η(s) = 0, ∀s) solutions with the pure Newton routine, see (B), (C), and (D) in figure . These solutions are saddle points in the energy landscape and are therefore unstable: we had to use the pure Newton routine, solving for roots of the gradient of the energy instead of minimizing the energy itself.

. . Comparison with experiments

Currently we are developing a systematic protocol to validate simulators of slender structures. Generally speaking, by means of computer vision techniques, we are able to capture and reconstruct D deformed physical objects and project them on images with simulated results. For this work, we use ribbons sagging under their own weight, and the range in which such a problem is well defined, w/h 1 and 200 > Γ = 12(1 -ν 2 )ρgL 3 /(Y h 3 ) > 10 and L/w > 10, limits the dimensions of model specimens and puts constraints on materials manufacturing and characterization. We have found that the best option is to work with commercial polymeric plastics such as BoPP (Bi-Oriented Polypropylene), BoPET (Bi-Oriented Polyethylene Terephthalate), Acetate (butylene adipateco-terephthalate), or PVC (PolyVinyl Chloride) with Elastic moduli in the 2 -4 GPa range, and density around 10 3 kg/m 3 . These polymeric plastics are commercially available with thicknesses in the range [50 -500]µm. Such materials are elastic within the deformation we are expecting, but they are not completely isotropic. Furthermore their mechanical parameters vary from sample to sample or are difficult to measure (in particular Poisson's ratio). Hence it proves difficult, for a given experiment, to accurately simulate complex configurations, i.e. with large torsion.

Nevertheless, there are some basic configurations where we are able to quantitatively or/and qualitatively compare numerics with experiments. In this section we begin by showing a perfect agreement between our simulations, experiments, and theory for the simple cantilever case with no torsion and no natural curvature. Next, we compare a clamped naturally curved ribbon with the output of our simulator. In this case we found that, despite the match is not perfect, most of the geometrical characteristics are well generated by our synthetic output.

Cantilever experiment

A variation of the well-known Galileo cantilever is commonly used in mechanics of soft elastic objects to estimate material parameters. It consists of a horizontally clamped elastic object with a free end, deformed by the action of gravity, as depicted in the inset of figure . The global aspect ratio of the deformed shape (vertical divided by horizontal positions of free end) follows a master curve when plotted against the non-dimensional gravito-bending parameter Γ. To test our code, we set simulations of naturally flat ribbons with fixed width (w = 5cm), density (ρ = 1200kg/m 3 ), thickness (h = 1mm), and Poisson's Ratio (ν = 0.5), and span the Elastic modulus in the range Y ∈ [0.5 -100] MPa and the length within a range L ∈ [5 -90] cm. For each of the simulations we measure the horizontal (x) and vertical (y) components of the position of the ribbon free end, and plot the aspect ratio (y/x) as a function of Γ in fig. . Next, we numerically solve the equations for a planar elastica with bending rigidity Y h 3 w/(12(1 -ν 2 )) and show the result as a solid red line in figure . Finally, we run a series of experiments (green squares) with nominally flat ribbons. We first independently estimate Y /(1 -ν 2 ) for each ribbon as described in appendix B in [ ]. We measure the aspect ratios (y/x) of the deformed shapes and plot them as a function of Γ in green in fig. . A perfect agreement between experiments and ribbon and elastica simulations is observed for the large range of studied Γ values.

D deformations for naturally curved ribbons

We now turn to the more interesting case of naturally curved ribbons. The main challenges here are manufacturing such ribbons and accurately estimating both their Young modulus and Poisson's ratio. Inspired by the method described in [ ], we fabricate curly ribbons by winding flat specimens around a metallic cylinder, keeping them at degrees Celsius for hours, and finally letting them cool down while still attached to the cylinder. We tested several materials and observed that only PVC and BoPET adopt an acceptable circular shape, but we find that the resulting radius of curvature is in general larger than that of the metallic cylinder.

In figure we show two perpendicular views for one experimental realization with a sample made out of BoPET (right panels for each view). We measured the ribbon's thickness to be 100µm, its length 290mm, its natural curvature 26.7m -1 , its width 30mm, and density 1253kg/m 3 . For the ribbon's elastic parameters we estimated that the value of Y /(1-ν 2 ) is in the range of 2-4MPa. Since we are dealing with high torsion, Poisson's ratio plays an important role in the final shape, but we are not able yet to obtain a precise estimation of it, thus we use the classical value ν = 0.4. For capturing the images of the ribbon, we use a rig with a camera and a mirror. We calibrate the corresponding views as pinhole model cameras, as described in [ ]. This technique, together with a marker, allows us to measure experimentally the material frame orientation at the clamp (red, green and blue arrows in the experimental images represent {d 1 (0), d 2 (0), d 3 (0)} respectively). Thus the initial frame is fed into the simulator in addition to the elastic and geometrical parameters. As a result, we obtain a synthetic deformed shape for the ribbon (blue), see figure .

It is important to mention that for computing the simulated ribbon, we did fine tune the value of Y /(1 -ν 2 ) to be 4MPa, which is in the high limit of the estimated range. For a full validation of our model, in the future we will need to figure out a more accurate way to estimate the elasticity of the experimental ribbons.

. Conclusion and future work

We have proposed a robust and accurate computational approach dedicated to the statics of inextensible elastic ribbons, relying on the recent Wunderlich model. The key of our discrete model is to combine curvature-based elements with reduced parameterization, analytic computation of all kinematic and static terms based on power series, including the Hessian of energy, and robust secondorder minimization of energy under the Wunderlich constraint. We have tested our method on various ribbon configurations, from straight to highly curly, and have always observed excellent convergence in reasonable timings (a few seconds for 20 elements). Moreover, we have highlighted the value of carefully treating the Wunderlich energy together with the smoothness constraint. Less versatile approaches limited to the simpler Sadowsky energy, like current continuation methods, may lead to wrong ribbon configurations and/or serious numerical issues when encountering configurations close to the Wunderlich bounds. Finally, we have shown in our comparisons that a numerical model dedicated to ribbons is perfectly justified: rod models are inappropriate to capture the specific twist-curvature coupling of ribbons as they do not preserve inextensibility of the surface, while a generic shell model may reach the limit case of a thin ribbon but at the price of much heavier computations.

All of this makes our discrete model a practical tool suitable for exploring a wide range of new ribbon equilibrium configurations. In the future we hope to enrich the space of configurations by treating clamped-clamped boundary conditions, as well as self-contact and friction. We would also like to extend our experimental validation setup to a wider range of ribbon configurations. tionship written in the Cosserat frame,

v 1 = n 1 H 1 , v 2 = n 2 H 2 , v 3 = 1 + n 3 H 3 , (B. )
where we take H 1 = Y 2(1+ν) hw = H 2 and H 3 = Y hw. Note that the precise definition of the stiffness coefficient H 1,2 may differ in the literature where (near unity) prefactors are sometimes included [ ]. We recall that the Kirchhoff model is the limit case

1 H 1 → 0 , 1 H 2 → 0 , 1 H 3 → 0. (B. )
Next, the bend-twist constitutive relation is the same as in Kirchhoff's model,

m 1 = K 11 (κ 1 -κ N ) , m 2 = K 22 κ 2 , m 3 = K 33 κ 3 (B. )
with the stiffness coefficients given in ( ). Classical balance equations for linear and angular momentum read, in component form,

n 1 = n 2 κ 3 -n 3 κ 2 -f 1 m 1 = m 2 κ 3 -m 3 κ 2 + n 2 v 3 -n 3 v 2 (B. a) n 2 = n 3 κ 1 -n 1 κ 3 -f 2 m 2 = m 3 κ 1 -m 1 κ 3 + n 3 v 1 -n 1 v 3 (B. b) n 3 = n 1 κ 2 -n 2 κ 1 -f 3 m 3 = m 1 κ 2 -m 2 κ 1 + n 1 v 2 -n 2 v 1 . (B. c)
The system of equations is closed to the Darboux equations, listed in (A. ). The set of equations is then made dimensionless using L as unit length and K 11 /L 2 as unit force. Dimensionless equations are readily obtained by performing the following replacements,

K 11 = 1 K 22 = w h 2 K 33 = 2 c J 1 + ν (B. ) 1 H 1 = 1 + ν 6 h L 2 1 H 2 = 1 + ν 6 h L 2 1 H 3 = 1 12 h L 2 (B. )
where c J is defined in ( c). In typical situations where L w, the nondimensionalized coefficients H 1,2,3 1 and hence v 1 and v 2 are small and v 3 near unity, unless internal forces n 1,2,3 are large. As the typical scale of the internal force in our clamped-free situation is given by the weight, we conclude that as long as Γ L h 2 , the Antman and Kirchhoff models yield similar results.

In order to illustrate this fact, we compute the equilibrium solution for Γ ∈ (0; 100), L = 30 w, w = 10 h, and a dimensionless natural curvature K N = κ N L = 10 for the three models: Sadowsky, Kirchhoff, and Antman. We plot in figure B. the end-to-end distance in the plane orthogonal to the gravity direction, y 2 (1) + z 2 (1), as a function of Γ. We note that the curves from the Kirchhoff and Antman models are virtually indistinguishable, while they stand clearly apart from the curve of the Sadowsky model. In addition, we show in Figure B. the shapes given by the three models for Γ = 80 and observe again that shapes from the Kirchhoff and Antman models are indistinguishable.

Figure :

 : Figure : Stable equilibrium of a curled ribbon made of 5 elements, computed in 0.55 seconds with our method.

  Figure illustrates the geometric configuration of a ribbon.

  Figure : Base parameterization of a ribbon consisting of the centerline r, the Cosserat material frame R, and the scalar field η = tan θ.

Figure :

 : Figure :On the left, parameterization of the ribbon with s, v and η, on the flatten configuration. On the right, illustration of a singular point. Rulings passing at the points of arc lengths sa and s b are drawn in orange. The surface of the ribbon being developable and smooth, rulings are not allowed to intersect inside the ribbon.

  Figure : Left: a trapezoidal element (in white) shown on the flat configuration. Note that boundaries between the element and its neighbors coincide with rulings. Right: (top) Taking a constant curvature vector per element leads to discontinuous rulings, thus disconnected elements; (bottom) in contrast, a linear normal curvature and quadratic geodesic torsion per element ensures continuity of the rulings, and results in a smooth ribbon surface.

  deg n=0 λ n s n ,where deg is the degree of the polynomials. We instantiate this upper-bound with deg = 2 and find that R n s n ≤ 2 m if C(s) ≤ m ln 2 2 . In double precision (m = 52), this gives us s max =

  : E energy function, x starting point while convergence criterion not reached do A ← sysMat(E, x); d ← solve d {A • d = -∇E(x)} // Descent direction ; α ← stepSize(E, x, d)) // Step ; x ← x + α d; end Algorithm : Generic minimization algorithm

a

  x , d 1 (s = 0) = cos α e x +sin α e y , Iterations # {E} a # {∇E} b # {HE} c time (sec) Number of function evaluations b Number of gradient evaluations c Number of Hessian evaluations d Algorithm was stopped at 1000 iterations, even if the convergence was not reached.

Figure :

 : Figure : Three test cases of increasing complexity: (left) K N = 10 and Γ = 100, simulated with p = 10 elements, (middle) K N = 40 and Γ = 1200, simulated with p = 20 elements, and (right) K N = 100 and Γ = 7500, simulated with p = 20 elements. The natural shape is depicted in purple, and energy minimization algorithms are all initialized with the natural shape. For comparison purposes, the (simpler) Sadowsky energy is considered. Timing comparisons are provided in Table .

  Figure :A ribbon cannot be properly simulated using a rod model with a rectangular crosssection. Here we compare the equilibria generated by our discrete ribbon model (in blue) with that of the super-clothoid rod model with the same rectangular cross-section (in orange), when increasing their length. Both models coïncide when the configuration is planar (until L ≈ 29 cm), but substantially diverge when torsion enters the game (see accompanying video).

Figure :

 : Figure : Comparison of our method (in blue) with the FEniCS-Shell finite element code implementing the Naghdi shell model (six increasing width/thickness ratios depicted with different levels), under three different viewpoints (see accompanying video).

Figure :

 : Figure :Exactly same configuration found by our method (blue) and the shooting method (green).

Figure :

 : Figure : Configuration where bounds are reached. Left and middle: results obtained using the Sadowsky model, (left) without and (middle) with the Wunderlich constraint ( ). Right: result using the Wunderlich model and the constraint ( ), and zoom on the region where the constraint is active (convergence of the rulings close to the boundary). Without the Wunderlich constraint, simulating the Wunderlich model leads to fatal numerical issues.

Figure :

 : Figure : Five concurrent equilibria with K N = -8 and Γ = 85, using p = 20 elements.The natural state is shown in purple, the initial state for the solving routine in green, and the equilibrium state in navy blue/turquoise. Please note that for visualization purposes, ribbons are displayed here with an increased thickness.

Figure :

 : Figure :Aspect ratio y/x of the deformed shape of a sagging cantilever, as a function of the gravito-bending parameter Γ, evaluated in four ways: linear theory, y/x = Γ/8 valid for small Γ, planar elastica computations, elastic ribbons computations, and experiments.

Figure :

 : Figure : Side by side comparison between simulations (blue) and experiments (gray). Left pair shows a front view and right pair a side view.

  Figure B. : Comparison of three models for the statics of an elastic ribbon with L = 30 w, w = 10 h, K N = κ N L = 10, and ν = 0.4. We plot the end-to-end distance in the plane orthogonal to the gravity direction, y 2 (1) + z 2 (1), as a function of Γ = 12(1-ν 2 )ρgL 3 /(Y h 3 ).

Figure

  Figure B. : Comparison of three models for the statics of an elastic ribbon with L = 30 w, w = 10 h, K N = κ N L = 10, and ν = 0.4. We show the deformed shape of the structure at Γ = 12(1 -ν 2 )ρgL 3 /(Y h 3 ) = 80. We use the same color conventions as for figure B. , and, we draw the natural shape of the ribbon in purple. The shapes given by the Kirchhoff and Antman models are virtually the same, but clearly differ from the shape given by the Sadowsky model.

  [ ] R. Fosdick, E. Fried (Eds.), The Mechanics of Ribbons and Moebius Bands, ] M. Borri, C. Bottasso, An intrinsic beam model based on a helicoidal approximation-Part I: Formulation, International Journal for Numerical ] F. Bertails, B. Audoly, M.-P. Cani, B. Querleux, F. Leroy, J.-L. Lévêque, Super-Helices for Predicting the Dynamics of Natural Hair, ACM Transac-] Z. Pan, J. Huang, H. Bao, Modelling Developable Ribbons Using Ruling Bending Coordinates, Computing Research Repository. URL https://arxiv.org/abs/ . ] R. Casati, Quelques contributions à la modélisation numérique de structures élancées pour l'informatique graphique, Ph.D. thesis, Université ] J. Simo, D. Fox, On a stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization, Computer Methods in ] E. Doedel, H. B. Keller, J. P. Kernevez, Numerical Analysis and Control of Bifurcation Problems (I) Bifurcation in Finite Dimensions, International Journal of Bifurcation and Chaos ( ) ( ) -. [ ] E. Doedel, H. B. Keller, J. P. Kernevez, Numerical Analysis And Control of Bifurcation Problems (II): Bifurcation In Infinite Dimensions, International ] V. Duclaux, Pulmonary occlusions, eyelid entropion and aneurysm: A physical insight in physiology, Theses, Université de Provence -Aix-Marseille I (Nov. ). URL https://tel.archives-ouvertes.fr/tel-[ ] O. A. Arriagada, G. Massiera, M. Abkarian, Curling and rolling dynamics

	Springer,	, previously published in Journal of Elasticity Volume	,
	. Methods in Engineering Applied Mechanics and Engineering ( ) (	) ( ) (	-	)	. doi: . -.	/nme.
	[ ] B. Audoly, Y. Pomeau, Elasticity and Geometry: From hair curls to the nonlinear response of shells, Oxford University Press, . . URL http://dx.doi.org/ . [ ] R. Casati, F. Bertails-Descoubes, Super Space Clothoids, ACM Transac-[ ] D. N. Arnold, F. Brezzi, Locking-free finite element methods for shells, /nme. tions on Graphics ( ) ( ) . doi: . / . . Mathematics of Computation ( ) -.
	[ ] J. Chopin, A. Kudrolli, Helicoids, wrinkles, and loops in twisted ribbons, Phys. Rev. Lett. ( ) . [ ] F. A. Moebius, Über die bestimmung des inhaltes eines polyeders (On the [ ] M. A. Crisfield, G. Jelenić, Objectivity of strain measures in the geomet-[ ] A. Wächter, L. T. Biegler, On the implementation of an interior-point filter [ ] J. S. Hale, M. Brunetti, S. P. Bordas, C. Maurini, Simple and extensible rically exact three-dimensional beam theory and its finite-element imple-mentation, Proc. Royal Society of London, Series A ( ) ( line-search algorithm for large-scale nonlinear programming, Mathematical plate and shell finite element models through automatic code generation ) -. Programming ( ) ( ) -. doi: . /s ---y. tools, Computers & Structures ( ) -.
	determination of the volume of a polyhedron), Ber. Verh. Sächs. Ges. Wiss. ( ) -, see also Gesammelte Werke, Band II (Collected Works, vol. II), p. . Hirzel, Leipzig ( ). Computer Graphics Forum (Proc. Eurographics' ) ( ) ( ) . -. -ential Equations of the Möbius Band", Journal of Elasticity ( ) ( ) [ ] D. Pai, Strands: Interactive simulation of thin solids using Cosserat models, [ ] D. F. Hinz, E. Fried, Translation of Michael Sadowsky's paper "The Differ-[ ] Wolfram Research, Mathematica, Version , Champaign, IL, USA, .
	[ ] D. F. Hinz, E. Fried, Translation of Michael Sadowsky's paper 'An ele-
	mentary proof for the existence of a developable Möbius band and the
	attribution of the geometric problem to a variational problem', Journal of Elasticity ( ) ( ) -. [ ] W. Wunderlich, Über ein abwickelbares Möbiusband, Monatshefte für Mathematik ( ) ( ) -. doi: . /BF . tions on Graphics ( ) -. doi: . / . Grenoble Alpes ( ). . tact forces for deformable objects, Journal of WSCG ( ences, . ) -. [ ] J.-P. Demailly, Analyse numérique et équations différentielles, EDP Sci-[ ] J. Spillmann, M. Becker, M. Teschner, Non-iterative computation of con-URL http://www.theses.fr/ GREAM /document Journal of Bifurcation and Chaos ( ) ( ) -.
	Nature Materials ( [ ] A. Moore, T. Healey, Computation of elastic equilibria of complete möbius ) . doi: . /nmat . Mechanics ( ) ( [ ] L. Y. Adrianova, Introduction to Linear Systems of Diferential Equations, ) -. doi: . /j.ijnonlinmec. . . Vol. of Translations of Mathematical Monographs, American Mathe-. URL http://www.sciencedirect.com/science/article/pii/ matical Society, .
	bands and their stability, Mathematics and Mechanics of Solids S [ ] F. Bertails, Linear time super-helices, Computer Graphics Forum (Proc. ( ) of naturally curved ribbons, Soft Matter ( ) -.
	( [ ] E. J. Doedel, B. E. Oldeman, Auto-p: Continuation and bifurcation soft-) -. doi: . / . ware for ordinary differential equations, Tech. rep., Concordia University ( ). [ ] M. A. Dias, B. Audoly, "Wunderlich, Meet Kirchhoff": A General and Unified Description of Elastic Ribbons and Thin Rods, Journal of Elasticity ( ) ( ) -. doi: . /s --[ ] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, E. Grinspun, Discrete Eurographics' ) ( ). [ ] Z. Zhang, A Flexible New Technique for Camera Calibration., IEEE trans-elastic rods, ACM Transactions on Graphics (Proc. ACM SIGGRAPH' ) ( ) ( ) -. doi:http://doi.acm.org/ . / . [ ] M. Neher, An enclosure method for the solution of linear ODEs with poly-actions on pattern analysis and machine intelligence ( ) ( ) -. nomial coefficients, Numerical Functional Analysis and Optimization . URL http://www.cs.columbia.edu/cg/rods/ ( ) -. [ ] E. Reissner, On one-dimensional finite-strain beam theory: The plane [ ] H. Lang, J. Linn, M. Arnold, Multi-body dynamics simulation of geo-metrically exact cosserat rods, Multibody System Dynamics ( ) ( [ ] D. Goldberg, What every computer scientist should know about floating-problem, Zeitschrift für Angewandte Mathematik und Physik (ZAMP) ) -point arithmetic, ACM Comp. Surveys ( ) -. ( ) -. . doi: . /s ---x. -. URL https://doi.org/ . /s ---x [ ] J. Nocedal, S. Wright, Numerical Optimization, nd Edition, Springer, [ ] G. R. Cowper, The Shear Coefficient in Timoshenko's Beam Theory, Jour-
	[ ] B. Audoly, K. A. Seffen, Buckling of Naturally Curved Elastic Strips: The Ribbon Model Makes a Difference, Journal of Elasticity ( ) ( ) . doi: . /s --Methods in Applied Mechanics and Engineering ( ) -. algorithms . General considerations, IMA Journal of Applied Mathematics -y. element formulated on the special Euclidean group SE( ), Computer [ ] C. G. Broyden, The convergence of a class of double-rank minimization -[ ] V. Sonneville, A. Cardona, O. Brüls, Geometrically exact beam finite New York, NY, USA, . nal of Applied Mechanics ( ) ( ) -.
	[ ] D. E. Moulton, P. Grandgeorge, S. Neukirch, Stable elastic knots with no doi:http://dx.doi.org/ . /j.cma. . . . ( ) ( ) -. doi: . /imamat/ . . .
	self-contact, Journal of the Mechanics and Physics of Solids -. doi: . /j.jmps. . . . URL http://www.sciencedirect.com/science/article/pii/ ( ) [ ] S. S. Antman, Nonlinear problems of elasticity, nd Edition, Springer-S Verlag, New York, .
	[ ] T. G. Sano, H. Wada, Twist-induced snapping in a bent elastic rod and ribbon, Phys. Rev. Lett. ( ) . [ ] F. Renda, F. Boyer, J. Dias, L. Seneviratne, Discrete cosserat approach for [ ] A. I. Lurie, Theory of Elasticity, Foundations of Engineering Mechanics, multisection soft manipulator dynamics, IEEE Transactions on Robotics ( ) ( ) -. doi: . /TRO. . Springer, . . [ ] J. C. Simo, A finite strain beam formulation. The three-dimensional dy-namic problem. Part I, Computer Methods in Applied Mechanics and En-gineering ( ) ( ) -. [ ] S. S. Antman, Kirchhoff's problem for nonlinearly elastic rods, Quarterly [ ] Z. Shen, J. Huang, W. Chen, H. Bao, Geometrically exact simulation of inextensible ribbon, Computer Graphics Forum ( ) ( ) -. of Applied Mathematics ( ) ( ) -.

[ ] E. L. Starostin, G. H. M. van der Heijden, The shape of a Möbius strip, [ [ [ ] S. Goyal, N. Perkins, C. Lee, Non-linear dynamic intertwining of rods with self-contact, International Journal of Non-Linear [ [ [ ] P. M. Naghdi, Foundations of elastic shell theory, in: I. N. Sneddon, R. Hill (Eds.), Progress in Solid Mechanics, Vol. , North-Holland, , pp. -. [ [ [

Acknowledgments

The authors would like to thank Corrado Maurini for his help with FEniCS-Shell, Basile Audoly for fruitful discussions about the ribbon clamping conditions, Thibaut Métivet for interesting comments regarding our comparisons to thin elastic rods and shells, Mickaël Ly for rendering some of the figures and making the accompanying video, and José Bico and Benoit Roman for their support with materials and manufacturing. We used the free software Blender for rendering all our tridimensional figures. This work was supported in part by the ERC grant GEM (StG--).

Appendix A. Strong formulation of the equilibrium of a Sadowsky ribbon

We list the system of differential equations describing the equilibrium of a ribbon with natural curvature κ N in the approximation of the Sadowsky energy [ , ]. We normalize the physical quantities with the length L of the ribbon and its bending rigidity Dw. Besides the Poisson ratio ν, only two dimensionless parameters remain: Γ = ρhgL 3 /D and K N = κ N L. Vector components in the canonical basis are written with (x, y, z) indices and components in the material frame are written with (1, 2, 3) indices. The strong formulation of the equilibrium equations introduces the internal force n and moment m. We recall that we are in the case where

and that we replace κ 3 (s) by η(s) κ 1 (s). The nonlinear constitutive relation reads

As explained in [ ] using (A. ) would result in having to deal with a differentialalgebraic system, which we wish to avoid. Consequently we differentiate (A. )

and solve for κ 1 and κ 3 to write

These two equations are coupled with the classical (Kirchhoff-like) equations for position, equilibrium of force and moment

and material frame

The boundary condition at the clamp, s = 0, reads . ) and the boundary conditions are the free end, s = 1, are

Appendix B. Comparison with a shearable rod model

In this Section, we compare the Sadowsky model for elastic ribbon to Antman's rod model, which includes shear and extensional deformations [ ]. We consider a structure of length L with rectangular cross-section of width w and thickness h and we compute the equilibrium of the structure in the clamp-free setup used throughout this report. The Antman model is the D version of the Reissner model [ ] for elastic beams bent in the plane, and is an extension of the Kirchhoff model in the sense that both shear and extensional deformations are possible. The shear v 1 (s), v 2 (s) and extension v 3 (s) components are introduced to express the tangent to the centerline in the Cosserat frame as


S Perret

Y Curé

J Bleuse
email: joel.bleuse@cea.fr

L Grenet

R André

H Mariette

M A Green

Y Hishikawa

E D Dunlop

D H Levi

J Hohl- Ebinger

M Yoshita

Progress in Photo- A W Y Ho-Baillie

Optical Determination of the Band Gap and Band Tail of Epitaxial Ag 2 ZnSnSe 4 at Low Temperature

Keywords: numbers: 78.40.Fy, 78.20.-e, 78.55.-m, 78.47.jd

We report on the precise determination of both the band gap E g and the characteristic energy U of the band tail of localized defect states, for monocrystalline Ag 2 ZnSnSe 4 . Both photoluminescence excitation and timeresolved photoluminescence studies lead to E g = 1223 ± 3 meV, and U = 20 ± 3 meV, at 6 K. The interest of the methodology developed here is to account quantitatively for the time-resolved photoluminescence and photoluminescence excitation spectra by only considering standard textbook density of states, and state filling effects. Such an approach is different from the one most often used to evaluate the energy extent of the localized states, namely by measuring the energy shift between the photoluminescence emission and the excitation one -the so-called Stokes shift. The advantage of the present method is that no arbitrary choice of the low power excitation has to be done to select the photoluminescence emission spectrum and its peak energy.

I. INTRODUCTION

The occurrence of a band tail just below the band gap in the absorption spectra, when adding impurities to a pure crystal was first reported a long time ago by F. Urbach et al. 1 in AgCl crystals containing copper as an impurity. The observed band tail was well described by a single exponential and attributed to localized states created by the Cu-Ag cation disorder. Such analytic contribution to account for optical properties of doped semiconductors was then found to be verified and useful in a huge number of cases, especially when dealing with localized, substitutional impurities. More recently, as an example, the optical signature of impurity-impurity interactions in copper-containing II-VI alloy semiconductors has been identified by systematically comparing optical absorption and emission spectra 2 . We show here such an experimental, optical approach; we quantitatively analyze photoluminescence excitation data (reflecting the density of states) and time-resolved emission spectra (reflecting the evolution of the electron-hole population), using a single exponential shape for the localized states, as proposed by Urbach. This allows us to accurately determine the band gap and the energy extent of the band tail in the case of Ag 2 ZnSnSe 4 , a new promising material for solar energy conversion.

Indeed, by contrast to solar cells based on chalcopyrites such as Cu(In,Ga)(S,Se) 2 and zinc blende structures (CdTe), for which a power conversion efficiency (PCE) above 20% has been achieved 3 , kësterite-based solar cells such as Cu 2 ZnSn(S,Se) 4 (CZTSSe) only reached 12.6% PCE 4 . This limited performance is mainly related to a large open-circuit voltage (V OC ) deficit, as evidenced by the above CZTSSe record cell, which attains only 59% of the Shockley-Queisser limit for V OC , whereas it reaches 80% of this limit for both the short circuit current and the fill factor 4 .

Several hypotheses have been drawn to explain such a V OC limitation. Among them, the presence of a large band tail of localized states below the extended state bands is considered to be one of the main culprits, as it diminishes the effective bandgap, and consequently the available V OC 5 ; it also induces the localization of charge carriers, and thereby reduces their collection efficiency in devices.

The presence of large densities of intrinsic native point defects, as predicted by Density Functional Theory (DFT) calculations [6][7][8] , can directly influence the generation, separation and recombination of photo-generated electron-hole pairs. For CZTSSe, a large density of acceptor Cu Zn antisites and Cu vacancies V Cu has been calculated, given that both these point defects have the lowest formation energies 6 . Moreover, these point defects are stabilized by the formation of self-compensated neutral defect complexes made of triplets of exchange atoms such as [2Cu Zn + Sn Zn ]. In other words, the disorder in the cation sublattice, especially between Cu and Zn, strongly influences the occurrence of such defect clusters 7,8 , and is then directly responsible for the presence of potential fluctuations 9 : large amounts of these defects are likely to form band tailing, and believed to be one of the main reasons for the V OC deficit 5,10 . This issue is of particular importance, which is why we pursued our optical investigation in depth.

Such band tails have been studied qualitatively by optical spectroscopy: mainly photoluminescence (PL) as a function of temperature and excitation power [11][12][13][14][15][16][17] , but also by PL excitation spectroscopy (PLE), which reveals the density of states in the presence of potential fluctuations 18,19 , and by time-resolved photoluminescence spectroscopy (TR-PL), which gives information about the localization and transfer of carriers between the band tail states, as was done on CZTS single crystals 16,19 : all these data have been analyzed as a function of the order/disorder degree in the quaternary structure 19,[START_REF] Timmo | 13th E-MRS Spring Meeting[END_REF] .

The substitution of copper by silver is one of the proposed methods to reduce the density of defect states. The motivation for Ag 2 ZnSnSe 4 (AZTSe) or (Cu 2-x Ag x )ZnSnSe 4 derives from the reduced Ag Zn anti-sites concentration that was theoretically predicted by DFT calculations, as a consequence of a much larger formation energy [START_REF] Nakamura | 26th European Photovoltaic Solar Energy Conference (PV-SEC)[END_REF][START_REF] Yuan | [END_REF][23] for Ag Zn than for Cu Zn anti-sites. Then, if the anti-site defects between atoms of column I and II of the periodic table are indeed the main cause of band tailing, Ag 2 ZnSnSe 4 has the potential to significantly reduce the problem and to allow for larger V OC .

Such a substitution was successfully used [24][25][START_REF] Hages | Solar Energy Materials and Solar Cells[END_REF] to observe an improved efficiency, mainly due to a 8-10% improvement of V OC as compared to the baseline CZTSSe device. More precisely, it has been estimated, from room temperature optical data, that the band tailing effect is dramatically suppressed for AZTSe (pure-Ag) samples 25 .

In this work, we study single crystal, stoechiometric Ag 2 ZnSnSe 4 epilayers, as they offer a way to reduce the ambiguity on the nature of defects by eliminating grain boundaries. Complementary, low temperature optical experiments are compared : (i) the variation of the PLE intensity as a function of excitation photon energy; (ii) the evolution in time and energy of the TR-PL emission after a short pulse excitation; (iii) the variations of PL spectra with excitation power density. All these data are analyzed with a quantitative model that enables us to accurately determine both the band gap and the energy extent of the localized states band tail, which appears just below the band gap. These experimental data also evidence the transfer mechanisms that occur between the localized states, within this band tail.

II. EXPERIMENTAL MEASUREMENTS

The samples are grown in a MBE chamber using four Knudsen cells for Ag, Zn, Sn, and Se, without post treatment of the layer. The beam equivalent pressures of the four elements were measured by a gauge pressure meter positioned at the location of the sample. The AZTSe layers are grown on a (001) InAs substrate kept a temperature of 460 • C with a growth rate of 4 nm per minute. X-ray diffraction and Raman spectroscopy are used to check the absence of secondary phases, as was done on CZTSe [START_REF] Curé | [END_REF] .

For optical spectroscopy, both steady-state and timeresolved micro-photoluminescence experiments are carried out at cryogenic temperatures (6 K) in a helium-flow, optical cryostat. Electron-hole pair injection in the AZTSe layer is provided by a Ti:sapphire laser (Coherent Mira), operating either in steady-state mode or in pulsed mode with 200 fs-long pulses, and a repetition rate set by a cavity dumper at either 246 kHz (extracting one pulse every 222) or 501 kHz (every 109). The laser excitation is focused down to a 1.5-µm diam-eter spot on the sample with the same microscope objective (0.4 numerical aperture) that collects the luminescence signal. For PLE experiments, the excitation source is a 1 kW, halogen lamp coupled to a monochromator. The emission of the AZTSe layer is then spectrally dispersed by a 640 mm-focal length monochromator and detected by a silicon CCD (Andor DU420A-BRDD, for steady-state PL and PLE) or a silicon avalanche photodiode (Perkin-Elmer SPCM AQR-15, for TR-PL).

Low-temperature TR-PL spectra are measured at a series of about 20 different photon energies, evenly spaced from 1100 to 1370 meV, over a time window of about 600 ns, with 50-ps time bins. This gives a full view over the relaxation and recombination processes involved in the Ag 2 ZnSnSe 4 sample. For each of the nearly 12000 time bins, a lowresolution PL spectrum is therefore available for analysis. The temporal instrument response of the TR-PL setup, measured with the excitation laser pulses, presents a time resolution of about 400 ps, which is associated with the time jitter of the avalanche photodiode. To provide a clean measurement of the PL decay, the TR-PL measurements are conducted at a 3-µW average pumping power.

III. PHOTOLUMINESCENCE EXCITATION

PLE experiments allow us to measure the threshold of the optical absorption, for the direct evaluation of the AZTSe band gap at low temperature. Fig. 1 shows the PLE integrated emission intensity for a detection energy window between 1.05 and 1.09 eV. Knowing that such PLE data do not simply reflect the absorption spectrum but include also the relaxation processes to the detection energy, we measured these PLE spectra over various detection energy windows (not shown here): the spectra are similar whatever the detection energy with an uncertainty of ±3 meV, showing that the PLE spectra can be, in this instance, considered to be essentially proportional to the absorption coefficient [START_REF] Klingshirn | Experimental methods[END_REF] , itself proportional to the density of states.

The bandgap is therefore estimated from fitting the PLE data with a square root density of states in the band, above the band gap, typical of direct band gap semiconductors, and a density of band tail states in the gap -related to the localized states -described by a single exponential, as first proposed by Urbach 1 , which appears to be well adapted to the case at hand.

More precisely, we model the density of states of such a direct band gap semiconductor as:

d(E) =    A 2 E-E g U if E ≥ E g +U/2
A exp

E-E g U -1 2 if E ≤ E g +U/2
(1)

where E g is, as usual, the width of the band gap, and U is the characteristic energy that scales the extent of the Urbach tail density of localized defect states in the gap; A is a global amplitude factor. With the above definition, d(E) and its derivative are both continuous functions at the connection energy

IV. TIME-RESOLVED PHOTOLUMINESCENCE

Another way to estimate the band gap, and the extent of the band tail of localized states just below, is to perform a study of TR-PL spectra 5,19 . These spectra are reported in figure 2 for various time delays after a 1.55 eV, ≈ 200-fs long, pulsed excitation 29 . In figure 2(a), the rise of the TR-PL signal is plotted on a linear scale, which shows a maximum intensity at a 0.90-ns delay: this should correspond to the fast relaxation of the energy of hot carriers within the bands. Then, in figure 2(b), the subsequent signal decay exhibits a red energy shift that increases with time. This reveals the transfer of carriers within a band tail of states below the band gap, from the highest states to the lowest ones, on a much longer time scaleof the order of 100 ns.

These transfers between states of differing energies explain why, in such a case, a decay curve measured at a single photon energy can only accidentally give rise to a single exponential, and would account for only meaningless characteristic transi-tion times.

In the present case, the time-dependent spectra are plotted, on a semi-logarithmic scale, for time delays in a geometric series [START_REF]Figure 2(b) shows 9 spectra for delays from 1 to 200 ns, geometrically "equidistant[END_REF] . In such a scheme, if the PL decay were following a power law dependence on time [START_REF] Mourad | [END_REF]32 , the intensity maxima between two successive curves should be visually equidistant on this semi-logarithmic plot: it is not the case here, at least for times up to ≈ 25 ns, as the larger the time, the closer the spectra look, which evidences that the PL intensity peak decrease is slower than any power law, therefore slower than any exponential. For times larger than ≈ 25 ns, the data is not so clear, and a power law behaviour cannot be fully excluded.

From this, we can then infer that the carriers at the origin of the emission in the band tail are not fully independent, since the power law decay inherent to bimolecular recombination [START_REF] Mourad | [END_REF] is absent for delays below 25 ns. It is also worth mentioning here that, at room temperature, such a carrier redistribution in CZTSSe was found to occur over only a much shorter, 1 ns time scale, and on an energy range of only 10 meV 32 .

To get a more quantitative account of our data, we use the following analysis methodology: each of these lowresolution, TR-PL spectra is fitted with a function that is the product of the above density of states d(E) (eq. 1) with a Fermi occupation function:

f (E) = 1 + exp E -E F k B T -1 (2)
Since the band gap E g and the Urbach tail extent U do not depend on time, the fitting procedure is two-step: the nearly 2000 spectra, from 0.50 ns to 100 ns after the excitation pulse, are fitted [START_REF] Levenberg-Marquardt | algorithm is used. 34 The minimization search goes over the rectangle[END_REF] with the function f (E)d(E), keeping the values of E g and U constant, with only 3 parameters varying with time: A, E F , and k B T . The cumulated errors between data and fitted values over these 2000 spectra then represent a global error for the chosen (E g , U) pair. This is then repeated to find the optimal E g and U values that minimize this global error 34 .

The results of this two-step procedure are given in figure 2, showing the fits (full lines) together with experimental data (points): a bandgap E g = 1226.3 ± 3 meV, and an Urbach tail characteristic energy U = 16.8 ± 1 meV are precisely found, along with the time variation of the parameters A, E F , and k B T . This determination of E g and U is in good agreement with the one deduced above from the PLE experiments. As in figure 1, the grey vertical line in figure 2 covers the energy interval from 1220.4 to 1226.3 meV.

One must however take note that we do not currently understand the exact nature of the emission that is always observed on the high energy side of the spectra (E > 1.3 eV) ; this emission is also observed in steady-state PL spectra (see figure 3). As it is clearly distinguished from the emission peak, it is not taken into account in our data analysis, .

V. DISCUSSION

It is clear from this analysis that the maxima of the TR-PL peak, for very short time delays after the pulse, are at energies larger than E g (figure 2(a)), and correspond to the recombination of hot carriers, before and along energy relaxation from phonon emission. On the other hand, for delays larger than ≈ 3 ns, the maxima of the TR-PL peak occur at energies lower than E g , and are related to carrier recombination involving localized states in the band tail.

Considering the parameters extracted from the model, the Fermi energy E F is -for delays larger than 1 ns -a decreasing function of time that closely follows the peak energy maxima, as expected from the characteristics of our model function f (E)d(E). This superposition of E F with the energy variation of the peak maxima is confirmed with steady-state PL spectra as a function of power excitation (see fig. 3). There, by varying the power density over five orders of magnitude, the PL maxima are shifted by ≈ 100 meV, from 1145 to 1245 meV, as are the values of E F (see figure 4).

In both sets of experimental data -steady-state PL as a function of power excitation and TR-PL, the emission spectra can be related to the steady-state or instantaneous carrier density. PL peak and Fermi energies, as well as carrier temperatures should therefore correspond when the carrier densities match. This correspondence is attempted in figure 4 where we adjust the excitation power and time scales to show that peak energy and full width at half maximum (FWHM) pairs exhibit similar trends between steady-state and time-resolved PL. This shows that, for delays larger than ≈ 5 ns, the carrier populations in TR-PL measurements are probably close to thermal equilibrium.

Our second fit parameter, the "temperature" k B T , which is directly related to the peak FWHM, decreases as a function of time and saturates, at delays larger than 5 ns, to a value of ≈ 6 meV, which is equivalent to T ≈ 70 K: this temperature, much larger than the bath one (6 K), shows the limitation of our model to account for weakly-populated discrete defect states in the band tail. Nevertheless, these values of k B T be-tween 1 and 200 ns can be, used to estimate the FWHM (full green line in figure 4), which shows values between 35 and 70 meV, fully accounting for the TR-PL experimental data. These are in very good agreement with the peak FWHM deduced from steady-state PL spectra (green squares), and a factor of 2 smaller than the FWHM usually measured on pure copper kësterites compounds 19,[START_REF] Timmo | 13th E-MRS Spring Meeting[END_REF] , which lie between 70 and 110 meV, even for the most ordered CZTS monocrystalline samples.

VI. CONCLUSION

This work reports on a precise determination of the band gap of monocrystalline Ag 2 ZnSnSe 4 at low temperature. Both photoluminescence excitation and time-resolved photoluminescence lead to a value of E g at 6 K equal to 1223 ± 3 meV. It is worth noting that such values are much lower than the ones reported at room temperature for polycristalline thin films: figuring diffuse reflectance data in a Tauc plot, W. Gong et al [START_REF] Gong | 19th International Conference on Ternary and Multinary Compounds (ICTMC)[END_REF] obtained a value of E g = 1.34 eV, whereas Gershon et al. [START_REF] Gershon | [END_REF] deduced a value of 1.35 eV from External Quantum Efficiency data, in agreement with their room-temperature PL emission maximum. Variations of the gap values can however be strongly dependent on the selenization temperature, i.e. Se content, as pointed out by Jiang et al. 37 who reported band gap values between 1.33 eV and 1.59 eV. Our results nevertheless indicate that the usual behavior of semiconductor band gapa decrease with temperature because of electron-phonon coupling, does not hold here: further studies as a function of temperature are under way to elucidate this unexpected behavior.

This work also presents a quantitative methodology to precisely deduce the band gap E g and the characteristic energy U for the extent of the tail of localized states, from the PLE and TR-PL data. Such an approach is different from the one usually used, namely measuring the energy shift between the PL emission and the absorption PLE threshold (the Stokes shift). The advantage of the present method is that no arbitrary choice of the low power excitation has to be done to select the PL emission spectrum and its peak energy (see figure 4).

Finally, with this systematic optical spectroscopy study and the above analysis methodology, the characteristic energy of the band tail in this new Ag 2 ZnSnSe 4 material is found to be 20 ± 3 meV, which leads to a photoluminescence FWHM that is half those usually reported for even the most-ordered Cu-based kësterites compounds.

FIG. 1 .

 1 FIG. 1. PLE experimental data (dots) at 6 K. The integrated emission intensity is detected in the 1.05-1.09 eV energy range as a function of the excitation energy. Least squares fitting of the model density of states d(E) (solid line -see text) yields values for the band gap E g = 1220.4 meV and band tail characteristic energy U = 23 meV. The fit deduced from the TR-PL data is also shown (dashed line).

FIG. 2 .

 2 FIG. 2. Low temperature TR-PL spectra. Experimental data (dots) and fitting curves (solid lines) for various time delays after the pulse excitation. (a) Signal rise for short delays, on a linear scale. The TR-PL signal peak is maximum at a 0.90 ns delay. (b) Signal decay, on a semi-logarithmic scale, for time delays that follow a geometric series from 1 to 200 ns. TRPL spectra exhibit red energy shifts that increase with time, revealing the transfer of carriers within the band, up to ≈ 2 ns, and between the localized states in the band tail, for longer delays.

FIG. 3 .

 3 FIG. 3. Steady-state PL spectra as a function of excitation power -over 5 orders of magnitude, on a semi-logarithmic scale. The dashed lines corresponds to fits with the model function f (E)d(E), taking the E g and U values obtained from the TR-PL data; the E F parameters extracted from these fits are reported in figure 4. The exponential decrease on the low energy side of the spectra correspond to the characteristic energy U of the Urbach tail density of localized defect states in the gap.

FIG. 4 .

 4 FIG.4. Steady-state PL data as a function of excitation power (lower horizontal, logarithmic axis): the PL peak energy positions (blue dots, left vertical axis) are closely followed by the values of E F (open red circles), which are deduced from the spectral fits in figure3; the experimental values of the FWHM are also reported (green squares, right vertical axis). Superimposed are TR-PL data as a function of time delay (upper horizontal, logarithmic axis): the TR-PL peak energy positions (blue curve, left vertical axis), and peak FWHM (green curve, right vertical axis) are taken from the fitted curves. The two horizontal, logarithmic axes give the correspondance between delay and excitation power, as both correspond to variations of the electronhole density. The "jump" in TR-PL curves around 36 ns delay comes from a weak, parasitic re-excitation.

ACKNOWLEDGMENTS

This research was supported by the french ANR program "Carnot Energies du Futur", under the project name CAZTS.