Agniva Sengupta

Alexandre Krupa

Eric Marchand

Visual Tracking of Deforming Objects Using Physics-based Models

In this paper, we propose a framework for tracking the deformation of soft objects using a RGB-D camera by utilizing the physically-based model of the considered object. A coarse, 3D template of the object being tracked is the only prior information required by the proposed method. The proposed approach does not rely on the accurate knowledge of the material properties of the object being tracked. In this paper, we integrate computer vision based tracking methodology with physical model based deformation representation without requiring expensive numerical optimization for minimizing nonlinear error terms. The proposed approach enables deformation tracking by joint minimization of a geometric error and a direct photometric intensity error while utilizing co-rotational Finite Element Method (FEM) as the underlying deformation model. The proposed method has been validated both on synthetic data (with groundtruth) and real data.

I. INTRODUCTION

Physically based models have been used to analyze deformations of objects of complex shapes since the inception of structural mechanics. Some of the recent computer-vision based approaches have tried to leverage the advantages of the physics based models for deformation tracking. Surprisingly, these physical-model based approaches often fall short of other competing deformation models in terms of accuracy [START_REF] Hoppe | A benchmark and evaluation of non-rigid structure from motion[END_REF]. Many deformation tracking methodologies used by the computer vision and robotics community utilizes geometric regularizers or assumption of isometry or conformality as a deformation model. However, structural mechanics offer many highly evolved methodologies for modelling deformation in complex shaped objects. In this paper, we propose a method of non-rigid object tracking using a RGB-D camera that combines the minimization of a geometric error that penalizes the model to pointcloud distance and a photometric error that penalizes inter-frame photometric intensity difference. A coarse, 3D template of the object being tracked is the only prior input necessary for the proposed approach. The entire observable surface of the deforming object is tracked in 3D using the data from the RGB-D camera. The deformation model used for representing the underlying physics is based on FEM. Thanks to the physics based model, it is possible to track large deformations that do not preserve the volume of the object. The proposed method is robust to occlusion and the tracking accuracy is not dependent on apriori knowledge of physical properties of the deforming object. However, if the physical parameters are known apriori, we demonstrate that the approach can also be utilized for estimating the external contact forces applied on the surface of the deforming object. The following are the main contributions of this paper:

• We propose a new methodology for deformation tracking using joint optimization of a geometric and photometric error defined on the surface model of the deforming object. A non-linear least squares optimization based method with an analytic expression of the Jacobian relating the variation of geometric and visual error to displacement of the vertices of the mechanical mesh is proposed to track the deformation of objects • The method proposed here has also been utilized to track external forces applied on the surface of deforming objects. Such applications are an important advantage of using physically-based models for deformation tracking.

II. BACKGROUND

A method to track deformable objects by applying virtual forces on the simulation of a physically based model was proposed in [START_REF] Schulman | Tracking deformable objects with point clouds[END_REF]. The physics of the deforming object was modelled as a collection of linked rigid bodies or particles and the tracking was formulated as an Expectation-Maximization (EM) problem. The qualitative results in [START_REF] Schulman | Tracking deformable objects with point clouds[END_REF] are impressive, but there are a few sequences where the tracking suffers noticeable drift in heavily occluded regions. It is not trivial to modify the approach of [START_REF] Schulman | Tracking deformable objects with point clouds[END_REF] to utilize a Newtonian optimization scheme for minimizing the visual error. [START_REF] Wang | Deformation capture and modeling of soft objects[END_REF] improved the approach of [START_REF] Schulman | Tracking deformable objects with point clouds[END_REF] by considering co-rotational FEM. But, in [START_REF] Wang | Deformation capture and modeling of soft objects[END_REF], no attempt was made to tackle the deformation tracking problem using conventional optimization since the gradient estimation of the high dimensional optimization step was considered to be intractable. In this context, it must be noted that physical model free, RGB-D based deformation tracking approaches such as [START_REF] Richard A Newcombe | Dynamicfusion: Reconstruction and tracking of non-rigid scenes in real-time[END_REF] depend heavily on the assumption that the deformation is isometric. Consequently, these approaches are not preferable when there is a significant change in the volume of the deforming object. In [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF], the authors utilized co-rotational FEM to track deformation by purely geometric matching of pointcloud with the 3D object model, producing highly accurate tracking results at frame-rate. However, the approach of [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF] was also sensitive to occlusion, large interframe motion and demonstrated a tendency to drift along the tangent of large planar surfaces, which can be intuitively explained as an inherent drawback of tracking solely using point-to-point or point-to-plane correspondences. In [START_REF] Sengupta | Tracking of non-rigid objects using rgb-d camera[END_REF], the deformation tracking using physically based model has been formulated as a conventional optimization problem, but only using depth information. Moreover, the optimization in [START_REF] Sengupta | Tracking of non-rigid objects using rgb-d camera[END_REF] requires expensive simulation of object deformation using FEM for computing the gradient of the visual error, which is unsuitable for deformation tracking at frame-rate. [START_REF] Sengupta | Simultaneous tracking and elasticity parameter estimation of deformable objects[END_REF] proposes a force tracking methodology based on the approach of [START_REF] Sengupta | Tracking of non-rigid objects using rgb-d camera[END_REF] which is slow and does not utilize the photometric information for tracking.

III. METHODOLOGY

The proposed method tracks deformable objects using a RGB-D camera with the help of a coarse 3D template of the deforming object. The proposed approach can be summarized by the following iterative steps:

• Rigidly track the deforming object w.r.t the camera • Analyze the residual error after rigid tracking to determine the region which is likely to have deformed • Minimize a combination of geometric and photometric error terms to track the deformation using a set of selected vertices of the model in those deformed region We now describe the entire approach in details.

A. Notation

Throughout the rest of the paper, the notation Q is used to denote the vertices of the volumetric mechanical tetrahedral mesh K, while P S = (P S X , P S Y , P S Z ) denotes the vertices of the surface model S. The mapping S → K is done via barycentric mapping. Following standard convention, the object to camera homogeneous transformation matrix is given as

C T O = C R O C t O 0 1
and this is utilized to align K and S to the actual object, where C R O and C t O are the rotation matrix and translation vector respectively. Between two consecutive data frames, the pose of the object in current frame (p) with respect to pose in the previous frame (p -1) is q = ( p t p-1 , θu), where θ and u are the angle and axis of rotation p R p-1 . The time derivative of q is given as p v p-1 = δq, where v ∈ se(3) is the velocity screw. The color image for the p-th data frame is given by I p . The perspective projection function Π(•) is used to project the 3D point P to its corresponding 2D coordinates (u, v) in the image plane, such that

  u v 1   = Π P 1 = K p Π p c T o P 1
