Simon Apers
email: smgapers@gmail.com

Laurent Miclo

Quantum walks, the discrete wave equation and Chebyshev polynomials

A quantum walk is the quantum analogue of a random walk. While it is relatively well understood how quantum walks can speed up random walk hitting times, it is a long-standing open question to what extent quantum walks can speed up the spreading or mixing rate of random walks on graphs. In this expository paper, inspired by a blog post by Terence Tao, we describe a particular perspective on this question that derives quantum walks from the discrete wave equation on graphs. This yields a description of the quantum walk dynamics as simply applying a Chebyshev polynomial to the random walk transition matrix. This perspective decouples the problem from its quantum origin, and highlights connections to earlier (non-quantum) work and the use of Chebyshev polynomials in random walk theory as in the Varopoulos-Carne bound. We illustrate the approach by proving a weak limit of the quantum walk dynamics on the lattice. This gives a different proof of the quadratically improved spreading behavior of quantum walks on lattices.

Introduction

Random walks (and Markov chains, more generally) show up in many areas of science. A particularly relevant question is how fast random walks spread on graphs. This is related to topics such as the thermodynamic equilibration of physical systems, the time complexity of Markov chain Monte Carlo algorithms, or the diffusion speed of diseases and rumours in networks. It is typically quantified by the random walk mixing time [START_REF] Levin | Markov chains and mixing times[END_REF], which captures the time that it takes for a random walk to converge to its stationary distribution.

Quantum walks are the quantum mechanical analogues of random walks. They are a key tool in the design of quantum algorithms [START_REF] Ambainis | Quantum walks and their algorithmic applications[END_REF], and they have been a source of interesting mathematics. Building on a long line of work, it was recently established [START_REF] Ambainis | Quadratic speedup for finding marked vertices by quantum walks[END_REF][START_REF] Piddock | Quantum walk search algorithms and effective resistance[END_REF] that quantum walks can quadratically speed up the random walk hitting time, i.e., the expected time to hit a subset of vertices when starting from the stationary distribution. This can be thought of as a graph-theoretic version of the quadratic speedup in Grover's quantum search algorithm [START_REF] Lov | A fast quantum mechanical algorithm for database search[END_REF]. The random walk mixing time is somehow dual to the hitting time, requiring instead to converge to the stationary distribution starting from a single state. A significantly smaller body of work has been devoted to proving a quadratic quantum walk speedup in the mixing time, paralleling the quadratic speedup in the hitting time. Such a speedup was proven in a number of relevant but very restricted settings, including settings where the underlying graph is highly symmetric [START_REF] Aharonov | Quantum walks on graphs[END_REF][START_REF] Peter | Almost uniform sampling via quantum walks[END_REF][START_REF] Peter | Quantum speedup of classical mixing processes[END_REF], very poorly connected [A19], or when we have access to a simulated annealing schedule [START_REF] Wocjan | Speedup via quantum sampling[END_REF][START_REF] Harrow | Adaptive quantum simulated annealing for Bayesian inference and estimating partition functions[END_REF].

Despite advances in these restricted settings, it seems that progress on a general quantum walk speedup in the random walk spreading rate or mixing time has largely stalled. We revisit this question using an alternative perspective on the quantum walk dynamics, based on Chebyshev polynomials. 1 Rather than following the standard exposition that is usually found in the quantum algorithms literature (see e.g. [START_REF] Childs | Lecture notes on quantum algorithms[END_REF]), we present an alternative perspective based on the discrete wave equation on graphs. This approach is inspired and largely follows a blog post on the discrete wave equation by Terence Tao [START_REF] Tao | Discretised wave equations[END_REF] (which, on its turn, was inspired by discussions with Yuval Peres). It also connects to earlier work on the discrete wave equation on graphs by Friedman and Tillich [START_REF] Friedman | Wave equations for graphs and the edge-based Laplacian[END_REF]. Surprisingly, the connection with quantum walks seems to not have been made earlier, and formalizing this is one of the motivations for writing this paper. As the main technical contribution and proof-of-concept, we use the analysis of Chebyshev polynomials to prove a weak limit on the quantum walk spreading behavior. In particular, this implies that quantum walks have a ballistic spreading rate on lattices, in correspondence with earlier findings [START_REF] Mackay | Quantum walks in higher dimensions[END_REF][START_REF] Baryshnikov | Two-dimensional quantum random walk[END_REF] that used very different techniques.

Throughout the manuscript, we put care in giving a self-contained exposition, with the hope of rejuvenating interest in the problem and attracting a broader (and potentially non-quantum) audience.

Random walks from the diffusion equation

As a warm-up, consider the random walk dynamics

∀ t ∈ Z + , u(t + 1) = P u(t) (1)
where P is the random walk transition matrix on a graph. We will argue how these dynamics correspond to a discrete diffusion equation on graphs.

To this end, consider the diffusion equation

u = 1 2d ∆u (2)
over R d with ∆ the Laplacian operator. A discrete space and time analogue of (2) is given by (1) with P the random walk transition matrix on the lattice Z d . Indeed, one gets an approximation of (2) by introducing a small parameter ε > 0 and replacing the time set Z by εZ and the lattice Z d by √ εZ d . Consider then the evolution

∀ t ∈ εZ + , u ε (t + ε) = P √ ε [u ε (t)] (3)
where P √ ε is the random walk transition matrix on √ εZ d . This evolution coincides with (1) for ε = 1. Now assume e.g. that the initial condition u(0) of ( 2) is continuous and has a compact support and that the initial condition of (3) is just the restriction of u(0) to √ εZ d . Then the following well-known approximation result holds, as a consequence of Donsker's invariance principle [START_REF] Donsker | Justification and extension of Doob's heuristic approach to the Kolmogorov-Smirnov theorems[END_REF][START_REF] Mörters | Brownian motion[END_REF]: for any given (t,

x) ∈ R + × R d , we have lim ε→0 + u ε (t ε , x ε ) = u(t, x)
where for any ε > 0, t ε ∈ εZ + and x ε ∈ √ εZ d are such that lim ε→0 + t ε = t and lim ε→0 + x ε = x. This motivates identifying a random walk on the lattice Z d with a discrete diffusion equation. Now, through a small leap of faith, we can generalize this by letting P be the random walk transition matrix on a general graph. We identify the corresponding dynamics in (1) with the discrete wave equation on that graph.

Quantum walks from the wave equation

Now we consider the wave equation instead of the diffusion equation, given by

ü = 1 4d ∆u over R d . Equivalently, u = v, v = 1 4d ∆u. ( 4
)
We will consider a particular2 discrete space and time approximation of the diffusion equation, given by

∀ t ∈ Z + , u(n + 1) = P u(n) + v(n) v(n + 1) = P v(n) -(I -P 2 )u(n) (5) for u(n), v(n) ∈ ℓ 2 (Z d ).
The following exposition largely follows [START_REF] Tao | Discretised wave equations[END_REF].

Discrete approximation. We can again motivate this approximation by considering the case with P the random walk transition matrix on Z d . We introduce a small parameter ε > 0 and replace the time set Z by εZ and the lattice Z d by εZ d . Consider then the evolution

∀ t ∈ εZ + , u ε (n + 1) = P ε u ε (n) + εv ε (n) v ε (n + 1) = P ε v ε (n) - 1 ε (I -P 2 ε )u ε (n) (6)
where P ε is the random walk transition matrix on εZ d . This evolution coincides with (5) for ε = 1. Now assume e.g. that the initial condition (u(0), v(0)) of ( 4) is continuous and has a compact support and that the initial condition of ( 5) is just the restriction of (u(0), v(0)) to εZ d . Then the following again holds: for any given (t,

x) ∈ R + × R d , we have lim ε→0 + u ε (t ε , x ε ) = u(t, x), lim ε→0 + v ε (t ε , x ε ) = v(t, x)
where for any ε > 0, t ε ∈ εZ + and x ε ∈ εZ d are such that lim ε→0 + t ε = t and lim ε→0 + x ε = x. This argument motivates identifying the dynamics (5) on Z d with the discrete wave equation on Z d . More generally, we identify (5) with P the random walk transition matrix on a general graph with the discrete wave equation on that graph.

Unitary dynamics. It is easy to check that the following "energy" quantity

∥ I -P 2 u(n)∥ 2 + ∥v(n)∥ 2 is conserved. Indeed, if we define w(n) = √ I -P 2 u(n)
, then we can equivalently write (5) as

w(n + 1) = P w(n) + I -P 2 v(n) v(n + 1) = P v(n) -I -P 2 w(n).
The corresponding propagator

U = P √ I -P 2 - √ I -P 2 P (7)
is then unitary, and hence conserves the norm or energy ∥w(n)∥ 2 + ∥v(n)∥ 2 . Such a unitary that encodes an operator P as one of its blocks is called a "unitary dilation" in the mathematics literature.

The n-step propagator U n can be neatly expressed as

U n = T n (P ) √ I -P 2 U n-1 (P ) - √ I -P 2 U n-1 (P ) T n (P ) . ( 8
)
where (T n ) n∈Z + and (U n ) n∈Z + are respectively the families of Chebyshev polynomials of the first and second kind.

Quantum walks and block encodings

Arguably, such unitary dynamics are a natural map between quantum states, in which case we propose to call the resulting dynamics a quantum walk on the underlying graph. In fact, the unitary dilation in (7) is bread and butter in the quantum computing literature -there it is usually called a block encoding [START_REF] Chakraborty | The power of block-encoded matrix powers: improved regression techniques via faster Hamiltonian simulation[END_REF] of the operator P . Such block encodings are central to our current understanding of quantum algorithms -see e.g. the recent "grand unification of quantum algorithms" in the form of the quantum singular value transformation [START_REF] Gilyén | Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics[END_REF]. One of the first (if not the first) demonstrations of such block encoding was in the seminal work on quantum walks by Szegedy [START_REF] Szegedy | Quantum speed-up of Markov chain based algorithms[END_REF]. Szegedy constructed a unitary quantum walk operator W , different from the propagator U in (7), but resembling it in that 1. W is a unitary dilation of the random walk operator P , and 2. W n is a unitary dilation of T n (P ).

Szegedy used some additional properties of their particular construction of W . However, we now understand that any operator satisfying the properties 1. and 2. can be used as a quantum walk operator, allowing one to reproduce the algorithmic results on quantum walks referenced earlier [START_REF] Apers | A unified framework of quantum walk search[END_REF]. This motivates our interpretation of the discrete wave dynamics as a quantum walk -in fact, it is surprising that this connection seems to not have been made earlier.

We note that, while the u-and v-components (or w-and v components) could be interpreted as position and momentum in the original wave equation, they lose that meaning in the quantum walk interpretation. Rather, the two components should be interpreted as the components of a state with two chiralities.

Varopoulos-Carne bound and quantum fast-forwarding

A nice demonstration of this connection is the use of a particular identity for Chebyshev polynomials in random walk theory on the one hand (also discussed in [START_REF] Tao | Discretised wave equations[END_REF]), and in quantum algorithms design on the other hand. The identity is the following: for any self-adjoint contraction P it holds that

P n = k∈Z q n (k)T |k| (P ), (9)
where q n (k) is the probability that a random walk on the line, starting from the origin, is at position k after n steps.

In random walk theory, this identity is used to prove the Varopoulos-Carne bound (see [START_REF] Varopoulos | Long range estimates for Markov chains[END_REF][START_REF] Keith | A transmutation formula for Markov chains[END_REF] and a nice discussion in [START_REF] Lyons | Probability on trees and networks[END_REF]). The bound states that for any pair of vertices x, y it holds that

P n (x, y) ≤ 2 exp(-d(x, y) 2 /(2n)),
where d(x, y) is the shortest path distance between vertices x and y. The main idea of the bound is to apply Hoeffding's bound to the coefficients q n (k), showing that

P n ≈ |k|≤O( √ n) q n (k)T |k| (P ). ( 10
) Now notice that [T |k| (P )](x, y) = 0 if d(x, y) > k, and so P n (x, y) ≈ 0 when d(x, y) > O( √ n).
In quantum walk algorithms, the identities (9) and (10) were used recently in a technique called quantum fast-forwarding [START_REF] Apers | Quantum fast-forwarding Markov chains and property testing[END_REF]. Combining the expression (8) for the quantum walk propagator with the approximation (10), we see that

|k|≤O( √ n) q n (k)U |k| ≈ P n ⋆ ⋆ P n , (11)
where the ⋆ represent non-relevant entries. Through an additional quantum algorithmic technique called "linear combination of unitaries" [START_REF] Childs | Hamiltonian simulation using linear combinations of unitary operations[END_REF], it is possible to algorithmically implement a unitary dilation of (11) (and hence of P n ) while effectively only implementing O( √ n) quantum walk steps. This explains the terminology "quantum fast-forwarding".

Quantum walks on lattices

The preceding discussion shows that we can understand quantum walks through the lense of Chebyshev polynomials. We demonstrate this connection by studying such Chebyshev polynomials for the particular case of lattices. Specifically, we are interested in the dynamics

U n e 0 0 = T n (P ) e 0 - √ I -P 2 U n-1 (P ) e 0 ,
where e 0 is the indicator on some vertex of the lattice. This corresponds to the discrete wave equation with initial conditions u(0) = e 0 and v(0) = 0. Equivalently, it corresponds to an n-step quantum walk starting from e 0 . We will be interested in the spreading behavior of the quantum walk, and will focus on the u-component of the resulting state, i.e., the vector T n (P )e 0 . This vector induces a (subnormalized) measure ν on the lattice by setting ν x = [T n (P )] 0,x 2 , and this measure will be the object of our study. Operationally, it corresponds to the probability of measuring the quantum walk in a certain chirality at vertex x (see discussion at the end of Section 3.1).

Quantum walk on the line

On the infinite line, the T n (P ) component of the quantum walk takes a particularly simple form. Specifically, consider the line with set of vertices Z and commuting operators

Q 1 and Q -1 1 defined by Q 1 e x = e x+1 and Q -1 1 e x = e x-1
, where for any x ∈ Z, e x is the function on Z taking the value 1 at x and 0 elsewhere. The random walk transition matrix is

P = (Q 1 + Q -1
1 )/2, and the random walk distribution [P n ] 0,• after n steps essentially corresponds to a binomial distribution around 0 with standard deviation O( √ n). In contrast, if we use the formula T n ((x + x -1 )/2) = (x n + x -n )/2 we can rewrite the Chebyshev polynomial T n (P ) as

T n (P ) = T n Q 1 + Q -1 1 2 = Q n 1 + Q -n 1 2 . ( 12
)
This implies that [T n (P )] 0,x = e n (x)/2 + e -n (x)/2, and so the quantum walk shows two peaks precisely at a distance n from the origin. These dynamics represent the ballistic nature of the quantum walk dynamics, as opposed to the diffusive nature of the random walk dynamics. The hope is that this improved spreading rate persists in more general graphs as well.

If instead of the simple random walk we consider the lazy random walk on the cycle, P L = (I + P )/2, there is no longer a simple expansion of T n (P L ). Instead, we retrieve the fickle (but still ballistic) behavior usually associated to the quantum walk on the line. The behaviors of P n L and T n (P L ) are illustrated in Fig. 1.

Quantum walk on the lattice

As a more involved example, consider a random walk on the 2-dimensional lattice with vertices

V = {(x, y) ∈ Z 2 } and commuting generators Q 1 , Q -1 1 , Q 2 , Q -1 2 defined by Q 1 e
x,y = e x+1,y and Q 2 e x,y = e x,y+1 , and similarly for the inverses. If we set

P 1 = (Q 1 + Q -1 1 )/2 and P 2 = (Q 2 + Q -1
2 )/2, then the random walk transition matrix is given by

P = P 1 + P 2 2 .
In Fig. 2 we compare the measure induced by the random walk transition matrix P n and that of the Chebyshev polynomial T n (P ). To get a handle on the operator T n (P ), we will use the following expansion of a degree-n bivariate Chebyshev polynomial into univariate degree-n Chebyshev polynomials:

T n X + Y 2 = (p,q)∈ 0,n a n,p,q T p (X)T q (Y ), (13)
for a set of real coefficients (a n,p,q ) p,q∈ 0,n . Since T n (P ) = T n ((P 1 + P 2 )/2) is a degree-n polynomial in the commuting operators P 1 and P 2 , we can use the same expansion to rewrite

T n P 1 + P 2 2 =
(p,q)∈ 0,n a n,p,q T p (P 1 )T q (P 2 )

= (p,q)∈ 0,n a n,p,q Q p 1 + Q -p 1 2 Q q 2 + Q -q 2 2 ,
where we reused the fact (Eq. ( 12)) that

T p (P 1 ) = (Q p 1 + Q -p 1 )/2. Now notice that [T n (P )] (0,0),(p,q) = a n,|p|,|q| /2 k(|p|,|q|) ,
where k(|p|, |q|) ∈ {0, 1, 2} denotes the number of nonzero entries among |p|, |q|. We can hence lower bound the probability for a quantum walk, starting from the origin (0, 0), to reach a vertex (p, q) by

[T n (P )] (0,0),(p,q) 2 = a 2 n,|p|,|q| /2 2k(|p|,|q|) ≥ a 2 n,|p|,|q| /16 (14)
(since k(p, q) ≤ 2). These coefficients are fully determined by the Chebyshev expansion in (13), and so the problem reduces to bounding these coefficients. More specifically, we will consider the induced measure µ n on [-1, 1] 2 given by

µ n p,q a 2 n,|p|,|q| δ (p/n,q/n) p,q a 2 n,|p|,|q|
where δ (p/n,q/n) is the Dirac mass at (p/n, q/n) ∈ [0, 1] 2 . As our main technical contribution, we prove a weak limit theorem on µ n . Theorem 1. There exists a continuous measure µ on [-1, 1] 2 such that µ n → µ.

As the most importance consequence of this weak limit, it proves that the quantum walk spreads ballistically on the lattice -after n steps, the expected distance from the origin is Ω(n). Using very different machinery, similar conclusions were found in earlier work [START_REF] Mackay | Quantum walks in higher dimensions[END_REF][START_REF] Baryshnikov | Two-dimensional quantum random walk[END_REF].

Proof of weak convergence on the lattice

Our goal is to investigate the weak convergence for n large of µ n toward a law µ on [0, 1] 2 . To do so, it will be more convenient to consider the measure on [0, 1] 2 given by γ n p,q∈ 0,n a 2 n,p,q δ (p/n,q/n)

We will show there exists a probability measure γ ̸ = 0 on [0, 1] 2 such that for large n, γ n weakly converges toward γ. Since for any n ∈ Z + ,

µ n = γ n γ n ([0, 1] 2 ) we deduce the weak convergence of µ n toward µ γ/γ([0, 1] 2 ) = γ.
The method will consist in showing the convergence of the (bivariate) even moments of γ n toward those of γ. More precisely, for any K, L ∈ Z + , consider the polynomial function

φ K,L : [0, 1] 2 ∋ (x, y) → x 2K y 2L
We are to find a measure γ ̸ = 0 on [0, 1] 2 such that for any

K, L ∈ Z + , lim n→∞ γ n [φ K,L ] = γ[φ K,L ]
Taking into account the Stone-Weierstrass theorem (applied on the compact set [0, 1] 2 with respect to the algebra generated by the φ K,L , for K, L ∈ Z + ), this convergence is sufficient to deduce the weak convergence of γ n toward γ.

The proof is split in three main parts corresponding to the following subsections.

Second moment as an integral

We start by approximating the second moments of µ n by some simpler integral. Consider the function h(x, y) := T n cos(x) + cos(y) 2 = p,q∈ 0,n a n,p,q cos(px) cos(qy),

where the last equality follows from Eq. ( 13) and the definition of the Chebyshev polynomials. We show the following approximation.

Proposition 1. We have

lim n→∞ γ n [φ K,L ] - 1 π 2 n 2(K+L) [0,2π] 2 ∂ K x ∂ L y h(x, y) 2 dxdy = 0.
The proof consist of the two following results.

Lemma 1. We have

1 π 2 n 2(K+L) ∂ K x ∂ L y h(x, y) 2 dxdy = 1 n 2(K+L) p,q∈ 0,n a 2 n,p,q p 2K q 2L w K (p)w L (q) = γ n [φ K,L w K (n•) ⊗ w L (n•)]
where

∀ k ∈ Z + , ∀ r ∈ Z + , w k (r)     
2 , if k is even and r = 0 0 , if k is odd and r = 0 1 , otherwise Proof. We rewrite

∂ K x ∂ L y h(x, y) = p,q∈ 0,n a n,p,q p K q L ξ K (px)ξ L (qy) (16)
where ξ r (z) denotes the r-th derivative cos (r) (z). Now recall the orthogonality relations:

∀ p, p ′ ∈ Z + , 1 2π 2π 0 ξ r (px)ξ r (p ′ x) dx =      1 , if p = p ′ = 0 and r ≡ 0[2] 1/2 , if p = p ′ ≥ 1 0 , otherwise
Taking the square of ( 16) and integrating, we get

1 (2π) 2 ∂ K x ∂ L y h(x, y) 2 dxdy = p,q,p ′ ,q ′ ∈ 0,n a n,p,q p K q L a n,p ′ ,q ′ (p ′ ) K (q ′ ) L 1 (2π) 2 ξ K (px)ξ L (qy)ξ K (p ′ x)ξ L (q ′ y) dxdy = p,q,p ′ ,q ′ ∈ 0,n a n,p,q a n,p ′ ,q ′ p K (p ′ ) K q L (q ′ ) L 1 2π ξ K (px)ξ K (p ′ x) dx 1 2π ξ L (qy)ξ L (q ′ y) dy = p,q∈ 0,n a 2 n,p,q p 2K q 2L 1 2π ξ 2 K (px) dx 1 2π ξ 2 L (qy) dy = 1 4 p,q∈ 0,n a 2 n,p,q p 2K q 2L w K (p)w L (q) = 1 4 n 2(K+L) γ n [φ K,L w K (n•) ⊗ w L (n•)].

■

We deduce:

Corollary 1. We have for any n ∈ Z + , p,q∈ 0,n a 2 n,p,q ≤ 4 Proof. Considering the case K = L = 0 in Lemma 1, we get

1 π 2 [0,2π] 2 T 2 n cos(x) + cos(y) 2 dxdy = p,q∈ 0,n w 0 (p)w 0 (q)a 2 n,p,q ≥ p,q∈ 0,n a 2 n,p,q
To conclude to the announced bound, note that the l.h.s. is bounded above by 4, since

T n (z) takes its values in [-1, 1] for any z ∈ [-1, 1]. ■
The next step enables to end the proof of Proposition 1.