, where K p is the camera calibration matrix in the standard notation and Π p = I 3×3 0 3×1 . The function I(u, v) gives the grayscale intensity at that pixel position and (∇I u p , ∇I v p ) gives the corresponding image gradient. (f x , f y ) are focal lengths in pixel units.

B. Deformation Model

We employ FEM [START_REF] Nealen | Physically based deformable models in computer graphics[END_REF] to model the underlying physics of the non-rigid object. Let the four vertices of a particular face of the tetrahedral mechanical mesh be denoted by the vertices Q 1 , Q 2 , Q 3 and Q 4 and the centroid be a point Q C . To obtain the displacement of vertices as a result of application of external force on the object, a second order differential equation of the following nature needs to be solved:

M Q + D Q + KQ = F ext + F e (1)
where M and D represent the mass and damping matrices of the model respectively, and K is the global stiffness matrix.

F ext denotes the external forces acting on the vertices and F e are the internal elastic force vector acting on the vertices of the tetrahedral elements. A linear solver based on conjugate gradient descent [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF] is used to solve Eqn. [START_REF] Hoppe | A benchmark and evaluation of non-rigid structure from motion[END_REF]. Co-rotational FEM [START_REF] Muller | Interactive virtual materials[END_REF] is a modification of linear FEM that allows robust handling of strong deformations in the mesh. If the model undergoes a large deformation, let the centroid of the tetrahedral element at this deformed configuration be denoted by a point Q R C . An arbitrary vertex Q x of the mechanical mesh in the undeformed frame gets deformed to Q in the new frame (variables with the • notation are expressed in the corotated frame Q R C ). With the co-rotational formulation [START_REF] Felippa | A systematic approach to the element-independent corotational dynamics of finite elements[END_REF], it can be shown that

d = R U -Q C + Q x -Q x given that R is the rotation matrix, U is the displacement of Q x , expressed in the object frame of Q C , and d is the displacement of Q x in the corotated frame Q R C .
Applying similar deformation to the vertices of the tetrahedron, the internal, elastic forces can be represented by:

F e = R e K e ŨR (2)
Here, ŨR = ( Q1 , Q2 , Q3 , Q4 ), K e represents the stiffness matrix and R e is the 12 × 12 block diagonal matrix of four 3 × 3 rotation matrices R stacked diagonally.

C. Visual Tracking

We track the deforming object by minimizing two error terms defined on the surface of S. The minimization is done using a set of control handles C on K, such that the entire deformation becomes optimizable by regulating the displacement of C. C are those vertices of K that shall be controlled to minimize the error terms defined on S. Any displacement of C gets propagated to all vertices of K using the deformation model. We are interested in minimizing a geometric error function E depth and a photometric error function E photo . The combined error term can be represented as:

E S = E depth E photo (3)
E S is minimized using a set of displacement vector on C. For the sake of simplicity, we describe the proposed mechanism of minimizing E S with the help of a single control handle, i.e., the u-th control handle of C, denoted by u C . However, to optimize in a non-linear least squares fashion, we need to obtain the gradient of the stacked cost function. Numeric estimation can be computationally expensive as determining the node displacements typically involves multiple iteration of conjugate gradient descent for solving an equation of the nature of Eqn. [START_REF] Hoppe | A benchmark and evaluation of non-rigid structure from motion[END_REF]. One of the key contribution of this paper is the proposed mechanism for minimizing E S with non-linear least squares optimization along with a strategy for analytically estimating the Jacobian.

Let us analyze the case of an arbitrary 3D point P of the pointcloud, which lies on or near the triangle

P S i P S i+1 P S i+2 of S.
Let us assume that a small displacement ∆u C = ∆u CX ∆u CY ∆u CZ produces the displacement ∆x, ∆y and ∆z on the j-th surface plane of S (due to the underlying deformation of the mechanical mesh) comprising of three vertices P S i , P S i+1 and P S i+2 respectively. This vector of vertex positions is given by:

ϑ j = (P S i ) (P S i+1 ) (P S i+2 ) (4)
Subsequently, the gradient of the error E S (P) w.r.t u C is given by the Jacobian:

J S (P) = ∂E S (P) ∂u C = ∂E S (P) ∂ϑ j 1×9 ∂ϑ j ∂u C 9×3 (5)
The term ∂ϑj ∂uC is somewhat similar to the classical straindisplacement matrix, as expressed conventionally in FEM literature [START_REF] Friswell | Finite element model updating in structural dynamics[END_REF]. The details of how ∂E S (P) ∂ϑj can be derived in our case is described in Sec. III-E.

Computing ∂ϑj ∂uC can be done by estimating the deformation of the mechanical model K, for every new control handle explored by our proposed approach, only once. For the j-th control handle C j , this is done by successively displacing the vertex in K corresponding to C j by a small distance along the positive direction of X, Y and Z axis. This estimation step produce three deformed meshes per control handle, which can be directly utilized for determining ∂ϑj ∂uC using forward finite differences. Assuming that u C points to the c-th vertex of K, the values of

∂P S i ∂uC can be mapped to a matrix Γ O by a mapping function C such that Γ O = C ∂P S i ∂uC where: ∂ϑ j ∂u C = ∂P S i ∂uC , ∂P S i+1 ∂uC , ∂P S i+2 ∂uC (6)
Here, C represents the mapping

∂P S i ∂uC → Γ O .
Assuming that M denotes the number of vertices in K and S, C can be expressed as:

Γ O = C ∂P S i ∂u C = M i=0 M c=0 ∂P S i ∂u C ⊗ J i,c M×M (7)
where ⊗ gives the Kronecker product and J i,c M×M denotes the conventional single-valued matrix of dimension [M × M] for the c-th control handle, such that:

J i,c M×M x,y = 1 if x = i and y = c, ∀x, y ∈ M 0 else ( 8
)
and the x and y in the subscript of (•) x,y denotes the row and column index of an element of the matrix. Based on the matrix indices, Γ O can be represented as:

Γ O aM+i,bM+c = ∂P S i ∂u C a,b | ∀i, c ∈ M ∧ a, b ∈ {0, 1, 2} (9)
We term the matrix Γ O as influence matrix (since it denotes the influence of the displacement of a vertex on its neighbors) and it can be conveniently computed offline once per every object model and stored. During tracking, Γ O is loaded from memory when required.

For estimating Eqn. ( 5), we derive a new matrix R Γ such that:

R Γ 3M×3M = C R O ⊗ I M×M ( 10
)
where

I M×M is an identity matrix of dimension [M × M].
R Γ allows us to rotate the influence matrix to the current camera reference frame. This transformed matrix Γ C is given by Γ C = R Γ Γ O . We then derive a new matrix Γ by thresholding Γ C such that its value at an arbitrary index m, n associated with the i-th vertex P S i is given by:

Γ m,n = 1 if ε > 3α||∆u C || 2 0 else (11)
given that ε =

2 a=0 2 b=0 Γ C m+a,n+b 2
where ∆u C is the calibrating deformation (implying that ||∆u C || is a constant for a particular Γ O ). α is a tunable parameter which regulates the area for which the value of Γ will be 1.

This thresholding of Eqn. ( 11) is done because we intend to utilize the point-to-plane distance as the error term E depth . It is known classically that point-to-plane distance minimization using planar (or nearly planar) surface model can result in the model 'sliding' over the 3D data [START_REF] Chen | Object modeling by registration of multiple range images[END_REF]. The calibrating deformation ∆u C produces larger deformation in vertices closer to the vertex index c, i.e., Γ C highly prioritizes vertices close to c. This, in turn, converts the Jacobian in Eqn.

(5) to effectively consider only a small patch of surface near the vicinity of the c-th vertex, while neglecting the small magnitude of deformation observed in the vertices farther away from the c-th vertex. To overcome this problem, we propose to binarize Γ C to Γ for obtaining the gradient of E depth . However, the Jacobian for E photo can be obtained directly from Γ C without modification.

D. Determining the Control Handles

It is possible to estimate suitable positions for the control handles by analyzing the cost function. Given the two objective functions E depth and E photo , we define two vectors Z depth and Z photo , such that their values at the index corresponding to the 2D image point Π(P l ) is given by:

Z depth Π(P l ) = (W depth ) l,l E depth (P l ) Z photo Π(P l ) = (W photo ) l,l E photo (P l ) (12)
for the l-th point P l of the pointcloud. Z depth and Z photo forms a vectorized image from the pixel intensities corresponding to the error values of E depth and E photo respectively. W depth and W photo are the diagonal weighting matrix derived from E depth and E photo using Tukey based m-estimator [START_REF] Da Freedman | On inconsistent m-estimators[END_REF]. The combined error matrix Z Π on the image plane is obtained by:

Z Π = vec -1 H×W Z depth vec -1 H×W Z photo (13)
where denotes the Hadamard product and • gives the normalized matrix. vec -1 H×W (•) denotes the vector to matrix map R HW → R H×W , where H and W are the height and width image respectively. A median blur with 3 × 3 sized kernel is applied on Z Π followed by a linear thresholding and the output of this operation is clustered into multiple clusters. The centroid of each cluster is associated with the nearest projection of the visible vertices of K on the image plane. These associated vertex indices of K are identified as the control handles for that particular frame.

E. Non-rigid Error Minimization

Assuming that an estimate of p T p-1 has already been obtained (using the method described in Sec. III-F), a combination of depth based geometric error and direct photometric error is minimized to track the deforming object. We now define the error term that needs to be minimized for a single point P of the pointcloud at the p-th data frame. The pointto-plane distance based geometric error is given by:

E depth (P p ) = n j • P p -d j (14)
assuming that the j-th surface plane of S (where n j is the normal and d j is distance to origin) corresponds with the point P on the image plane. We propose a photometric error term defined on the (p-1)-th frame by:

E photo (P p-1 ) = I p (P p-1 e ) -I p-1 (P p-1 ) (15)
The updated point position P p-1 e is determined by a barycentric map P p-1 e = (P S i ) (P S i+1 ) (P S i+2 ) B such that P p-1 corresponded with the triangle P S i P S i+1 P S i+2 and (P S i ) gives the updated vertex position of P S i when subject to the update ϑ, given that B is a column vector denoting the barycentric coordinates of P p-1 w.r.t P S i , P S i+1 and P S i+2 . Eqn. 3 can be slightly modified to:

E S = E depth µE photo ( 16
)
as the combined cost function, where µ = β E depth E photo is used for bringing the geometric and photometric error terms to the same scale and 0.9 ≤ β ≤ 1.2 is a tunable parameter for weighting the relative influence of the photometric cost function.

The combined cost function that we seek to minimize is Eqn. ( 16) w.r.t u C , the displacement of the control handle. The Jacobian relating the change of u C to the change of E S is obtained by utilizing Eqn. [START_REF] Muller | Interactive virtual materials[END_REF] and Eqn. [START_REF] Felippa | A systematic approach to the element-independent corotational dynamics of finite elements[END_REF] in Eqn. [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF], such that:

J = ∂E S ∂u C = ∂E depth ∂ξ Γ µ ∂E photo ∂ξ Γ C ( 17
)
where:

∂E depth ∂ξ =    -n j -n l A P S i+2 -P S i+1 × ) -n l (A P S i -P S i+2 × -n l A P S i+1 -P S i ×    ( 18
)
given that:

A = 1 n j (I 3×3 -n j n j ) (19)
provided (•) × gives the skew-symmetric matrix and:

n j = (P S i+2 -P S i ) × (P S i+1 -P S i ) (20)
and n l = P p l -P S i . On the other hand:

∂E photo ∂ξ = ∇I u p ∇I v p        b1 0 0 b1 b2 0 0 b2 b3 0 0 b3          G i P S G i+1 P S G i+2 P S   (21)
where: The optimization method of our choice is Levenberg-Marquadt (LM) like [START_REF] Moré | The levenberg-marquardt algorithm: implementation and theory[END_REF], and the update is given by:

G q P S =   fx (P S q )Z 0 -f x (P S q )X (P S q ) 2 Z 0 fy (P S q )Z -f y (Pq)Y (P S q ) 2 Z   ∀q ∈ i, (i+1), (i+2) (
δu C = -(H + λI N×N ) -1 J W WE S ( 23
)
where H = J W WJ is the approximation of the Hessian and W = diag W depth , W µ photo where W depth and W µ photo is obtained using the Tukey based m-estimator on E depth and µE photo respectively. λ is a scaling factor and N is the size of the error vector E S . The new position u C + δu C of C j is to be applied on the mechanical model such that it deforms to minimize the errors from ( 14) and [START_REF] Moré | The levenberg-marquardt algorithm: implementation and theory[END_REF]. Apart from α in Eqn. [START_REF] Felippa | A systematic approach to the element-independent corotational dynamics of finite elements[END_REF], β related to Eqn. ( 16) and λ (and excluding the clustering mechanism of Sec. III-D), there are no parameters involved in the non-rigid visual tracking approach presented here.

F. Initialization

The rigid pose of the object in the camera's reference frame is updated at the beginning of each frame using an approach similar to [START_REF] Trinh | A modular framework for model-based visual tracking using edge, texture and depth features[END_REF], by jointly minimizing a geometric error and a sparse feature based error, given by:

e D ( p v p-1 ) = p R p-1 P p-1 + p t p-1 • n j -d j (24) e K ( p v p-1 ) = Π(P p-1 ) x -x * Π(P p-1 ) y -y * (25)
where P p-1 represents an arbitrary 3D point from the last data frame, which has been matched to the j-th plane of the object's 3D model, denoted by the normal vector n j = (n X j , n Y j , n Z j ) and distance to origin d j . Π(P p-1 ) x , Π(P p-1 ) y are the projection of the same arbitrary point in the (p-1)-th image while (x * , y * ) are the same image points matched in the p-th image using Harris corner features. p v p-1 ∈ se(3) gives the velocity twist between the (p -1)-th frame and p-th frame.

G. Force Tracking

Next, the method of Sec. III-C -III-E is used to estimate the deforming forces acting on an object. Two additional information are required for force tracking: the material properties of the object being tracked, i.e., Young's modulus, Poisson's ratio, Rayleigh stiffness etc. and the approximate point of contacts on the object. The force tracking experiments have been performed by tracking the tool applying the force using fiducial markers. Once a point P P have been identified as the point of contact, its nearest mechanical mesh vertex Q P is obtained using a nearest neighbor search.

Thereafter, the force applied on P P is given as i F Qi , for all indices i such that Q i is a neighbor of Q P , wherein the force vectors F Qi is derived using (2).

IV. RESULTS

The results can be roughly divided into three categories:

• For the simulated objects with groundtruth, we focus on comparing the fundamental concept of the non-rigid tracking approach proposed here with comparable stateof-the-art methods. We also use challenging simulated sequences to quantitatively validate the robustness of our approach; • For deformation tracking on real data, we emphasize the capability of our algorithm to accurately track large volumetric deformations in soft objects; • We provide experimental validation of the force tracking approach proposed in Sec. III-G

A. Simulated Results: Validation with Groundtruth

We base our simulated results on two objects 1 , a cuboid and an icosphere, as shown in Fig. 3 and Fig. 2. The simulated data is generated using the Blender software [START_REF]Blender -a 3D modelling and rendering package[END_REF] and the object deformations are generated manually using harmonic coordinates with B-spline basis [START_REF] Joshi | Harmonic coordinates for character articulation[END_REF]. Some quantitative results on synthetic sequences are expressed as a percentage of the largest diagonal of the bounding box of the objects (2.884m for the cuboid and 3.408m for the icosphere), and hence the values are unitless.

First, we demonstrate that the difference in accuracy between maintaining the influence matrix Γ constant (C. Γ) and re-computing and updating Γ (U. Γ) at each frame numerically is not significant. To establish this proposition, we run tests on the synthetic data with and without holding Γ constant for both standard, linear FEM (SF) and corotational FEM (CRF). The results, summarized in Table I, are expressed in terms of the Hausdorff distance (H) from groundtruth (GT). It can be clearly seen that the variation in accuracy of tracking is between 0 to 0.59 % for both the sequences tested here. This comes in exchange for a large improvement in runtime (upto > 53% improvement in per frame time requirement, see Sec. IV-D for time requirements of the implementation) of the entire algorithm, since multiple FEM simulations per iteration of LM is highly expensive. We also use these sequences for comparing the proposed approach to [START_REF] Petit | Tracking elastic deformable objects with an rgb-d sensor for a pizza chef robot[END_REF] and [START_REF] Sengupta | Tracking of non-rigid objects using rgb-d camera[END_REF]. The corotational FEM based version of the proposed approach outperforms [START_REF] Petit | Tracking elastic deformable objects with an rgb-d sensor for a pizza chef robot[END_REF] and [START_REF] Sengupta | Tracking of non-rigid objects using rgb-d camera[END_REF] in both the sequences. Fig. 2 shows the visual comparison of all the approaches tested on the synthetic sequences. Using a specific control handle and maintaining all other parameters constant, the variation of the tracking accuracy (in terms of H) of the method proposed in this paper w.r.t the Young's modulus (Y) and Poisson's ratio (σ) is shown in Fig. 1. It is clear that the tracking accuracy do not vary significantly with change in Young's modulus. This observation is in accordance with [START_REF] Sengupta | Tracking of non-rigid objects using rgb-d camera[END_REF]. However, the tracking accuracy varies significantly with change in Poisson's ratio. This is because materials with Poisson's ratio in the range of 0.35 -0.45 are highly ductile in nature and hence significantly harder to track with. Fig. 3 shows the output of the proposed approach on two challenging deformation sequences on cuboid and icosphere. The two objects were subjected to large volumetric deformation and were occluded using large floating objects in the synthetic scene. The mean value of H is 0.0696m for the cuboid and 0.0836m for the icosphere, while the value of H for every frame of the sequence if shown in Fig. 4.

SF CRF Input [19] [6] U. Γ C. Γ U. Γ C. Γ

B. Experiments on Real Data

The real data has been captured using an Intel RealSense D435 camera. Three objects were used, a block of sponge, a rugby ball (cut in half) and a soft dice. These objects were strongly deformed from the top. The output of tracking has been demonstrated visually2 in Fig. 5 while the evolution of point-to-plane distance across the two sequences have been logged in Fig. 6. The sponge and ball demonstrates highly accurate tracking with a mean geometric error (from Eqn. 14) of -0.81 mm and -0.18 mm respectively, with a standard deviation of 0.47 mm and 0.62 mm. This accuracy is highly robust to occlusion and demonstrates the suitability of the proposed approach to robotic manipulation of soft objects. The dice is a slightly smaller deformation, but shows a high mean accuracy of ∼ 0.61 mm.

C. Force Tracking

To validate the force tracking methodology, we use a 6-DOF anthropomorphic robot arm (Viper 850 from ADEPT) fitted with a ATI Gamma IP65 force/torque sensor and a 3Dprinted stylus as an end-effector distal tool. The robot is used only to utilize its force sensor to obtain a GT of the force. We use the sponge and the ball to validate the force tracking method. The Young's modulus of the sponge and the ball is determined by repeated indentation tests and is found to be 460 kPa and 160 kPa respectively. Next, these objects are subjected to a strong deformation using the robot's 3D printed end-effector, and the point of contact on the object is tracked using pre-trained markers (fitted to the probe) from the image data. The results of the force tracking approach, as summarized in Fig. 7, shows an accuracy of ∼ 97% for the sponge and of ∼ 90% for the ball when the deformation is optimally observable.

D. Implementation

The approach proposed in this paper has been implemented in C++ on an Intel Xeon CPU working at 3.70 GHz. Utilizing only a single core of the computer (without using GPU), the un-optimised code was able to achieve a runtime of 350 ms -550 ms per frame while tracking deformations using the proposed approach, showing that it can be possible to achieve real-time performance at frame rate.

No. of data frame The paper presented here describes a new method for combining geometric and photometric visual error minimization with FEM to create an accurate and fast deformation tracking method. The algorithm has been tested on synthetic and real data and has been shown to outperform state-of-the-art methods in tracking accuracy for generic deforming objects. The algorithm proposed here can be extended to other, more complex physical models without loss of generalization. The method proposed here has been demonstrated as a reliable visual force tracking system, when the material properties of the object being tracked is available.

  22) and B = b 1 , b 2 , b 3 are the barycentric coordinates.

Fig. 1 :

 1 Fig. 1: Variance of H with Y and σ for the two simulated objects, cuboid and icosphere (color coded with the value of H)

Fig. 2 :

 2 Fig. 2: Comparison between GT and tracking output for the two synthetic sequences cube and icosphere. The red edges and vertices shows the object model from the GT, the black ones are from the tracking output.

Fig. 3 :

 3 Fig. 3: The cuboid and icosphere are subjected to large deformation ((a) & (b)) with large, floating objects flying across the scene (as occlusion). The tracking output is shown in (c) and (e) while E depth is shown on the surface of S in (d) and (f)

Fig. 4 :

 4 Fig.4: The value of H (in m) across all the frames of the two sequences shown in Fig.3

Fig. 5 :Fig. 6 :

 56 Fig. 5: The color and depth image input for the ball, sponge and dice sequences are shown in (a), (b), (f), (g), (k) & (l). The tracked model are shown in (c), (h) & (m). (d), (i) & (n) shows the object model overlaid on the image, while (e), (j) & (o) shows the geometric error for the visible surfaces, derived using[20]

Fig. 7 :

 7 Fig.7: Setup and results for force tracking. The red plot is the force measurement from the robot's force sensor and the blue plot is the estimated force reported by the proposed approach

1

  All data available at: github.com/icra2021/VisualDeformationTracking

			Cube	Icosphere
	SF	Constant Γ Updated Γ	5.76 % 5.75 %	4.83 % 5.42 %
	CRF	Constant Γ Updated Γ	2.14 % 2.14 %	5.25 % 5.34 %
		[19]	2.93 %	7.93 %
		[6]	3.28 %	6.67 %

TABLE I :

 I Tracking accuracy for two synthetic sequences in terms of Hausdorff distance between tracking output and GT

Univ Rennes, Inria, IRISA, CNRS, Rennes, France e-mail: {agniva.sengupta, alexandre.krupa, eric.marchand}@irisa.fr.

For a detailed video of the results, please visit: youtu.be/ScJnz j4-cs

ACKNOWLEDGEMENT

We gratefully acknowledge the support from The Research Council of Norway through participation in GentleMAN (299757) project.


Volodymyr M Hiiuk

Karl Ridier

A Il'ya

' Gural

Alexander A Skiy

Igor O Golub

Gábor Fritsky

William Molnár

Azzedine Nicolazzi

Bousseksou

Gural'skiy

Alexander A Golub

Igor O Fritsky

Gábor Molnár

William Nicolazzi

Azzedine Bousseksou
email: azzedine.bousseksou@lcc-toulouse.fr

Influence of the ultra-slow nucleation and growth dynamics on the room-temperature hysteresis of spin-crossover single crystals

Keywords: Spin crossover, Nucleation and growth dynamics of the phase, Hysteresis, Defects, Single crystals

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

Introduction

Iron (II) spin-crossover (SCO) complexes constitute one of the most interesting representatives of switchable molecular materials [START_REF] Gütlich | Spin Crossover in Transition Metal Compounds I-III[END_REF][START_REF] Halcrow | Spin-Crossover Materials: Properties and Applications[END_REF][START_REF] Bousseksou | Spin crossover phenomenon/Phénomène de transition de spin[END_REF]. The SCO phenomenon in these coordination complexes is related to the ability of the central metal ion to change its spin state, in a reversible manner, under the influence of external stimuli [START_REF] Gütlich | Thermal and Optical Switching of Iron(II) Complexes[END_REF]. These complexes, which have the particularity of exhibiting bistability properties, characterized by the possible emergence of hysteretic phenomena at technologically relevant temperatures, constitute appealing molecular systems for practical applications, for example, as active elements of information processing, energy transducing and display devices [START_REF] Kahn | Spin-Transition Polymers: From Molecular Materials Toward Memory Devices[END_REF].

Due to the drastic change of different physical (magnetic, optical, mechanical, electrical, etc.) properties during the SCO, this phenomenon can be monitored by various techniques (magnetometry, X-ray diffraction, Mössbauer spectroscopy, differential scanning calorimetry, optical microscopy, etc.) [START_REF] Gütlich | Thermal and Optical Switching of Iron(II) Complexes[END_REF]. Among these techniques, optical microscopy (OM) occupies a special place because it allows a real-time, non-destructive and non-contact detection of the SCO in individual single crystals, combining high temporal and spatial resolutions. (N.B. Here we primarily refer to OM in bright-field, transmission mode, which probes local absorbance changes and allows for quantitative spatiotemporal studies of the spin-state transformation. It should be noted, however, that OM has been also used to probe refractive index and luminescence intensity changes, enabling to detect the spin transition in single SCO micro-nano-particles [START_REF] Arnaud | Observation of an asymmetry in the thermal hysteresis loop at the scale of a single spin-crossover particle[END_REF][START_REF] Quintero | Soft lithographic patterning of spin crossover complexes. Part 1: fluorescent detection of the spin transition in single nano-objects[END_REF][START_REF] Liu | Imaging the Thermal Hysteresis of Single Spin-Crossover Nanoparticles[END_REF]).

OM studies on large single crystals (typically 50-100 m) help to better understand the fundamental aspects of the spin-state switching process (nucleation and growth phenomenon, spatiotemporal properties), which have an influence on the key characteristics of the spin transition (switching temperatures, sharpness, completeness and hysteretic behaviour).

The first (visible-range) OM studies on SCO single crystals were performed by Jeftic and Hauser [START_REF] Jeftic | Patterns During Photoexcitation and High-Spin → Low-Spin Relaxation in [Fe(ptz)6](BF4)2 Spin Transition Crystal[END_REF]. These authors observed the spatiotemporal development of the high-spin (HS) to low-spin (LS) relaxation in a single crystal of the SCO complex [Fe(ptz)6](BF4)2 (ptz = 1-propyltetrazole). The observed inhomogeneous features were explained due to structural modifications and internal stresses in the crystal associated with the spin-state transition. In the past decade, numerous OM studies on different SCO single crystals allowed significant advances in the understanding of the spatiotemporal properties of the spin transition with the identification of some general features, which have been reviewed by different groups involved in this research [START_REF] Varret | Real-Time Observation of Spin-Transitions by Optical Microscopy[END_REF][START_REF] Ridier | Hysteresis, nucleation and growth phenomena in spin-crossover solids[END_REF]. In summary, in highly cooperative compounds exhibiting hysteretic behaviour, heterogeneous nucleation followed by the formation of a well-defined HS/LS interface and its slow propagation (typically between 1-10 m s -1 , but sometimes up to a few hundreds of m s -1 ) was observed in both cooling and heating modes. Moreover, it was shown that the nucleation and growth phenomenon significantly depends on the crystal morphology, the local microstructure and the presence of internal defects. Several studies have indeed pointed out the prominent role of crystal edges and defects on the nucleation process and the phase boundary propagation [START_REF] Fourati | Interplay between a crystal's shape and spatiotemporal dynamics in a spin transition material[END_REF][START_REF] Traiche | Spatio-temporal aspects of the domain propagation in a spin-crossover lattice with defect[END_REF],

while the possibility of modulating the spatiotemporal dynamics of the nucleation and growth phenomenon by the creation of artificial micrometric defects has been experimentally demonstrated [START_REF] Bedoui | Laser-Induced Artificial Defects (LIADs): Towards the Control of the Spatiotemporal Dynamics in Spin Transition Materials[END_REF]. Another general observation is that the SCO is often accompanied by irreversible damage (cracks, fractures, dislocations, etc.) upon subsequent switching cycles, which mainly originate from the internal stresses generated by the transition process -due to the relatively large volume change (and lattice mismatch) between the HS and LS phases. Overall, the substantial volume change in SCO solids appears as a key parameter at the origin of the slow phase-boundary propagation velocity and the rather low mechanical resilience of SCO crystals [START_REF] Ridier | Hysteresis, nucleation and growth phenomena in spin-crossover solids[END_REF].

In the present work, we report on OM investigations of SCO single crystals of [Fe(1,6naphthyridine)2(Ag(CN)2)2] (1). This two-dimensional (2D) coordination network exhibits a reproducible, isostructural spin transition characterized by a wide thermal hysteresis loop centred near room temperature (297 K) [START_REF] Hiiuk | Room temperature hysteretic spin crossover in a new cyanoheterometallic framework[END_REF]. Here we show that this hysteretic region is associated with singular kinetic effects and exceptionally slow switching rates, whose origin can be traced back to the highly defective nature of the crystals.

Experimental

Single crystals of [Fe(1,6-naphthyridine)2(Ag(CN)2)2] (1) were obtained by a slow diffusion method as previously described by some of us [START_REF] Hiiuk | Room temperature hysteretic spin crossover in a new cyanoheterometallic framework[END_REF] (see also the Supporting Information, SI for the synthesis details). OM studies have been performed on several single crystals of 1. The data shown in the manuscript refers to a crystal with dimensions of 44 µm × 39 µm × 10 μm. The crystal was placed on a glass slide and then enclosed in a temperature-controlled heating-cooling stage (THMS350V, Linkam Scientific Instruments). Before the experiments, the stage chamber was purged with dry nitrogen gas for ca. 3 min at 373 K. Photographs were recorded in brightfield transmission mode using an Olympus BX51 upright microscope equipped with a ×50 long working-distance objective (numerical aperture, NA = 0.5) and a CCD camera (Clara, Andor Technology). Image acquisition was conducted either isothermally or at fixed heating/cooling rates (from 0.1 to 10 K min -1 ). A band-pass filter centred at λ = 532 nm was used to obtain a maximum optical contrast between the LS and HS phases. Image processing was performed using the ImageJ software.

Results and discussion

Recurrent pitfalls for spatiotemporal studies on single crystals is the sample-to-sample variation of the SCO properties due to different sample morphologies and microstructures and, for a given crystal, the progressive formation of defects and associated drift of the measured properties upon thermal cycling. These issues often make difficult the identification of the transition mechanism as well as the quantification of the associated kinetics. As a first step, we have thus carried out a series of preliminary measurements on different crystals of 1. The most important findings are summarized in Figure 1. Overall, it appears that minor details of the thermal SCO curves differ both from crystal-to-crystal and from cycle-to-cycle, but the SCO behaviour appears sufficiently robust and reproducible to conduct more extensive investigations on single crystals of 1.

The hysteretic behaviour stems from the existence of a macroscopic free energy barrier (arising from cooperative intermolecular interactions), which defines a region of thermodynamic meta-stability between the LS and HS phases. Close enough to the equilibrium temperature, the macroscopic energy barrier is sufficiently high to make the lifetime of metastable states virtually infinite (compared to the experimental time scale), thus defining a quasi-static hysteretic region.

However, as the lifetime of the metastable states progressively decreases as moving far from the equilibrium temperature, a kinetic dependence of the thermal hysteresis can be observed, characterized by an artificial increase of the hysteresis width. Indeed, numerous studies clearly demonstrated the strong dependence of the shape, sharpness, width and/or position of the hysteresis loop upon the temperature scan rate [START_REF] Brooker | Spin crossover with thermal hysteresis: practicalities and lessons learnt[END_REF]. A spectacular example was reported on the mixed crystal series of the two-dimensional network {[Zn1-xFex(bbtr)3](ClO4)2}∞ (bbtr = 1,4di(1,2,3-triazol-1-yl)-butane). In this compound, temperature scan-rate studies have revealed that the apparent hysteresis loop emerging at low temperature (ca. 60-70 K) is of purely kinetic origin, owing to the slow HS-to-LS relaxation rate at this temperature [17]. Another remarkable example is that of the complex [Fe II (nBu-im)3(tren)](PF6)2 ((nBu-im)3(tren) = tris(n-butylimidazol(2-ethylamino))amine), which is characterised by an unusual increase of the thermal hysteresis width upon reducing the temperature scan rate (from 14-K-wide at 4 K min -1 to 41-Kwide at 0.1 K min -1 ) [START_REF] Seredyuk | Unprecedented Multi-Stable Spin Crossover Molecular Material with Two Thermal Memory Channels[END_REF]. Such an effect is explained by the fact that two distinct LS phases, with different structural organizations, can be kinetically formed from the same HS phase, depending on the temperature sweep rate. Moreover, an important experimental fact is that the width of the thermal hysteresis may also simply depend on the thermal exchange coupling with the sample environment. As shown for single crystals of the SCO compound [Fe(HB(1,2,4-triazol-1yl)3)2], the transition temperatures, as well as the width of the thermal hysteresis loop, are independent of the temperature scan rate (up to at least 5 K min -1 ) at atmospheric pressure. However, by decreasing the air pressure down to P = 0.1 mbar, noticeable changes were observed both on the thermal hysteresis width and on the switching kinetics of the crystal, due to the reduced efficiency of convective heat exchange between the sample and its environment [START_REF] Ridier | Scan-rate and vacuum pressure dependence of the nucleation and growth dynamics in a spin-crossover single crystal: the role of latent heat[END_REF]. The problem of heat exchange on the SCO hysteresis properties was also identified and investigated in single crystals of the SCO compound [{Fe(NCSe)(py)2}2(m-bpypz)] (py = pyridine and bpypz = 3,5-bis(2-pyridyl)pyrazolate) [START_REF] Traiche | Spatiotemporal Observation and Modeling of Remarkable Temperature Scan Rate Effects on the Thermal Hysteresis in a Spin-Crossover Single Crystal[END_REF].

Figure 2 depicts the typical spin-transition curves of a single crystal of 1, obtained by varying the temperature scan rate from 0.1 to 10 K min -1 . The values of the spin-transition temperatures in both heating and cooling modes (T1/2 ↑ and T1/2 ↓ , respectively) are summarized in Table S1 for the different scan rates used. As shown in Figure 2a, the LS-to-HS transition (heating mode) remains abrupt (ΔT < 10 K) whatever the scan rate. In contrast, the transition curves in the cooling mode exhibit a more noticeable scan-rate dependence. Indeed, while the HS-to-LS transition is rather abrupt at low scan rates (ΔT ~ 10 K at 0.1-1 K min -1 ), much more gradual transitions are observed for high rates (ΔT ~ 30 K at 5-10 K min -1 ). Thus, the increase of the temperature scan rate induces a pronounced enhancement of the asymmetry of the thermal hysteresis loop. This observation may suggest the existence of different driving/frictional forces for the spin transition in the two branches of the thermal hysteresis loop. As an example, the HS-to-LS and LS-to-HS transitions might involve distinct frictional forces due to the difference of stiffness between the LS and HS states. On the other hand, as depicted in Figure 2b, the shift of the transition temperatures is similar in the two branches: T1/2 ↓ (respectively T1/2 ↑ ) varies from 290.5 K (respectively 298 K) to 280 K (respectively 308 K) for scan rates of 0.1 K min -1 and 10 K min -1 , respectively. As a result, the width of the thermal hysteresis goes from 7.5 K up to 27.5 K when increasing the scan rate from 0.1 to 10 K min -1 .

Obviously, the symmetric widening of the hysteresis loop might be indicative of thermalization problems at high scan rates. However, our previous studies [START_REF] Ridier | Scan-rate and vacuum pressure dependence of the nucleation and growth dynamics in a spin-crossover single crystal: the role of latent heat[END_REF] using the same experimental setup discard this "extrinsic" origin of the hysteresis broadening in the present study. To investigate the plausible origin of such kinetic effects in 1, isothermal relaxation experiments were performed. The sample was cooled/heated at 1 K min -1 to various given temperatures within the hysteresis loop and the optical transmission across the crystal was monitored as a function of time at constant temperature (Figure 3). Further details on these isothermal experiments, including the relaxation curves, are provided in the SI. It is worth noting that a similar approach has been already used for other SCO systems to distinguish the quasi-static and "apparent" hysteresis loop [START_REF] Stassen | Physical properties of the spin-crossover compound hexakis(1-methyltetrazole-N 4 )iron(II) triflate, steady state and relaxation studies, X-ray structure of the isomorphic Ni(II) compound[END_REF]. It is important to stress, however, that in most cases the "apparent" hysteresis broadening has been observed at low temperatures (typically below ca. 100-120 K), mostly on the cooling branch, due to the slow spin-state interconversion rates at those temperatures.

In the case of 1, as depicted in Figure 3a for the cooling mode, no relaxation is observed (hour time scale) down to a temperature of 291.8 K, this temperature thus defining the descending branch of the quasi-static hysteresis loop. For lower temperatures (for example, 291.0 K) a slow relaxation is observed to the full LS state. One should note that the long relaxation times encountered here (up to 1 hour) definitely exclude the extrinsic origin (sic thermalization) of the observed phenomena. The same procedure in the heating mode allows to unveil a quasi-static 8-K-wide hysteresis loop (Figure 3b). The width of this quasi-static hysteresis turns out to be approximately the same as that observed at the slowest scan rate (ΔT = 7.5 K at 0.1 K min -1 ). Such a quasi-static hysteresis implies the existence of a "real" thermodynamic bistability region (characterized by long-lived metastable states) near room temperature in this SCO material, which is an interesting asset for practical applications.

The spatiotemporal observation of the thermally induced spin transition in 1 also provides crucial information for the rationalization of the slow transition kinetics associated with the hysteresis. Figure 4 displays a set of selected images and movies showing the nucleation event(s) and the real-time evolution of the phase transition in the cooling and heating modes (see also the SI). These visual demonstrations reveal that the presence of numerous defects in the crystal undoubtedly impacts the spin-transition process. Indeed, the crystals of 1 appear systematically imperfect and can be better represented as a system consisting of many, small structural domains, which are more or less randomly oriented. This high degree of mosaicity was observed both via optical microscopy and single-crystal X-ray diffraction. † These different observations illustrate the so-called multi-scale causal chain of the SCO process proposed recently [START_REF] Guionneau | Multiscale Approach of Spin Crossover Materials: A Concept Mixing Russian Dolls and Domino Effects[END_REF].

Indeed, it appears that microscopic details (defects, mosaicity) can have a huge impact on the macroscopic behaviour of the material, in particular here concerning the spin-state switching dynamics.

In the cooling regime, inhomogeneous nucleation of the LS phase simultaneously occurs in different places of the crystal, mostly near the edges, and several, non-consecutive phase boundaries propagate across the entire sample. On the contrary, in the heating mode, the nucleation of the HS phase always takes place in the middle of the crystal and a better-defined phase boundary propagates toward the edges of the crystal in a roughly circular manner. These features are found well reproducible for different crystals and for several consecutive thermal cycles. However, it should be noted that after typically more than 10 cycles some irreversible damage were observed in certain crystals (defects, cracks, etc.).

The singular spatiotemporal dynamics of the spin-state transformation in crystals of 1 can be qualitatively reproduced by kinetic Monte Carlo simulations, using a dynamic version of a spinphonon model. This class of models is known to be relevant for the theoretical investigations of the spatiotemporal development of the spin transition [START_REF] Enachescu | Elastic models, lattice dynamics and finite size effects in molecular spin crossover systems[END_REF], since the lattice degree of freedom and the structural differences between the HS and LS phases are explicitly considered (see details in the SI). By introducing structural inhomogeneities within the system (i.e., a local tensile stress in the middle of the crystal lattice), we show that we are able to qualitatively reproduce some key features of the nucleation and growth phenomenon experimentally observed in crystals of 1. In particular, distinct nucleation sites for the LS (edges) and HS (middle of the crystal) phases as well as the existence of a quasi-circular phase-boundary propagation across the crystal can be depicted in the simulated images (Figure S6). This however raises questions about the possible origin of such a residual elastic stress gradient within the crystals of 1, especially in view of the synthesis method (slow diffusion), which might possibly lead to the build-up of stress inhomogeneities within the crystals.

A distinctive feature of the spin transition in crystals of 1 is the exceptionally slow kinetics of the growth process of the nucleated phase. Table 1 displays typical switching times (t), temperature ranges (T) and phase-boundary propagation velocities recorded in a crystal of 1.

Remarkably, the switching time at low scan rates (0.1 K min -1 ) is ca. 30-50 min, which is ca. 4 orders of magnitude longer than the switching times (50-100 ms) previously reported for crystals of the SCO complex [Fe(HB(1,2,4-triazol-1-yl)3)2] (2) [START_REF] Ridier | Scan-rate and vacuum pressure dependence of the nucleation and growth dynamics in a spin-crossover single crystal: the role of latent heat[END_REF] -in closely comparable experimental conditions. The difference of switching times clearly parallels the difference of phase-boundary velocities typically observed in the two materials. In crystals of 1, we observe propagation velocities in the range of 10 -2 -10 -1 m s -1 , down to 3 nm s -1 (Table 1), whereas these velocity values fall in the 10 2 -10 3 m s -1 range in crystals of 2 [START_REF] Ridier | Scan-rate and vacuum pressure dependence of the nucleation and growth dynamics in a spin-crossover single crystal: the role of latent heat[END_REF]. Previous studies of the nucleation and growth phenomenon in various SCO systems indicate that the average propagation velocity of the LS/HS phase boundary is typically 1-20 µm s -1 , which is in between the two extreme cases represented by compounds 1 and 2 (Table S2) [START_REF] Varret | Real-Time Observation of Spin-Transitions by Optical Microscopy[END_REF][START_REF] Ridier | Hysteresis, nucleation and growth phenomena in spin-crossover solids[END_REF]. The exceptionally slow switching kinetics in crystals of 1 must certainly come from the existence of particularly strong resistive forces in this compound, whose effect is to hinder the

Conclusions

In conclusion, we have performed an optical microscopy study of the thermally induced spin transition in single crystals of the SCO compound [Fe(1,6-naphthyridine)2(Ag(CN)2)2], including temperature scan-rate measurements, the assessment of the quasi-static hysteresis loop as well as the spatiotemporal investigation of the spin-state switching phenomenon. Despite the high transition temperatures (around room temperature), the thermal hysteresis loop is marked by pronounced kinetic effects, which appear to be related to the exceptionally slow growth rate of the nucleated phases. This slow switching process, which is a well-reproducible, general feature of this compound, is characterized by extremely slow propagation velocities of the LS/HS phase boundary (down to 3 nm s -1 ), the slowest known to date. Our OM studies reveal the critical influence of the unusually high number of microstructural defects in the crystals, which constitute a myriad of energy barriers to overcome, making the transition properties strongly dependent on the temperature scan rate. An important perspective of this work will be the investigation of size reduction effects on the switching kinetics. Indeed, as the number of structural domains decreases one can expect a significant acceleration of the switching rateapproaching the intrinsic limit, represented by a perfect single crystal.

Notes

  phase-boundary motion. The main contribution of resistive forces commonly arises from the emergence of internal frictions and mechanical stresses during the phase-boundary propagation, due to the structural mismatch between the LS and HS phases. In addition, the role of crystal defects appears more dramatically in the present study. Cracks, failures and structural domain boundaries are numerous in crystals of 1 and obviously represent additional energy barriers for the propagation of the LS/HS phase boundary, making the domain growth mechanism a thermally-activated process at the origin of the "anomalous" hysteretic behaviour. The role of crystal defects is exemplified in Figure 4 (blue arrow), where a large failure is shown to open and close in a crystal during the spin transition.

  †For the single crystal described in Chem. Commun. 2019, 55 (23), 3359-3362 [15] the mosaicity parameters are 2.54 (e1 = 2.35, e2 = 1.10, e3 = 4.16) and 2.46 (e1 = 1.82, e2 = 1.09, e3 = 4.48) degrees in the LS and HS forms, respectively (CrysAlisPro).

Figure 1 .

 1 Figure 1. (a) Reproducibility of the SCO in four different single crystals of 1. (b) Reproducibility of the SCO in a single crystal of 1 for three subsequent thermal cycles. All data were acquired at 1 K min -1 . Heating and cooling are indicated by arrows.

Figure 2 .

 2 Figure 2. (a) Temperature scan-rate study (0.1-10 K min -1 ) of the spin transition (normalized optical transmission vs. temperature) in a single crystal of 1. (b) Scan-rate dependence of the spin-transition temperatures in the heating and cooling modes.

Figure 3 .

 3 Figure 3. (a) Isothermal HS → LS relaxation measurements on a crystal of 1 at selected temperatures. A reference curve acquired at constant cooling rate (1 K min -1 ) is included (violet line). (b) Quasi-static hysteresis loop (green region) extracted from the isothermal relaxation experiments. Two reference curves, recorded in heating and cooling modes at scan rates of 1 K min -1 (black line) and 0.1 K min -1 (grey line), are also shown.

Figure 4 .

 4 Figure 4. Snapshots of the spin-state transformation of a crystal of 1 in the cooling and heating modes at 0.1 K min -1 . Red arrows indicate the propagation of phase boundaries. The blue arrow shows the opening and closing of a failure in the course of the phase transition.

Table 1 .

 1 Switching temperature range (ΔT), switching time (Δt) and typical phase-boundary propagation velocities (v) measured for a single crystal of 1 at various scan rates (dT/dt) in the heating and cooling modes.

		dT/dt	ΔT	Δt	v		dT/dt	ΔT	Δt	v
		(K min -1 )	(K)	(min)	(nm s -1 )		(K min -1 )	(K)	(min)	(nm s -1 )
		0.1	3	30	13		0.1	5	50	3
	Cooling mode	0.5	9	18	23	Heating mode	0.5	11	22	14
		1.0	13	13	32		1.0	11	11	24
		5.0	35	7	160		5.0	10	2	140
		10	35	3.5	160		10	10	1	240

Acknowledgements

This work was supported by H2020-MSCA-RISE-2016 Project 734322 and Projects 19BF037-01M and 19BF037-04 of the Ministry of Education and Science of Ukraine. V. Hiiuk thanks Lozynskyj Foundation and Yuchymenko Family Endowment Fund for the financial support. Olesia I. Kucheriv is acknowledged for useful comments on the manuscript.

Conflicts of interest

The authors confirm that this article content has no conflict of interest.