Lemma 2. We have

lim n→∞ |γ n [φ K,L w k (n•) ⊗ w L (n•)] -γ n [φ K,L ]| = 0
Proof. To prove the above convergence, observe that for any given K, L ∈ Z + , the quantity

|γ n [φ K,L w k (n•) ⊗ w L (n•)] -γ n [φ K,L ]| = p,q∈Z + a 2 n,p,q p n 2k q n 2L (w K (p)w L (q) -1)
converges to zero for large n ∈ N.

Let us consider the term

p 2k q 2L (w K (p)w L (q) -1)
It vanishes if either p ≥ 1 and q ≥ 1 or p = 0 and K ≥ 1 or q = 0 and L ≥ 1. Thus it can be non zero only in the two following (non-disjoint) alternatives:

• If p = 0 and K = 0, then its value is

q 2L (2w L (q) -1) =      q 2L , if q ≥ 1 0 , if q = 0 and L ≥ 1 3
, if q = 0 and L = 0

• If q = 0 and L = 0, then its value is

p 2K (2w K (p) -1) =      p 2K , if p ≥ 1 0 , if p = 0 and K ≥ 1 3
, if p = 0 and K = 0

It follows that for any n ∈ N,

|γ n [φ K,L w k (n•) ⊗ w L (n•)] -γ n [φ K,L ]| = 3a 2 n,0,0 1 {(0,0)} (K, L) + q∈N a 2 n,0,q q n 2L 1 {0} (K) + p∈N a 2 n,p,0 p n 2K 1 {0} (L) = 3a 2 n,0,0 1 {(0,0)} (K, L) + q∈N a 2 n,p,0 p n 2K 1 {0} (L) + p n 2L 1 {0} (K) ≤ 3a 2 n,0,0 + 2 q∈N a 2 n,p,0 ≤ 3 p∈Z + a 2 n,p,0
(where the the symmetry in (p, q) of a n,p,q was taken in account in the second equality). Thus it is sufficient to show that lim

n→∞ p∈Z + a 2 n,p,0 = 0 (17)
From the definition of the coefficients (a n,p,q ) n∈Z + ,p,q∈ 0,n we get for any given n ∈ Z + and x ∈ [0, π], The change of variable formula implies that for any measurable function f ∈ L 1 (ℓ I ) (where ℓ I is the restriction of the Lebesgue measure on I),

1 2π [0,2π] T n cos(x) +
1 2 [0,π] f (θ(y))|θ ′ (y)| dy = I f (z) dz (19)
Taking into account that for y ∈ [0, π], we have

| sin(y)| = 1 -cos 2 (y) = 1 -(2 cos(θ(y)) -cos(x)) 2
we can write (19) as

1 2 [0,π] f (θ(y)) 1 -(2 cos(θ(y)) -cos(x)) 2 | sin(θ(y))| dy = I f (z) dz
Define the function g on I via

∀ z ∈ I, g(z) f (z) 1 -(2 cos(z) -cos(x)) 2 | sin(z)| to obtain 1 2 [0,π] g(θ(y)) dy = I g(z) | sin(z)| 1 -(2 cos(z) -cos(x)) 2 dz equality which is valid as soon as g • θ ∈ L 1 (ℓ [0,π] ).
In particular, taking g = T n , which satisfies this condition, we get

[0,π] T n (θ(y)) dy = 2 I T n (z) | sin(z)| 1 -(2 cos(z) -cos(x)) 2 dz namely [0,π] T n cos(x) + cos(y) 2 dy = 2 I cos(nz) | sin(z)| 1 -(2 cos(z) -cos(x)) 2 dz
Considering this identity with n = 0, we get that the mapping T n cos(x) + cos(y) 2 dy = 0

I ∋ z → | sin(z)| 1 -(2 cos(z) -cos(x)) 2
In addition, we have

∀ n ∈ Z + , ∀ x ∈ [0, π], [0,π] T n cos(x) + cos(y) 2 dy ≤ π
It follows that ( 18) is a consequence of the dominated convergence theorem. ■

Negligible boundary region

For technical reasons that appear later, we wish to limit the integral in Proposition 1 to a subdomain A ε ⊂ [0, 2π] 2 , which for any ε ∈ (0, 1) is defined by

A ε (x, y) ∈ [0, 2π] 2 : cos(x) + cos(y) 2 ≤ 1 -ε .
This is justified since we have lim

ε→0 + lim n→∞ γ n [φ K,L ] - 1 π 2 n 2(K+L) Aε ∂ K x ∂ L y h(x, y) 2 dxdy = 0.
This is a consequence of the following lemma, which shows that the contribution of the boundary region [0, 2π] 2 \A ε is negligible.

Lemma 3. For ε ∈ (0, 1), denote B ε [0, 2π] 2 \ A ε . We have lim ε→0 + sup n≥1 1 n 2(K+L) Bε   p,q∈ 0,n a n,p,q p K q L ξ K (px)ξ L (qy)   2 dxdy = 0 Proof. Note that B ε ⊂ B - ε ∪ B + ε with B - ε {(x, y) ∈ [0, 2π] 2 : cos(x) ≤ -1 + 2ε, cos(y) ≤ -1 + 2ε} B + ε {(x, y) ∈ [0, 2π] 2 : cos(x) ≥ 1 -2ε, cos(y) ≥ 1 -2ε}
We will prove the statement for the B + ε region, since the proof for B - ε is analogous. Define η as the unique solution in (0, π) of cos(η) = 1 -2ε, we have that η goes to 0 + when ε goes to 0 + , and we can rewrite

B + ε   p,q∈ 0,n a n,p,q p K q L ξ K (px)ξ L (qy)   2 dxdy = [-η,η] 2   p,q∈ 0,n a n,p,q p K q L ξ K (px)ξ L (qy)   2 dxdy = p,q,p ′ ,q ′ ∈ 0,n a n,p,q a n,p ′ ,q ′ p K q L (p ′ ) K (q ′ ) L η -η ξ K (px)ξ K (p ′ x) dx η -η ξ L (qy))ξ L (q ′ y) dy
Let us compute the integrals of the r.h.s. When K is even, we have

η -η ξ K (px)ξ K (p ′ x) dx = 2 η 0 cos(px) cos(p ′ x) dx = η 0 cos((p + p ′ )x) + cos((p -p ′ )x) dx = - sin((p + p ′ )η) p + p ′ - sin((p -p ′ )η) p -p ′
with the usual convention that sin((p-p ′ )η)

p-p ′ = η when p -p ′ = 0. When K is odd, a similar computation shows that η -η ξ K (px)ξ K (p ′ x) dx = sin((p + p ′ )η) p + p ′ - sin((p -p ′ )η) p -p ′ .
It follows that the main integral is a sum of terms of the form p,q,p ′ ,q ′ ∈ 0,n ±a n,p,q a n,p ′ ,q

′ p K q L (p ′ ) K (q ′ ) L sin((p ± p ′ )η) p ± p ′ sin((q ± q ′ )η) q ± q ′ ,
for different choices of the signs. We consider the case with all minus signs, as the other cases follow by essentially the same argument. By a triangle inequality and the Cauchy-Schwarz inequality we can bound the sum by p,q,p ′ ,q ′ a n,p,q a n,p ′ ,q ′ sin((p

-p ′ )η) p -p ′ sin((q -q ′ )η) q -q ′ ≤ p,q,p ′ ,q ′ a 2 n,p,q sin((p -p ′ )η) p -p ′ sin((q -q ′ )η) q -q ′ 2 p,q,p ′ ,q ′ a 2 n,p ′ ,q ′ sin((p -p ′ )η) p -p ′ sin((q -q ′ )η) q -q ′ 2 = p,q,p ′ ,q ′ a 2 n,p,q sin((p -p ′ )η) p -p ′ sin((q -q ′ )η) q -q ′ 2 ≤ p,q a 2 n,p,q r,s∈Z sin(rη) r sin(sη) s 2 = p,q a 2 n,p,q   r∈Z sin(rη) r 2   2 14 = p,q a 2 n,p,q   η 2 + 2 r∈N sin(rη) r 2   2 ≤ 4   η 2 + 2 r∈N sin(rη) r 2   2
where we used Corollary 1. Dominated convergence now implies that lim

η→0 + r∈N sin(rη) r 2 = 0,
which proves the claimed convergence. ■

Simplifying the integral

We will now approximate the integral

Aε ∂ K x ∂ L y h(x, y) 2
dxdy by a (arguably) simpler integral.

First simplification

We use Faà di Bruno's formula to rewrite

∂ K x h(x, y) = ∂ K x T n ((cos(x) + cos(y))/2) as m 1 +2m 2 +•••+Km K =K K! m 1 !m 2 ! • • • m K ! T ( m ℓ ) n cos(x) + cos(y) 2 k∈ K ξ k (x) 2k! m k
, where the sum is over all nonnegative integers (m 1 , ..., m K ) satisfying the constraint

m 1 + 2m 2 + • • • + Km K = K. Applying Faà di Bruno's formula again we can rewrite ∂ K x ∂ L y h(x, y) as m 1 +2m 2 +•••+Km K =K n 1 +2n 2 +•••+Ln L =L K! m 1 !m 2 ! • • • m K ! L! n 1 !n 2 ! • • • n L ! T ( m ℓ + n ℓ ) n cos(x) + cos(y) 2 k∈ K ξ k (x) 2k! m k l∈ L ξ l (y) 2l! n l (20)
How terrible this expression may seem, as n → ∞ the final integral will be dominated by a single term in the summand. To see this, we use the following bound.

Lemma 4. For any ε ∈ (0, 1) and r ∈ N, there exists a constant C(ε, r) > 0 such that

∀ n ∈ Z + , ∀ x ∈ [-1 + ε, 1 -ε], |T (r) n (x)| ≤ C(ε, r)n r
Proof. Recall that T n (x) = cos(n arccos(x)) for all x ∈ [-1, 1]. The desired result follows from the fact that any given derivative of arccos is bounded on [-1 + ε, 1 -ε]. ■

It follows that generically in (20), as n → ∞, the main term corresponds to m 1 = K, n 1 = L, and all other coefficients equal to zero. We can prove the following lemma, where we define the integral

I K,L (ε, n) 1 π 2 n 2(K+L) Aε T (K+L) n cos(x) + cos(y) 2 -sin(x) 2 K -sin(y) 2 L 2
dxdy.

Lemma 5. For any ε > 0 we have

lim n→∞ 1 π 2 n 2(K+L) Aε ∂ K x ∂ L y h(x, y) 2 dxdy -I K,L (ε, n) = 0.
Proof. Using the expression (20) and expanding the square, up to some factor not depending on n, we end up with products of the form

T ( m ℓ + n ℓ ) n cos(x) + cos(y) 2 T ( m ′ ℓ + n ′ ℓ ) n cos(x) + cos(y) 2 (21)
where K ℓ=1 ℓm ℓ = K ℓ=1 ℓm ′ ℓ = K and L ℓ=1 ℓn ℓ = L ℓ=1 ℓn ′ ℓ = L. According to Lemma 4, up to a constant depending on ε, K and L, the quantity (21) is bounded independently for (x, y) ∈ A ε by

n m ℓ +m ′ ℓ + n ℓ +n ′ ℓ
which is negligible with respect to n 2(K+L) unless the only nonzero coefficients are m 1 = m ′ 1 = K and n 1 = n ′ 1 = L. The only potentially non-vanishing term in the integral is hence

T (K+L) n cos(x) + cos(y) 2 -sin(x) 2 K -sin(y) 2 L 2
and this leads to the announced convergence. ■

Second simplification

We do a second simplification by approximating the integral I K,L (ε, n). Applying Faà di Bruno's formula to the composition of T n with cosine, we get for any θ ∈ [0, π],

∂ K+L θ T n (cos(θ)) = m 1 +2m 2 +•••+(K+L)m K+L =K+L (K + L)! m 1 !m 2 ! • • • m K+L ! T (m 1 +m 2 +•••+m K+L ) n (cos(θ)) k∈ K+L ξ k (θ) k! m k
By a similar argument as in Lemma 5, we have that when θ in (0, π) and K, L ∈ Z + are fixed then

lim n→∞ n -(K+L) ∂ K+L θ T n (cos(θ)) -T (K+L) n cos(θ) (-sin(θ)) K+L = 0 (22)
Alternatively, taking into account that T n (cos(θ)) = cos(nθ), we have that

n -(K+L) ∂ K+L θ T n (cos(θ)) = n -(K+l) ∂ K+L θ cos(nθ) = ξ K+L (nθ)
Using some uniformity in the behavior (22) with respect to θ away from 0 and π, as in the proof of Lemma 5, we deduce that for any fixed ε ∈ (0, 1), lim

In particular we have lim

ε→0 + lim sup n→∞ |I K,L (ε, n) -J K,L (ε, n)| = 0 with J K,L (ε, n) 1 π 2 2 2(K+L) Aε ξ 2 K+L (nθ(x, y)) sin(x) sin(θ(x, y)) 2K sin(y) sin(θ(x, y)) 2L dxdy = 4 π 2 2 2(K+L) Aε∩[0,π] 2 ξ 2 K+L (nθ(x, y))
sin(x) sin(θ(x, y)) 2K sin(y) sin(θ(x, y)) 2L dxdy where symmetries x → -x and y → -y were taken into account.

Combining all of the above claims, we get that lim

ε→0 + lim n→∞ γ n [φ K,L ] -J K,L (ε, n) = 0.

Integral as a second moment

In the final part of the proof we show that, as n → ∞, the quantity J K,L (ε, n) approximates the second moment of some continuous measure γ. More specifically, we prove the following lemma.

Lemma 6. It holds that

lim ε→0 + lim sup n→∞ J K,L (ε, n) = 1 π 2 [0,π] 2 sin(x) sin(θ(x, y)) 2K sin(y) sin(θ(x, y)) 2L dxdy.
It follows that we can define the measure γ as the image of the measure 1 π 2 dxdy on [0, π] 2 by the mapping Ψ defined by Ψ : [0, π] 2 ∋ (x, y) → sin(x) sin(θ(x, y)) , sin(y) sin(θ(x, y)) .

In particular it appears that γ is a probability measure, so that µ = γ, as announced previously.

Combining the previous arguments, we see that lim n→∞ γ n [ϕ K,L ] = γ[ϕ K,L ], as wanted.

Proof of Lemma 6

Let us rewrite the integral J K,L (ε, n) using the change of variables

F : [0, π] 2 ∋ (x, y) → (x, θ(x, y)).
This is a bijection on its image, which is the set

Ω {(x, θ) ∈ [0, π] 2 : 2 cos(θ) -cos(x) ∈ [-1, 1]}. ( 25
)
From the change of variables formula (see Lemma 7 at the end of this proof), we have that for any measurable f :

Ω → R + , [0,π] 2 f (F (x, y)) dxdy = Ω f (x, θ) 2 sin(θ) 1 -(2 cos(θ) -cos(x)) 2 dxdθ. (26)
We can write the integral J K,L (ε, n) in this form by noting that

π 2 2 2(K+L) J K,L (ε, n) = | cos(θ(x,y))|≤1-ε ξ 2 K+L (nθ(x, y)) sin(x) sin(θ(x, y)) 2K 1 -cos 2 (y) L sin 2L (θ(x, y)) dxdy = | cos(θ(x,y))|≤1-ε ξ 2 K+L (nθ(x, y)) sin(x) sin(θ(x, y)) 2K 1 -(2 cos(θ(x, y)) -cos(x)) 2 L sin 2L (θ(x, y)) dxdy = [0,π] 2 f n (F (x, y)) dxdy where ∀ (x, θ) ∈ Ω, f n (x, θ) 1 [0,1-ε] (| cos(θ)|)ξ 2 K+L (nθ) sin(x) sin(θ) 2K 1 -(2 cos(θ) -cos(x)) 2 L sin 2L (θ)
The change of variables (26) implies

π 2 2 2(K+L) J K,L (ε, n) = Ω f n (x, θ) 2 sin(θ) 1 -(2 cos(θ) -cos(x)) 2 dxdθ Note that ∀ r ∈ Z + , ∀ z ∈ [0, π], ξ 2 r (z)    1+cos(2z) 2 , if r ≡ 0[2] 1-cos(2z) 2 , if r ≡ 1[2]
Consider the function g defined on Ω via

∀ (x, θ) ∈ Ω, g(x, θ) sin(x) sin(θ) 2K 1 -(2 cos(θ) -cos(x)) 2 L sin 2L (θ) 2 sin(θ) 1 -(2 cos(θ) -cos(x)) 2 so that 2π 2 2 2(K+L) J K,L (ε, n) = Ω 1 + (-1) K+L cos(2nθ) g(x, θ)1 [0,1-ε] (cos(θ)) dxdθ
For ε ∈ (0, 1) still fixed, the mapping Ω ∋ (x, θ) → g(x, θ)1 [0,1-ε] (cos(θ)) is integrable with respect to the Lebesgue measure on Ω. We can apply the Riemann-Lebesgue theorem to get

2π 2 2 2(K+L) lim n→∞ J K,L (ε, n) = Ω g(x, θ)1 [0,1-ε] (cos(θ)) dxdθ,
and by monotone convergence we have lim

ε→0 + Ω g(x, θ)1 [0,1-ε] (cos(θ)) dxdθ = Ω g(x, θ) dxdθ.
Finally, note that the previous change of variable computations can be reversed, and yields

1 2π 2 2 2(K+L) Ω g(x, θ) dxdθ = 1 π 2 [0,π] 2 sin(x) sin(θ(x, y)) 2K sin(y) sin(θ(x, y)) 2L dxdy.
This proves Lemma 6. ■

The following result gives the details for the change of variables formula used in Formula (26).

Lemma 7. We have for any measurable f :

Ω → R + , [0,π] 2 f (F (x, y)) dxdy = Ω f (x, θ) 2 sin(θ) 1 -(2 cos(θ) -cos(x)) 2 dxdθ.
Proof. In addition to the set Ω defined in (25), for given θ ∈ [0, π], denote 6 Discussion and future work

We discussed how quantum walks on graphs can be interpreted as solutions of a discrete wave equation on the corresponding graph. This naturally gives rise to a description of the quantum walk dynamics using Chebyshev polynomials. We then used the analysis of Chebyshev polynomials to prove a weak limit on the spreading behavior of quantum walks on lattices.

The main open direction of this work is to use Chebyshev polynomials for studying quantum walk spreading and mixing on more general graphs. The analysis that we did for the lattice seems to naturally extend to lattices of higher dimensions, or Cayley graphs of Abelian groups more generally. An interesting example out of this category, while still bearing some similarity with the lattice, is the Cayley graph of the Heisenberg group (as studied e.g. in [BDH + 17]). Preliminary simulations indicate a faster spreading rate of quantum walks, yet we are not aware of a proof of this.

Exploring the connection with the wave equation on graphs by Friedman and Tillich [START_REF] Friedman | Wave equations for graphs and the edge-based Laplacian[END_REF] is another open direction. These authors associate to a Markov generator △ G (so that I -△ G is a transition matrix, where I is the identity matrix) on a discrete graph G a diffusion generator △ G on the geometric realization G of G, where the edges have been replaced by segments of length 1 on which usual one-dimensional Laplacian operators are acting. The domain of △ G consists of functions satisfying some conditions at the vertices. Friedman and Tillich then relate the eigenvalues of △ G and △ G , a feature which enables them to interpret T k (I -△ G ) as an operator related to the evolution of the wave equation associated to △ G (see their Theorem 4.9). We leave it as an open question whether this interpretation of the action of the Chebyshev polynomials on a transition matrix can be used to deduce our weak convergence result.

Figure 1 :

 1 Figure 1: (left) Probability measure induced by [P n L ] 0,• for n = 50, corresponding to an n-step lazy random walk on the line. (right) Probability measure induced by [T n (P L )] 2 0,• for n = 50, associated to an n-step quantum walk on the line. (left and right) Lines drawn continuously for clarity.

Figure 2 :

 2 Figure 2: (left) Probability measure induced by [P n ] 0,• for n = 50, corresponding to an n-step random walk on the lattice. (right) Probability measure induced by [T n (P )] 2 0,• for n = 50, associated to an n-step quantum walk on the lattice.

  x ∈ [0, π]. To investigate the integral [0,π] T n cos(x)+cos(y) 2 dy we consider the change of variable y → θ x (y) where θ x (y) is the unique solution belonging to [0, π] of the equation cos(θ x (y)) = cos(x) + cos(y) 2To simplify notations, from now on we remove the index x. The mapping θ is a bijection from [0, π] to

  belongs to L 1 (ℓ I ). Thus we can apply Riemann-Lebesgue theorem to deduce lim

  n→∞ Aε ξ K+L (nθ(x, y))(-sin(θ(x, y))) K+L 2 -n -2(K+L) T (K+L)n cos(θ(x, y)) 2 sin 2K (x) sin 2L (y) dxdy = 0 with θ(x, y) the unique solution in [0, π] of cos(θ(x, y)) = cos(x) + cos(y) 2 .

  Ω θ{x ∈ [0, π] : (x, θ) ∈ Ω} Thus for any measurable and non-negative mapping g : Ω → R + , we have by the change of variables formula[0,π] 2 g(F (x, y))|det(Jac(F ))|(x, y) dxdy = Ω g(x, θ) dxdθ(27)where Jac(F ) is the Jacobian matrix associated to F . We haveJac(F (x, y)) = 1 ∂ x θ(x, y) 0 ∂ y θ(x, y)Differentiating (23) with respect to y, we get, when θ(x, y) ̸ ∈ {0, π} (i.e. (x, y) ̸ ∈ {(0, 0), (π, π)}),∂ y θ(x, y) = sin(y) 2 sin(θ(x, y)) so that |det(Jac(F ))|(x, y) = |∂ y θ(x, y)| = | sin(y)| 2| sin(θ(x, y))| = 1 -cos 2 (y) 2 sin(θ(x, y)) = 1 -(2 cos(θ(x, y)) -cos(x)) 2 2 sin(θ(x, y))Formula (27) rewrites[0,π] 2 g(F (x, y)) 1 -(2 cos(θ(x, y)) -cos(x)) 2 2 sin(θ(x, (θ) -cos(x)) 2 dxdθwhere g was replaced by the non-negative and measurable mapping f given by∀ (x, θ) ∈ Ω \ {(0, 0), (π, π)}, f (x, θ) g(x, θ) 1 -(2 cos(θ) -cos(x)) 2 2 sin(θ) ■ Remark 1.The fact that γ is a probability implies (in fact, is equivalent to) the convergence lim n→∞ p,q∈ 0,n a 2 n,p,q = γ(1) = 1(compare with Corollary 1), which also amounts to lim n→∞ p,q∈ 1,n a 2 n,p,q = 1 according to (17).

The connection between quantum walks and Chebyshev polynomials has been used before (see e.g.[START_REF] Childs | Quantum algorithm for systems of linear equations with exponentially improved dependence on precision[END_REF][START_REF] Gilyén | Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics[END_REF]) but not in the context of quantum walk mixing.

There are alternative options, such as u(n+ 1) = u(n) + v(n), v(n + 1) = v(n) -(I -P 2 )u(n).

Acknowledgements

Simon Apers is partially supported by French projects EPIQ (ANR-22-PETQ-0007), QUDATA (ANR-18-CE47-0010) and QUOPS (ANR22-CE47-0003-01), and EU project QOPT (QuantERA ERA-NET Cofund 2022-25).


Clément Boudaud

Luc Davenne

Julien Baroth

Laurent Daudeville

MULTI-SCALE MODELING OF TIMBER-FRAME STRUCTURES UNDER SEISMIC LOADING

Keywords: Timber-frame structures, Finite element, Cyclic and dynamic loading, Constitutive model, Macro-scale

. It is now generalized to a large class of joints and validated at each scale, thanks to an important number of tests (more than 400 on joints, 25 on shear walls and 12 on roofs). Moreover, an improvement of the 1D hysteretic constitutive behaviour [4] is proposed to model dissymmetric behaviors and to model cumulative damages appearing under alternate loading. Eventually, a shake table test is performed on a 6 x 6 meters timber-frame house, and the results are confronted to the predictions of the numerical model of structure.

INTRODUCTION 123

The study presented in this paper is motivated by two facts. Firstly, timber-frame construction is becoming a common building system in Europe. It presents many qualities, one of which being a good earthquake resistance due to both the excellent strength-to-density ratio of timber and the ductility of joints with metal fasteners, leading to limited inertia forces and providing good energy dissipation. Secondly, the recent European code for design of earthquake resistant buildings is accompanied in some countries by a new seismic hazard map. In France, based on this revised map, earthquake resistance calculations are now mandatory in a much larger part of the territory. Therefore, the seismic behavior of timber-frame structures has to be studied, in order to better understand their global and local behaviors.

Because nonlinear dissipative phenomena in timber-frame structures are mainly concentrated into metal connectors, simplified force-displacement models can be derived by fitting results of tests performed on joints. The proposed approach is based on a multi-scale concept, as proposed 1 LUNAM Univ., Groupe Ecole Supérieure du Bois, Nantes, France. Email: clement.boudaud@ecoledubois.fr 2 Université Paris Ouest Nanterre la Défense, LEME, Ville d'Avray, France. 3 Univ. Grenoble Alpes, 3SR, F-38000 Grenoble, France. CNRS, 3SR, F-38000, France. Emails: firstname.name@3sr-grenoble.fr previously by Richard et al. [START_REF] Richard | Timber shear walls with large openings: experimental and numerical prediction of the structural behavior[END_REF] or Xu et al. [START_REF] Xu | Development of a wood-frame shear wall in Abaqus[END_REF]. Such an approach requires a behavior law to represent the forcedisplacement evolution at each scale.

Figure 1:

One-dimensional constitutive model used for the description of the force-displacement evolution at the scale of elementary joints or shear walls [START_REF] Humbert | Joints and wood shear walls modeling I: Constitutive law, experimental tests and FE model under quasi-static loading[END_REF] Numerous constitutive laws have been developed over the years (Ceccotti et al. [START_REF] Bertil | Modélisation de l'action sismique[END_REF], Folz et al. [START_REF] Ceccotti | A pinching hysteretic model for semirigid joints[END_REF]). A new model, developed by Humbert et al. [START_REF] Humbert | Joints and wood shear walls modeling I: Constitutive law, experimental tests and FE model under quasi-static loading[END_REF], can be considered as an improvement of Richard et al. [START_REF] Folz | Seismic analysis of wood-frame structures II: Model implementation and verification[END_REF] and Yasumura [START_REF] Richard | Prediction of the seismic behavior of wood-framed shear walls with openings by pseudodynamic test and FE model[END_REF] models. It fulfils the needs of modeling asymmetric behaviors and cumulative damages appearing under alternate loading, and also provides a strict numerical continuity by using Bézier polynomials instead of commonly used exponential functions. This onedimensional constitutive model is shown in Figure 1 and its parameters are listed in tables 1, 2 and 3.

Experimental tests

More than 400 experimental tests on joints with metal fasteners were performed to provide input data for the numerical models of structural elements. Tests were performed under monotonic and cyclic loading following the protocol of the EN 12512 [START_REF] Yasumura | Evaluation of the damping capacity of timber structures for seismic design[END_REF]. Tests were also repeated two or three times for monotonic loadings and five times 2.2 SCALE 2: STRUCTURAL ELEMENTS

Experimental tests

Several configurations of common 2.4 x 2.4 meters shear walls (Figure 4) were used to obtain the experimental behavior of shear walls under quasi-static (monotonic and cyclic) and dynamic (shake table) loadings, for a total of respectively 14 and 11 tests. Figure 6 shows a roof installed on the shaking table. As the table is only unidirectional, two tests were performed on each configuration of roof, with the roof trusses respectively parallel and perpendicular to the loading direction.

Numerical modeling

A detailed FE model of shear walls is developed using beam elements for the frame, plate elements for the panels and two-node spring-like finite elements for the joints. Figure 7 displays an exploded view of the mesh for a better understanding. The behaviors of joints have been identified at scale 1, k n , k s , k f represent the non-linear stiffness of the joints. A detailed FE model of the roof is developed using the same approach. A refined model of the roof was developed following the same modelling principles as the refined models of the shear walls.

Then, the shake table tests performed on the shear walls (11 tests) and the roofs (12 tests) were simulated. Figure 9 presents the comparison between experimental and numerical max and min relative displacements of the shear walls for the 11 dynamic tests. Numerical predictions are more accurate for some tests than others, but on the whole the model is able to predict the behavior for different configurations of shear walls and different levels of loading. Similar work was achieved for roofs.

Figure 9: Comparison of experimental results and numerical predictions for a shear wall with 9 mm OSB panels [4]

Note that the modeling of the dynamical tests performed on shear walls and roofs was carried out by considering a viscous damping. A Rayleigh damping was considered i.e. the damping matrix is as a linear combination of the mass and stiffness matrices. The damping matrix is built by considering the two first modes. The damping ratio is identified using low amplitude free vibration responses after white noise solicitations. It is assumed that for relative displacement of the order of 1 mm, hysteretic damping does not occur. Based on the conclusions of Dutil and Symans [START_REF] Dutil | Experimental investigation of seismic behavior of light-framed wood shear walls with supplemental energy dissipation[END_REF], the logarithmic decrement method is used to calculate the damping ratio, rather than the bandwidth method based on the Frequency Response Function (FRF). For shear walls and roofs, the two identified ratios were 5 % and 2 % respectively.

SCALE 3: STRUCTURE

Experimental tests

A dynamic test has been performed on a single-story 6 x 6 meters timber-frame house (Figure 10). All three translational DOF of the shake table are submitted to ground motions corresponding to a severe seismic hazard. The scenario for Guadeloupe (French Caribbean Island) has been chosen. The most probable magnitude-distance couple was identified for this scenario and for a return period of 475 years (probability of exceedance of the peak ground acceleration of 10 % in 50 years). Seismic ground motion histories are represented Figure 11. The resulting peak ground acceleration is 0.33 g in Y-axis. These histories have been applied to the structure, considering a growing amplification coefficient. Resulting PGA are between 0.33 g (100%) to 0.99 g (300%).

Figure 10: Overview of the house on the shaking table

Experimental observations show that there is no torsion of the structure during the test. It was expected that the torsion would be small but observable, as one of the walls includes a wide opening of 2.8 m for a wall length of 6 m. This result is most certainly explained by the reinforcement of the roof, in order to behave like rigid diaphragm, which was in accordance with the French building regulations.

Numerical modeling

A detailed FE modeling of the whole timber-frame house is too time-consuming. Because of the particular parallelogram-like kinematics of the shear wall, it is proposed to consider the shear wall deformation by means of a simplified FE model consisting of a one degree of freedom moving frame (Figure 12). Shear wall kinematic is a combination of parallelogram-like deformation, overturning (due to uplift of the anchorage) and shear (Gupta and Kuo [START_REF] Gupta | Behavior of wood-framed shear walls[END_REF]). The simplified FE model only has the parallelogram-like kinematic. Its nonlinear behaviour of the simplified FE model arises from the constitutive relationship used for the joints [START_REF] Humbert | Joints and wood shear walls modeling I: Constitutive law, experimental tests and FE model under quasi-static loading[END_REF]. The parameters of this non-linear model were identified from the results of the refined FE model under both push-over and cyclic loadings. Therefore, the effect of overturning and shear are indirectly taken into account.

Figure 12: Simplified FE model of shear wall

The simplified FE model is calibrated for quasi-static loading and then used in dynamic calculations. To assess the accuracy of its behaviour under dynamic loading, the simplified FE model predictions and the experimental results under dynamic loading are compared. The results show a good agreement, validating the use of a single degree of freedom model using the developed constitutive behaviour law.

Figure 13: Overview of the finite element model of the house

The model of the structure is then built by assembling the simplified FE models of shear walls with kinematic relations, and eventually the detailed FE model of roof (Figure 13). The roof was attached to the simplified shear Y X Z walls by spring-like elements that model the 3D brackettype connectors.

In the following, only the nominal ground motion presented in Figure 11 is applied with 100 % and 200 % amplification coefficients (PGA respectively equal to 0.33 g and 0.66 g).

Results

One should note that the FE model predicts relative displacements at the top of the shear wall with the wide opening more important than those of the opposite wall. This induces a torsional effect on the structure that is not observed experimentally. The cause of that error between the experimentation and the model are currently ongoing research. Two possible causes are that the reinforcement of the roof might not be rigid enough, and that the wall to wall connections in the angles should transfer a resistive moment (this DOF was free for these calculations). The following comparison between experimental results and numerical predictions is made at the top of the shear wall opposite to the one with the wide opening. Figure 15 presents the model predictions in terms of peak displacements for the same input ground motions, considering that 0.33 and 0.66 g ground motions are applied consecutively in the same calculation (meaning that the 0.66 g solicitation is applied on the already damaged model after the 0.33 g loading).

Experimental and numerical histories are rather similar, with the same maximum displacements at same times: respectively around 4 and 10 mm for 0.33 and 0.66 g ground motions.

Discussion

The presented results were obtained with a first blind simulation. Further analysis will be carried out.

The slight underestimation of displacements (Table 4) may be due to an overestimated damping (identified for both shear walls and roofs by means of shake table tests).

CONCLUSIONS

This study presents a numerical model of timber frame structure able to predict the nonlinear behavior under dynamic loadings. Using a multi-scale approach, from the elementary metallic joints to the structure, provides a computationally efficient and yet relatively accurate tool. More than 400 tests provide a large database of connectors' behavior and detailed FE models of structural elements are confronted to 14 quasi-static and 11 dynamic tests for shear walls and 12 dynamic tests for roofs. Ongoing research consists in improving numerical predictions at the scale of the house. In the future, it is intended to use this model as a tool to better understand the behavior of timber frame structures under seismic loading. It will also be possible to compare the results of a timehistory analysis to the simplified design methods of the Eurocode 8 [START_REF]EN 1998-1: Design of structures for earthquake resistance -General rules, seismic actions and rules for buildings[END_REF].

Figure 2 : 2 . 1 . 2

 2212 Figure 2: Experimental tests on metal fastened joints

Figure 3 :

 3 Figure 3: Direct calibration of a nailed joint (sheathing to framing connection)[START_REF] Humbert | Joints and wood shear walls modeling I: Constitutive law, experimental tests and FE model under quasi-static loading[END_REF]

Figure 4 :

 4 Figure 4: Shear wall description

Figure 5 :

 5 Figure 5: Shear wall setup for dynamic test

Figure 6 :

 6 Figure 6: Roof setup for dynamic test

Figure 7 :

 7 Figure 7: Finite element modeling of a shear wall

Figure 8 :

 8 Figure 8: Comparison of experimental results and numerical predictions for a shear wall with 9 mm OSB panels[START_REF] Humbert | Joints and wood shear walls modeling I: Constitutive law, experimental tests and FE model under quasi-static loading[END_REF]

Figure 11 :

 11 Figure11: Ground motion histories corresponding to a severe seismic hazard (Guadeloupe, French Caribbean island)[START_REF] Bertil | Modélisation de l'action sismique[END_REF]

FiguresFigure 14 :

 14 Figures 14 depicts experimental top wall displacements for input ground motions with respectively 0.33 g (a) and 0.66 g (b) maximum peak ground accelerations (PGA).

Figure 15 :

 15 Figure 15: Evolutions of top wall displacements from FE simulations for 0.33 g and 0.66 g seismic loads

Table 1 :

 1 Model parameters governing the constitutive behaviour under monotonic loading

	Parameter	Unit	Description
	K 0	N/m	Initial elastic stiffness
	d y	m	Yield limit
	d 1	m	Displacement at peak force
	F 1	N	Peak force
	K 1	N/m	Pre-peak tangent stiffness
	d 2	m	Intermediate displacement limit
	F 2	N	Force at intermediate limit d 2
	d u	m	Ultimate displacement
	F u	N	Force at ultimate displacement

Table 2 :

 2 Model parameters governing the shape of the hysteresis loops

	Parameter	Unit	Description
	C 1	-	Unloading stiffness
	C 2	-	Reloading stiffness
	C 3	-	Tangent stiffness at F=0
	C 4	-	Residual displacement

Table 3 : Model parameters governing the damage indicator calculation

 3

	Parameter	Unit	Description
	B C	-	Linear coefficient of the DLF
	B R	-	Power term of the DLF
	η	%	Damage proportion at constant amplitude cycles
	DLF: Damage Limit Function

2 COUPLED EXPERIMENTAL- NUMERICAL APPROACH 2.1 SCALE 1: JOINTS

Table 4

 4

presents relative displacements, between the top of walls and roofs and base of the structure. Computed and experimental values are compared. In this case, experimental values are slightly under-estimated by the numerical model.

Table 4 :

 4 Computed and experimental average displacements at the top of walls and roofs

	Peak ground	0.33 g	0.66 g
	acceleration	num exp num exp
	Wall/base (mm)	3.3 4.5 8.5	10
	Roof/wall (mm)	2.8	4	6	6.5

ACKNOWLEDGEMENT

The French Agency (ANR RiskNat 2008) is gratefully acknowledged. The authors would like to express their gratitude to all the members of the SISBAT research project, especially C. Faye (FCBA Technological Institut), P.E. Charbonnel, M. Chaudat (CEA) and their colleagues for the effort provided for the realization of the experimental campaign.


Yannick Le Cozler
email: yannick.lecozler@institut-agro.fr

Maxime Dumesny

Jean-Michel Delouard

Laurent Delattre

Luginbuhl Thibault

Philippe Faverdin

Philippe †

Faverdin

One-shot technology for 3D imaging of large animals: perspectives for ruminant management

Keywords: cattle, depth camera, one shot, 3D images

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Early attempts to monitor changes in livestock animals using imaging technologies were initially based on two-dimensional (2D) imaging (Schofield et al., 1999). Two D camera-based approaches gave significant results on estimating body area and /or other morphological traits (Tasdemir et al., 2011;Ozykaya et al., 2016;Yan et al., 2019, Weber et al., 2020a,b). However, the process of converting a three-dimensional object into a 2D image encountered significant challenges associated with image distortion and deformation. Research focus thus shifted to threedimensional (3D) imaging technologies (Gauthier, 2005), which are, like 2D imaging technology, non-intrusive and capture partially or totally the external envelope of a living animal, depending of the number of cameras used (Le Cozler et al., 2019a;Wang et al., 2018;Martins et al., 2020;Li et al., 2022). With this approach, it is possible to describe the morphological characteristics of an animal, or changes in such characteristics over time, without any risk to the handler or the animal, even if the animal is vigorous and/or dangerous. When used with appropriate software, 3D imaging has great potential for quantifying or estimating technical indicators of interest such as body length or body condition score, live weight, or lameness detection, among others [START_REF] Bewley | Potential for estimation of body condition scores in dairy cattle from digital images[END_REF][START_REF] Halachmi | Cow body shape and automation of condition scoring[END_REF]Fischer et al., 2015;Pezzuolo et al., 2018;Le Cozler et al., 2020b;Le Cozler et al., 2022a). Past work by our group has demonstrated the added value of 3D imaging for monitoring morphological traits to which it was rarely possible to have direct access, such as surface area and volume. From these indicators, it then becomes possible to estimate body weight in both adults and young animals (Le Cozler et al., 2019a;Le Cozler et al., 2022b;Wang et al., 2023). We have also shown that information from 3D imaging technology can be used to study the nature of weight gain during lactation (Xavier et al., 2022a), but also to estimate the chemical composition of animals (Xavier et al., 2022b). For these earlier studies, we used an imaging device Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024 A c c e p t e d M a n u s c r i p t called "Morpho3D", which is very accurate but requires specific conditions of use that are not always easy to obtain, including the absence of direct light and a complete lack of movement by the animals while recording images (5-s acquisition time).

To address these problems, and to improve the speed and ease of acquiring a large amount of phenotypic information, Li et al. (2022) developed a real-time point-cloud collection gantry system, which enables body measurements for active cattle on a ranch. Our group, instead, developed an approach based on the use of multiple depth cameras, triggered simultaneously ("one shot"), which can capture the external envelope of the animal via the coupling of the different image captures. With this technique, we hoped to develop a method that was less affected by animal movements during the acquisition time and was less dependent on direct light. To this end, we have created a new imaging device, called "Deffilait3D" that combines multiple simultaneously acquired images to create a 3D reconstruction of an animal. Here, we present the results of validation tests in which we compared values measured on Holstein cows using the Deffilait3D and Morpho3D systems, the measurements from the latter, previously validated, device used as reference (Le Cozler et al., 2019a). In addition to evaluating repeatability and reproducibility, we also assessed the accuracy of earlier equations developed with the Morpho3D device in estimating the body weight of animals (Le Cozler et al., 2019b) using values collected from Deffilait3D.

MATERIALS AND METHODS

The experiment was conducted within the ANR Deffilait project, with procedures approved by the local ethics committee and the French Ministry of Higher Education, Research, and Innovation (reference number APAFIS 3122-2015112718172611).

Experimental design

The study was performed from March to June 2021 at the Méjusseaume experimental station of INRAE, Dairy Nutrition and Physiology unit (IE PL, 35560 Le Rheu, France (https://doi.org/10.15454/yk9q-pf68) under agreement for animal housing no. C-35-275-23). This facility has a herd of approximately 180 Holstein cows and their offspring (e.g., replacement heifers), producing an average of 8,500 kg of milk per cow per year. All animals were reared and managed under routine farm procedures consistent with French animal-welfare regulations.

Data collection

Two devices for the acquisition and processing of 3D images of whole animals are present in this experimental installation, both developed and fine-tuned by the company 3DOuest (Lannion, France):

-Measurements obtained from the Morpho3D device were used as reference values ("gold standard"), since this system has been repeatedly tested, validated, and used since 2017 (Le Cozler et al., 2019a for more details; Figure 1). For this reasons, no manual measurements, as traditionally performed and considered as "gold standards" for morphological traits [START_REF] Heinrichs | Predicting body weight and wither height in Holstein heifers using body measurements[END_REF], were performed. Morpho3D device is a sliding acquisition system with five cameras distributed over the sides and top of the gate, each coupled to a laser projector. As the gate moves back and forth, each camera takes 80 images per second, giving a total of 2,000 images. The images of the laser strips projected onto the animals are captured by the corresponding camera and sent to a computer. During the image acquisition period (5 to 6 s), the animal must stand completely still.

For this reason, four stainless steel cables are also used within the Morpho3D Scanner, on both side of cow, to secure the scanner and restrict cow movement. Cows can also be restrained by a self-locking head gate if necessary but as it becomes accustomed to the device and its movement, A c c e p t e d M a n u s c r i p t restraints are no longer necessary. After processing, a 3D reconstruction is generated, with a unique point cloud representation of the whole animal. The surface normals are estimated from the point cloud and a screened Poisson surface reconstruction algorithm is applied to construct a triangulated mesh [START_REF] Kazhdan | Screened Poisson surface reconstruction[END_REF], using the open source software Meshlab (Cignoni et al., 2008).

-The Deffilait3D device was installed at the end of 2020. This prototype consists of a fixed metal frame, with 15 RGB-D sensors in total (Intel RealSense D415) (Figure 2). The Intel RealSense D415 has a compact design, low cost, and low power usage (USB-powered). It is usually used in indoor environments, as sunlight may interfere with the infrared (IR) intensity of projector, leading to holes in the depth map. The RealSence D415 is depth camera and consists of a pair of depth sensors, an RGB sensor and an infrared spotlight. This technology is a stereo system: it triangulates the depth (distance) for each pixel between the 2 sensors. The infrared projector helps the algorithm to find correspondence between the 2 images. For the animal reconstruction, in a passive stereo system, two RealSence D415 cameras, whose relative position and orientation are known, are used to capture the same scene. The two images are searched for corresponding pixels, and for each correspondence found, the depth can be estimated using camera parameters and triangulation.

Here, the active stereo is used and aims to simplify the search for corresponding pixels by adding data with a light source emitting a pattern. Depth information is represented using false colors images, blue and red indicate the distance from far (blue) to near (red). Each sensor is connected to a computer, which is itself connected to a terminal. The images taken simultaneously by these 15 devices are then combined to reconstruct a 3D image of the entire animal. A common coordinate frame is chosen (the upper camera frame) and coordinate frame transforms for every acquisition are established through a calibration procedure: a checkerboard pattern is shown to consecutive camera' pairs, the relative positions of the cameras are estimated using this pattern and OpenCV A c c e p t e d M a n u s c r i p t calibration module. Once the depth-maps are computed, a smoothing algorithm is applied (Gastal and Oliveira, 2011), this is mostly to compensate for the staircase effect caused by discretized depth values. The different steps of the procedure are presented in Figure 3. Unlike the Morpho3D device, this method does not require immobilization of the animal. However, four stainless steel cables are also used within the Deffilait3D Scanner, on both side of cow, to secure the scanner and restrict cow movement, but cows cannot be restrained by a self-locking head gate. If necessary, a bucket of pellets is sometimes used, especially during the first few passes, to reassure the animal (but also to reward it and encourage it to return). The speed of this device allows high-speed acquisition of 3D images, but can also be used to record almost continuously as the delay between photographic triggers can be reduced to a fraction of a second.

Data

In the IEPL facility, the Morpho3D and Deffilait3D devices are located in the same room. 3D images of a given animal were then collected under the same conditions, a few seconds or minutes apart. In both cases, the images were analyzed using Metrux8α® software (3D Ouest, Lannion, France). When using Morpho3D device, the 3D acquisitions were not always fully exploitable due to the animal's movements or the necessity to block the animal thanks to the self-locking head gate for example. As a result, images were sometimes cut off at the shoulder level and we developed then predicting model of total surface area (TS) or volume (TV) from partial surface area (PS) or volume (PV) (Le Cozler et al., 2019b). When the image of full animal is available, TS and TV are calculated using the automatic algorithm. When cutting off is necessary, PS and PV are measured by placing one plane at the tip of the shoulder blades and another at the rump (Le Cozler et al., 2019b).

A c c e p t e d M a n u s c r i p t

Metrux8α® software has functions to quickly measure interactively:

the partial or total volume (PV or TV, respectively) of an animal; its partial or total surface area (PS or TS, respectively); height at the withers (WH), hip width (HW), buttocks width (WB), diagonal length (DL, obtained by averaging the left and right diagonal lengths, i.e., between the tips of the buttocks and the tips of the shoulders), chest depth (CD), and heart girth (HG) (see Le Cozler et al., 2019a for more details on locations). To these initial measurements that have already been presented in previous publications, we also added the abdomen circumference (AC). These measurements are not automatic and the points of interest must be "manually clicked".

Thereafter, these measures are performed thanks to geometric computations (distances, curve lengths etc.) based on features points selected by the user. For the determination of the volume, the method presented by Mirtich (1996) is used.

The images acquired using the Deffilait3D and Morpho3D devices were of variable quality. Each image was therefore manually assigned a quality score varying from 0 to 4, with the following characteristics:

-0: unusable point clouds, i.e. when there was too much movement by the animal; -1: reconstructed image usable for some measurements but still of poor quality (e.g., totally deformed legs); -2, 3, and 4 correspond respectively to average, good, and very good images, on which most or all measurements were feasible.

We also added a "+" sign after each score to indicate whether the image of the animal was complete or not, e.g., with or without the head. Systematically, because of the differences between the A c c e p t e d M a n u s c r i p t gantries, images from the Deffilait3D device always had a "+" while those from the Morpho3D device did not. Examples of scoring are shown in Figure 4.

Validation of the Deffilait3D device

In order to test the validity and utility of the Deffilait3D device, we carried out a series of trials in which the repeatability and reproducibility of the device were assessed. Analysis were performed on 15 cows, in two trials -Repeatability is "the closeness of agreement between successive results obtained with the same method on the same material under test and under the same conditions" (ISO Standard 3534). In this case, "under the same conditions" meant that the measurements were made with the same operator, the same equipment, in the same place, and in a very short time interval. In the context of points, an analysis of repeatability makes it possible to know if, for a given acquisition, the identification of the points for predicting, e.g., the same partial volume, is correct. For this purpose, eight cows were chosen, these cows being considered as representative of the variability in BW and size observed in this herd. The morphology and the chosen indicators (HW, WH, WB, DL, PS, PV, CD, HG, AC) were determined from a single image acquisition, from which the measurements were repeated five times by the same operator. In this way, we tested the repeatability of the point identification.

-Reproducibility is "the closeness of agreement between individual results obtained with the same method on the same test material under different conditions" (ISO Standard 3534). Different conditions may be related to different operators or different locations, for example, and intermediate cases are possible, such as different operators working in the same location or measurements made with the same equipment but with a long interval between acquisitions. It is thus possible to define and calculate reproducibility between A c c e p t e d M a n u s c r i p t operators, between devices, or across days. Here, we tested intra-facility reproducibility, which evaluates how well the method is able to predict the same data in different acquisitions in a changing environment (i.e. with the animals in slightly different positions). The assessment takes into account the image capture by the scanner, the fusion of the 15 images, from the 15 pairs of sensors and laser used, image cleaning, reconstruction, and measurement. For the reproducibility test, eight cows were used, of which seven were different from those used in the repeatability tests. As in previous tests, they were considered to be representative of the variability observed in BW and size in this experimental herd. Each was subjected to five different sessions of image acquisition, taken at relatively short time intervals (10 to 20 min). For each image, the same measurements were captured as for the repeatability test, one time by the same operator.

Comparison of the Deffilait3D and Morpho3D devices

Once the repeatability and reproducibility of the Deffilait3D device were validated, its measurements were compared with those obtained from the Morpho3D device, which were considered the reference values ("gold standard"). Fifteen adult dairy cows were randomly selected and passed ("scanned") through each device. These cows differed from those used in previous repeatability and reproducibility tests. Following the acquisition and reconstruction of the 3D images from each device, selected morphological measurements were determined and compared.

For the sake of time, we limited our assessment to this panel of 15 cows and values of HW, WH, CD, HG, TV, and TS. In order to limit animal stress as much as possible and to avoid retaking measurements unnecessarily, the animals were physically restrained in a feed fence in the Morpho3D device. For the surface and volume measurements, TS and TV were then calculated x TV + 408 x HW + 271 x WB -199-was then used to predict BW. Other equations for predicting BW are available but are not presented here (this one was selected because it gave the most accurate results).

Statistical analysis

Data were first visualized using the ggplot2 package (v. 3.2.1) (Wickham et al., 2019) of R software (R Core Team, 2019). Other statistical analysis were realised under R version 3.5.2 with "car" package for analysis of variance and regression, to compare values from Morpho3D and Deffilat3D devices (Fox and Weisberg, 2011). In general, measurement values may vary among the different 3D images used (five images in total for each cow) based on slight differences among images as well as differences in point identification within a single image. By definition, then, the reproducibility error encompasses the repeatability error of point identification and will therefore always be the higher value of the two. The parameters used to quantify repeatability and reproducibility were the standard deviations of point repeatability, σrp, and intra-facility A c c e p t e d M a n u s c r i p t reproducibility, σRi, which were calculated from the residuals of the one-factor analysis of variance (ANOVA) model (the cow effect). The ANOVA, as well as other repeatability and reproducibility calculations, were performed with R software (R Core Team, 2019). With this approach, existing morphological variations between cows are removed in the analysis. The coefficients of variation of point repeatability, CVrp, and intra-facility reproducibility, CVRi, are then expressed as follows:

𝐶𝑉𝑟𝑝= 𝜎𝑟𝑥̅ 𝑟×100 𝐶𝑉𝑅𝑖= 𝜎𝑅𝑥̅ 𝑅×100
Where 𝑥̅ 𝑟 and 𝑥̅ 𝑅 are the respective means of the study variable in the populations used for the repeatability and reproducibility assessment.

The lower the coefficient of variation, the more repeatable or reproducible the method or part of the method is. In other words, r and R represent in each case the arithmetic difference between two determinations with a 95% probability of not being exceeded. Repeatability and reproducibility are considered good if the CVs are below 4%, which was used as the reference value for this study (Fischer et al., 2015;Le Cozler et al., 2019a).

RESULTS

Repeatability and reproducibility tests to validate the "Deffilait3D" device. The repeatability and

reproducibility tests used a total of 40 (8 x 5) and 200 (8 x 5 x 5) images, respectively, with all being of sufficient quality (score of 3+ minimum). For all measurements, the coefficients of variation of point-related repeatability (CVrp) in the images were low (below 1.8%), with three CVrp values below 0.5% (HW, WB, AC; AC; Table 1). Coefficients of variation of reproducibility varied from 1.12 to 3.26%, with an average CVRi of 2.07%. The least reproducible indicators were partial surface area (3.26%) and complete volume (2.71%), and these were probably due to imperfect images, including poorly reconstructed limbs or heads. Concretely, for a cow with a volume of 764 litres (L), the repeatability and reproducibility errors were 8.71 and 21.41 L, respectively. For the same cow with a total surface area of 7.20 m², the errors in repeatability and reproducibility of surface area were 0.13 m² and 0.23 m², respectively. Finally, for an average withers height of 1.45 m, the repeatability error was 0.019 m and the reproducibility error was 0.026 m. The full repeatability and reproducibility tables for each cow used in the present experiment are available in Supplementary File 1.

Morphological measurements confirm that the two devices are equivalent.

Selected traits of cows used in trials 1, 2, and 3 are presented in Table 2. All these average values were determined from 3D images recorded from Morpho3D device, as it was considered as the reference method in these trials. Surface areas and volumes were 7.20 m 2 and 764 L in trial 1, 7.18 m 2 and 736 L in trial 2 and 6.28 m² and 581 L in trial 3, respectively. The height at the withers in trials 1, 2 and 3 was 1.45, 1.50, and 1.38 m, respectively (Table 2). A notable difference (about 0.08 m) in height at the withers (WH) was observed between the two devices (Table 3). For the other measurements, the differences between Morpho3D and Deffilait3D were of lesser magnitude (0.013 m for HW, 0.004 m for WB, 0.025 m for CD, 0.007 m for AC), and represented less than a 4% difference. For total volume, only 1.2 L of difference was detected between the two devices, which is negligible compared to the volume of a cow (about 750 L). For the total surface area, the values obtained by the two devices differed by 0.6 m², i.e. about 8%.

To better illustrate the correlations between the measurements obtained using Morpho3D and Deffilait3D, a total of 200 measurements of HW, WB, CD, and AC were taken from images collected on the same cows with the two devices, and are compared in Figure 4. HW measurements were very similar between the two devices, with the R² and the coefficient of the equation (0.97) being close to one. For measurements of WB, CD, and AC, the R² values were lower (respectively, 0.70, 0.55, and 0.66), but with the coefficient of the equation again close to 1. In the case of WH, the R² coefficient was very good (0.99), but the coefficient of the equation was lower (Figure 5). This could be due to differences in the positions of animals during image acquisition (see discussion section). For TS and TV, for which comparisons were based on PS and PV, the values were largely similar between the two devices, but with some differences in TS. This was probably due to the smoothing procedure used for images from the Deffilait3D device, which requires further optimization (Figure 6).

Application of the device to determine animal weight.

When we compared measured values of body weight to those that were predicted by the equation of Le Cozler et al. (2019a) based on data from the Deffilait3D device, we found that the R² coefficient was excellent (0.99), as was the coefficient of the equation (Figure 7).

DISCUSSION

The device, Deffilait3D, can be used to accurately measure or estimate, in a single shot, most morphological traits commonly used on farms. These data can then be combined for use in predicting body 7 Deffilait3D was designed and developed for the high-throughput acquisition of both partial and whole 3D images, which allows work to be carried out on the entire animal or only a part of it (e.g., body condition score, at the back of the animal). Because our earlier device, Morpho3D, had been extensively validated in previous work (Le Cozler et al., 2019a), we used the measurements obtained from it as the reference values ("gold standard") for the assessment of the Deffilait3D system. It should be noted, however, that the notion of a "gold standard" or "reference measurements" can be open to discussion: the fact that a value is used as a reference does not necessarily mean that it is more accurate or authoritative in all cases. Newly developed approaches and methods can sometimes give values closer to the true value than the reference method. This is especially the case for measurements for which the true value is notoriously difficult to determine, such as volume or surface area. Since 3D imaging technologies have become available only recently, the methods used for morphological measurements have traditionally relied on the use of a tape (HG) and/or height gauge (WH, CD, etc.). In this trial, no manual measurements were performed because the Morpho3D measurement tool has been tested and validated for several years, and the resulting data are now considered as references.

For the majority of indicators, morphological measurements generated using the Deffilait3D device were repeatable and reproducible. For linear measurements (e.g., withers height, hip width), partial volume, and partial surface area, the repeatability coefficients obtained here ranged from 0.4 to 1.8% and were similar to those observed for the Morpho3D device (Le Cozler et al., 2019a). Even the highest CVs for repeatability and reproducibility, which were associated with complete volumes and partial areas, were still below the theoretical threshold of 4%, and should thus still be considered very good (Fischer et al., 2015). The slightly lower performance of the Deffilait3D device compared to its predecessor can probably be explained, at least in part, by the lower resolution of the one-shot acquisition strategy, which resulted in lower quality 3D image reconstructions. This resulted in more artifacts (defects, small blisters, etc.) on images. However, during the experiment, we were able to make some improvements that resulted in improved image quality. Further calibration and modifications will make it possible to obtain smoother and betterreconstructed 3D images and will undoubtedly reduce measurement variability in the future.

Numerous studies have been carried out to estimate or perform body measurements in cattle [START_REF] Auriol | Production laitière et conformation dans la race montbéliarde. I. -héritabilité de la production laitière et de différentes mesures corporelles. Corrélations phénotypiques et génétiques entre ces mesures et la production laitière[END_REF][START_REF] Heinrichs | Predicting body weight and wither height in Holstein heifers using body measurements[END_REF]Tozser and Bedo, 2000), but few measurements have been performed using a 3D imaging device. In general, the main morphological measurements (HW, WB, TV, and TS) were very similar between the two devices compared here, with occasional minor differences (e.g., +0.4 cm between Deffilait3D and Morpho3D for WB). Point identification errors (on 3D images, in software such as Metrux6α®) seemed equally likely regardless of the device used. Although the results obtained with the two systems were largely comparable, the Deffilait3D has one considerable advantage: it is able to record the overall shape of an animal in only a few tenths of a second, compared with several minutes for Morpho3D. This new device therefore seems poised to make a significant contribution to the goal of high-throughput and accurate phenotyping. In the future, it may be possible to use successive images of the same animal, separated by a few tenths of a second, to generate dynamic data similar to what is possible with a strobe effect, but this remains to be further investigated.

A c c e p t e d M a n u s c r i p t

The largest difference was found for measurements of withers height, for which the two devices differed by about 8 cm. There are several possible explanations, including a distorted reconstruction of the underside of the legs with Morpho3D or a difference in the positioning of animals between the two devices. To address this, a "standard position" probably needs to be determined in the future, with which it may be possible to reduce the variation between devices.

For example, Ling et al. ( 2022) highlighted the importance of a standard position in correctly determining body size measurements in pigs and Li et al. (2023) proposed a posture based measurement adjustment method to improve accuracy measurements based on point cloud data in beef cattle. In the Morpho3D device, the heads of animals are restrained thanks to the use of a self locking head gate and therefore always raised, while with the Deffilait3D device, the animal's head may move freely, which can affect legs position and thus, withers height. Here, the average withers height measured with the two devices was 154.0 cm and 145.9 cm, respectively. In an ongoing project using Deffilait3D device (Tiercin et al., unpublished data), wither height was corrected according to head and legs positions for a group of cows and the difference disappeared.

Considering animal posture while registering 3D images is then of importance, at least for wither height measurement.

Using previously published equations validated in both adult and growing animals (Le Cozler et al., 2019b;2022a), it was possible to determine the body weight of cows. Here, we present the results of only one of the published equations, but results were similar with the others (not shown).

In the future, the improved access to traits and/or indicators granted by the use of 3D imaging will likely make it possible to develop new equations that may be used even in the absence of such devices. In the meantime, the speed and ease of use of devices such as Deffilait3D may enable continuous monitoring of body weight of a large number of animals throughout growth and/or A c c e p t e d M a n u s c r i p t adult life. As a step toward this goal, we are developing automatic extraction techniques based on a robust geometric algorithm as well as machine learning to improve automation of certain steps that must still be performed manually (e.g., preparing and "cleaning" the figures, determination of some measurements). Considerable progress has been made in the use of computer vision and deep learning methods, which are increasingly widespread in the field of agricultural data analysis (Kamilaris and Prenafeta-Boldu, 2018;[START_REF] Espejo-Garcia | End-to-end sequence labelling via deep learning for automatic extraction of agricultural regulations[END_REF] and livestock production [START_REF] Bezen | Computer vision system for measuring individual cow feed intake using RGB-D camera and deep learning algorithms[END_REF]. In the particular case of using convolutional neural networks (CNNs), a dataset that includes both the image and its associated target variable is used to train the neural network. If the training is successful, the network can be used to create measurements from new data. CNNs can be used on depth maps but are unsuited for unordered point clouds (although networks for point clouds are emerging). Another approach is to provide the learning algorithm a normalized representation of the 3D data, as described by Fischer et al. (2015).

Finally, for an on-farm use, the number of cameras will decrease, as it has been done for the Pheno3D device, which is a mobile device being tested since the beginning of 2023, on around 2000 growing beefs, from16 farms (Lebreton et al., 2023). Pheno3D is based on Deffilait3D development, using only 10 cameras, and costs around 18,000 US $. Preliminary results confirmed that this movable3D scanner allows high-throughput phenotyping of the whole cattle body for animal in motion, but also development of models to predict body weight and the Linear scored used by the French stakeholders to do selection on morphological criteria (Lebreton et al., 2023).

A c c e p t e d M a n u s c r i p t CONCLUSION Based on our results, the Deffilait3D device appears to be as accurate and precise as its predecessor, Morpho3D, in generating morphological data from cows, but can acquire images considerably faster. Measurements obtained with Deffilait3D were highly repeatable and reproducible for all of the studied traits, with coefficients of variation of less than 1.8% for repeatability and less than 3.3% for reproducibility. In the future, the development of automated reconstruction and analysis tools based on machine learning methods will facilitate improved extraction and/or analysis of indicators of interest at a very high rate, which will greatly increase the utility of this technology and its potential for widespread adoption in livestock farming or even on large wild animals. Such an approach opens the door to high-throughput acquisition of traits based on the whole animal as well as other partial measurements or analyses, such as rumen development, animal mobility, lameness detection, or body condition score.

  Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024 A c c e p t e d M a n u s c r i p t

  Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024 A c c e p t e d M a n u s c r i p t from PS and PV, which correspond to the surface area and volume without the head and neck, using previously published prediction equations (Le Cozler et al., 2019a): TV (total volume) = 1.0704 × PV (volume cut off at the shoulder) + 0.015 TS (total surface area) = 1.07 × PS (surface area cut off at shoulder tip) + 0.94 Validation of published equations In order to determine whether the Deffilait3D device could provide reliable and useful values from prediction equations established with the Morpho3D device, we tested the equation predicting body weight (BW) on a group of 13 heifers, which represented a high degree of variation in BW measured on a weighing scale (440 to 550 kg). The morphological traits described above were determined for these animals, and the equation presented by Le Cozler et al. (2019a)-BW = 644

  Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024

  From a mathematical point of view, repeatability (r) and reproducibility (R) are each defined by the value below which the absolute value of the difference between two individual results (x1 and x2) obtained under either repeatability or reproducibility conditions lies with a specified probability. If not otherwise specified, the probability is 95%.i.e.: P (|x1-x2| ≤r) = 0.95 and P (|x1-x2| ≤R) = 0.95

  Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024 A c c e p t e d M a n u s c r i p t

  Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024 A c c e p t e d M a n u s c r i p t

  Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024 A c c e p t e d M a n u s c r i p t

  Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024 A c c e p t e d M a n u s c r i p t

Figure 1 .

 1 Figure 1. The Morpho3D device at the Experimental Milk Production Facility (IEPL)

Figure 6 .

 6 Figure 6. Comparison of measurements of withers height (WH) between Morpho3D and

Figure 7 .

 7 Figure 7. Comparison of measurements of total volume (TV) and total surface area (TS) between

Figure 8 .

 8 Figure 8. Comparison of body weight (BW) measured using a weighing scale with values

Figure 1

 1 Figure 1

Figure 2

 2 Figure 2

Figure 3

 3 Figure 3

Figure 4

 4 Figure 4

Figure 5

 5 Figure 5

Figure 6

 6 Figure 6

Figure 7

 7 Figure 7

Figure 8

 8 Figure 8

Table 1 .

 1 Errors and coefficients of variation of repeatability and reproducibility for morphologicalDownloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024

	measurements using the Deffilait3D device
		Coefficient of variation (%)	Error
		Repeatability Reproducibility Repeatability	Reproducibility
	Complete volume, L	1.18	2.71	8.71	21.14
	Partial surface, m²	1.77	3.26	0.13	0.23
	Hip width, m	0.48	1.12	0.0028	0.0064
	Withers height, m	1.29	1.78	0.0186	0.0264
	Chest depth, m	0.81	1.42	0.0064	0.0116
	Heart girth, m	0.46	2.16	0.0118	0.047
	Buttocks width, m	1.34	1.83	0.007	0.01
	Diagonal length, m	0.76	2.55	0.013	0.0447
	Abdominal	0.46	1.8	0.0118	0.0475
	circumference, m

Table 2 .

 2 Average values of selected traits measured on 3D images of animals recorded with used Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024

	Morpho3D
	Trial number TS (m²) TV (L) WH (m) HW (m) WB (m)	CD (m)	AC (m)
	1 (n=16)	7.202	763.9	1.454	0.575	0.537	0.811	2.632
	2 (n=15)	7.181	739.2	1.499	0.565	0.530	0.787	2.587
	3 (n=13)	6.278	581.2	1.376	0.507	0.507	0.755	2.324

TV: total volume; TS: total surface; WH: withers height; HW: hip width; WB: buttocks width; CD: chest depth; AC: abdominal circumference.

Table 3 .

 3 Comparison of the Morpho3D and Deffilait3D systems based on measurements obtained : withers height; HW: hip width; WB: buttocks width; CD: chest depth; HG: heart girth; AC: abdominal circumference; TV: total volume; TS: total surface area. Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024

	from 15 adult cows
		WH (m) HW (m) WB (m)	CD (m)	AC (m)	TV (L)	TS (m²)
	Morpho3D	1.540	0.578	0.526	0.813	2.580	739.8	7.481
	Deffilait3D	1.459	0.565	0.530	0.787	2.587	738.6	6.881
	Difference	0.0811	0.013.2	0.0037	0.0252	0.0067	1.2	0.6
	(in %)	(0.0053)	(0.0023)	(0.0070)	(0.0031)	(0.0026)	(0.0016)	(0.0080)

WH

Downloaded from https://academic.oup.com/tas/advance-article/doi/10.1093/tas/txae018/7604297 by INRA -DOCUMENTATION user on 13 February 2024

Available at

Deffilait3D project, especially the technicians at the Mejusseaume experimental station, who took excellent care of the animals. FUNDING INFORMATION: The authors are grateful to the national Deffilait (Agence Nationale de la Recherche, Paris, France, ANR-15-CE20-0014), and APIS-GENE (Paris, France) for funding this research. The authors also thank #DigitAg for the financial support during the internship of Maxime Dumesny.