Xiao Xie

Johann Herault

Étienne Clement

Vincent Lebastard

Frédéric Boyer

Quasi-static motion of a new serial snake-like robot on a water surface: a geometrical approach

This paper reports methods to compute the equilibrium stances of a new snake-like robot designed to stabilize its head on a free water surface. To adjust rapidly the stability of the robot, this bio-inspired robot can rotate independently each body-shell, and modify the level of immersion of each module. To predict the stable stance accessible by this additional degree of freedom, a model is developed to compute the equilibrium configurations of the robot from a given parametrization of the body shape. Then, an algorithm is introduced to compute a sequence of controlled body deformations, such that the head configuration relatively to the water surface remains unchanged. Finally, we explore in simulation stances and quasistatic gaits, and investigate to what extent the buoyancy and the body deformations can be used to stabilize the head of the snake-like robot.

I. INTRODUCTION

Marine robots inspired by eels and aquatic snakes has emerged in academic contexts [START_REF] Hirose | Snake-like robots[END_REF], [START_REF] Ijspeert | From swimming to walking with a salamander robot driven by a spinal cord model[END_REF], [START_REF] Liljebäck | Eelume: A flexible and subsea resident IMR vehicle[END_REF]. Still under development for industrial underwater applications [START_REF] Liljebäck | Eelume: A flexible and subsea resident IMR vehicle[END_REF], these hyper-redundant serial robots are particularly compact and manoeuvrable, and offer an alternative to standard unmanned water surface vehicles. Despite the high potential of these bio-inspired robots, these applications on water surface remains complex due to their precarious stability. Indeed, their slender morphology, which is the main asset of their announced performance, makes them particularly sensitive to rolling motions (solid rotation about the longitudinal axis). The net torque of the forces of buoyancy associated with certain stances can destabilize the robot, and produces high acceleration due to their small longitudinal moment of inertia. Hence, this rolling motion is considered to be the most dangerous degree of freedom. The problem of instability of planar swimming is particularly obvious when robots [START_REF] Porez | Improved Lighthill fish swimming model for bio-inspired robots: Modeling, computational aspects and experimental comparisons[END_REF], as well as real snakes, swim at low velocity. This parasitic rolling motion is a serious impediment of the visual guidance of these bio-inspired robots, as well as the strong effects on the swimming performances [START_REF] Porez | Improved Lighthill fish swimming model for bio-inspired robots: Modeling, computational aspects and experimental comparisons[END_REF]. However, stabilization of snake-like robot swimming on water surface is poorly documented, and most control-oriented models consider neutrally buoyant robot moving in a 2D plane, except in [START_REF] Kelasidi | Modeling of underwater snake robots moving in a vertical plane in 3D[END_REF]. To prevent the rolling motion, the robot can be equipped with floaters or attached to cables, but these passive systems reduces the manoeuvrability of the robots.

To overcome these issues, the present article proposes an additional degree of freedom that stabilizes the rolling body motion while preserving the performance of the planar anguilliform gait. This new device allows to modify rapidly the immersed volume by a local controlled rolling motion on each module having triangular-like cross sections (see Fig. 1). This stabilization strategy is inspired from semiaquatic snakes living on river surface, such as Cottonmouth [START_REF] Herault | Standing on the Water: Stability Mechanisms of Snakes on Free Surface[END_REF]. These snakes deform actively their body shape, thanks to complex combinations of flexions and torsion, to monitor the surrounding environment. This active stabilisation relies on a control of the distribution of buoyancy forces along their slender body. Taking inspiration from these aquatic snakes [START_REF] Herault | Standing on the Water: Stability Mechanisms of Snakes on Free Surface[END_REF], we aim to design a robot and its control to achieve stabilization on water surface via an active control of buoyancy forces. The long-term goal is to maintain a gaze control on water surface paving the way to visual servoingbased navigation of autonomous snake-like robots. However, this additional degrees of freedom makes the design and the control of the robot more complex. In practice, it raises the central issue: how can these internal degrees of freedom be combined to stabilize and control the head configuration while moving on a flat water surface? This paper is a first attempt to address this problem by investigating the effect of the body deformations on the stability. As a first step, we analyse the motion of the robot as a succession of static stable stances. This fictive quasistatic motion consists in neglecting all the inertial effects (fluid+body) as well as the viscous force produced by the fluid, so that no propulsion is allowed in our simplified model. The remaining forces are the gravity and the buoyancy. By artificially decoupling the locomotive gait and the buoyancy effects, we aim to identify the mechanism of stabilisation in the simplest framework. The quasi-static analysis will build the ground of more complex control laws accounting for the hydrodynamic and body inertia. Indeed, the stability of the robot can only be studied near a fixed point, which justifies the need for exploring all the possible stances at equilibrium. When inertia will be considered, this analysis will be helpful to plane real trajectories passing near this sequence of static stances, where the stability properties are already known.

To describe the configuration of the robot, we introduce the configuration space C = G × S ( [START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF], [START_REF] Kelly | Geometric phases and robotic locomotion[END_REF]) given by the direct product of the workspace G, which is defined as the space of the head configuration g 0 taking values in the Lie group SE(3), and the shape space S. Any pair (g 0 , q) is an element of the configuration space C = G × S, a subspace of SE(3) × S, which will be used to parametrize the configuration of the robot. In differential geometry, the space SE(3) × S is called a principal fiber bundle' [START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF] formed by the local direct product of a fiber SE(3) with the base manifold S (see Fig. 2). The forces of gravity and buoyancy being invariant by translation in the surface plane and rotations along the vertical direction, the system is unchanged if the head configurations are related by a rigid transformation on the water surface belonging to the group SE [START_REF] Ijspeert | From swimming to walking with a salamander robot driven by a spinal cord model[END_REF]. Therefore, the configuration of the robot's head is given by the quotient space SE(3)/SE [START_REF] Ijspeert | From swimming to walking with a salamander robot driven by a spinal cord model[END_REF], such that the work space reads G E(1)×S 2 . The Euclidean space E(1) stands for net vertical displacements (heaving), while S 2 represents the configuration space of rotations along vector tangent to the water surface (rolling and pitching). For the body shape, each segments having two revolute joints (rolling+flexural), the shape space is defined by its joint coordinates q ∈ S, with S = (S 1 ) 2n , n being the number of segments.

To maintain a gaze control, our objective is to find a shape control law preserving static equilibrium. Such equilibria are reached when the robot is in a local minimum of potential energy determined by its coordinates (g 0 , q) ∈ C [START_REF] Herault | A geometrically exact approach for floating slender bodies with finite deformations[END_REF]. A sequence of static equilibrium can be seen as a trajectory connecting two points in the configuration space C = G × S with the control q ∈ S (the base manifold) [START_REF] Kelly | Geometric phases and robotic locomotion[END_REF]. If g 0 is fixed in G, the set of all the body shapes q with the same equilibrium g 0 corresponds to a differentiable manifold M g0 embedded in the shape space S (see Fig. 2). Hence, the sequence of equilibria can be parametrized by a trajectory belonging to the manifold M g0 . Therefore, our goal is to construct this trajectory on M g0 by finding at each step the local minimum of potential energy. Rather than standard methods of constrained optimization problem embedded in a Euclidean space, we address this problem of optimization directly on the differentiable manifold, and construct this trajectory along "admissible" directions tangent to M g0 [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF]. Geometrically, our strategy consists in moving on the manifold M g0 thanks to the properties of its local tangent spaces and its normal complements (see Fig. 2). This strategy is very similar to the geometrical concept of "connexion", which associates unequivocally each infinitesimal displacement on the base S with a corresponding displacement from one fiber to another [START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF]. However, the force of buoyancy being not left invariant, no connexion can be defined properly [START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF], [START_REF] Kelly | Geometric phases and robotic locomotion[END_REF].

Here, we will rather seek the admissible displacements on M g0 that produces no variation of g 0 ∈ G to remain on the manifold.

To construct this path on the set M g0 , we report two meth- ods of optimization: a first method providing the equilibrium configuration of the head for a given body shape, and a second method to move on M g0 along admissible directions.

The main contribution of our work is to replace the timeintegration, which disappears for quasi-static motion, by a line search strategies based on geometrical methods. For the first method, we use an algorithm of optimization on the Lie group to find the minimum of potential energy.

For the second method, we use a local mapping from the tangent space of the locally smooth manifold M g0 to elements on the manifold itself thanks to the projection-like retraction associated with a control system [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF]. This paper is the continuation of our previous works on stability of continuous snakes [START_REF] Herault | Standing on the Water: Stability Mechanisms of Snakes on Free Surface[END_REF]- [START_REF] Herault | A geometrically exact approach for floating slender bodies with finite deformations[END_REF] based on the Cosserat beam model, an approach which is now extended to our new discrete robot. The paper is organized as follows: Section II presents briefly the features of the robot and its geometrical model. In section III, the geometrical algorithm to compute the direct and inverse problem are reported. Section IV presents the results of these algorithms, and Section V discusses and concludes on the contributions of this work.

II. THE ROBOT AND ITS GEOMETRICAL MODEL

The bio-inspired robot NATRIX. The robot is similar to other planar snake-like robots [START_REF] Hirose | Snake-like robots[END_REF], [START_REF] Ijspeert | From swimming to walking with a salamander robot driven by a spinal cord model[END_REF], [START_REF] Liljebäck | Eelume: A flexible and subsea resident IMR vehicle[END_REF] with serial revolute flexural joints. To control independently the swimming and the stability, an additional servomotor is fixed to each segment connecting two revolute joints, so that it can rotate the external shell of the module. This new local rolling motion is thus independent of the local bending, and preserves the planarity of the robot, as illustrated on Figs. 1 and3. A similar concept has been recently used for terrestrial locomotion based on Archimedes' Screw [START_REF] Schreiber | ARCSnake: An Archimedes' Screw-Propelled, Reconfigurable Serpentine Robot for Complex Environments[END_REF], but it has never been proposed for aquatic snake-like robot. This robot is named NATRIX, and its design is detailed in an upcoming paper. The averaged robot's density is given by ρ rob = 312 ± 2 kg.m -3 . Each module has a triangular cross-section with rounded corners (see Figs. 1 and5) to reproduce the snake morphology with a surface S rob = 11.3 × 10 3 mm 2 . The advantage of this triangular geometry is to decrease the height of center of Geometrical model. To describe the orientation and position of our robot, a moving frame F 0 is attached to the head. Its orientation and position are defined relatively to the water interface thanks to a Galilean frame F G placed onto the surface (see Fig. 3). The configuration of F 0 relative to F G is given by the matrix g 0 ∈ R 4×4 , the homogeneous representation of its transformation in SE(3), characterized by the rotation matrix R 0 ∈ SO(3) and the position vector p 0 ∈ E(3) [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF]. The vector q defining the body shape is composed of the following joint coordinates (see Fig. 3): {(q f ) k } [1,n] for the angles of the flexural joint, and {(q r ) k } [1,n] for the angles of the rolling motion. To locate each module, a frame F j is attached to the j-th flexural joint, the third axis being the axis of rotation. The configuration g j relatively to the water frame F G is given by

g j = g 0 0 g 1 . . . j-1 g j , (1)
with j-1 g j the relative configuration of the j-th joint frame with respect to the (j -1)-th joint frame. The configuration j-1 g j is defined by the position vector j-1 p j T = (L, 0, 0)

and the rotation matrix

j-1 R j =   cos ((q f ) j ) -sin ((q f ) j ) 0 sin ((q f ) j ) cos ((q f ) j ) 0 0 0 1   . (2)
Then, the configuration (g shell ) j of the shell of j-th module is deduced from the frame F j by a longitudinal translation to the center of the module, followed by a longitudinal rotation given by the rolling angle with

(g shell ) j = g j R shell (j) p shell (j) 0 1 , (3)
with the position vector (p shell ) T = (L/2, 0, 0), and the rotation matrix

R shell (j) =   1 0 0 0 cos((q r ) j ) -sin((q r ) j ) 0 sin((q r ) j ) cos((q r ) j )   . ( 4
)
Hydrostatic forces. In a previous paper [START_REF] Herault | A geometrically exact approach for floating slender bodies with finite deformations[END_REF], we have introduced a model to compute the immersed volume of a solid body based on its configuration g relative to the surface. We review here briefly the method. For one given cross-section of a module, the air-water interface is given by the parametric equation p z + R 32 Q 2 + R 33 Q 3 = 0, defining the waterline, with Q 1 , Q 2 , Q 3 the coordinates of a point in the local section frame (section 3.2 in [START_REF] Herault | A geometrically exact approach for floating slender bodies with finite deformations[END_REF]), and R i,j the component of the rotation matrix of the considered cross-section. According to the different cases of intersection between the water and the cross-section of the module, the immersed surface and its barycenter can be numerically computed for a triangular cross-section. Then, we compute the immersed volume of one module by integrating numerically the contribution of all the immersed surfaces. Finally, we can evaluate numerically the net potential energy U (g 0 , q 0 ) of the robot by summing the potential energy of each module, and the net wrench W, which has only three non-zero components in the frame F G : the vertical force F z and the transversal torque Γ x and Γ y .

We write WT = (F z , Γ x , Γ y ) the reduced wrench with these three components.

III. THE ALGORITHMS

Now, we introduce two methods to find and then move on the manifold M g0 corresponding to all the elements q ∈ S satisfying W(g 0 , q) = 0 for a given g 0 . First, we seek the root of W for a fixed body shape q thanks to a sequence of head configurations {(g 0 ) n } indexed by n so that Wn → 0 when n increases. To find the root, a gradient descent method adapted to the geometry of the Lie group SE(3) is used. We define an update formula to pass from a head configuration (g 0 ) n to (g 0 ) n+1 , such that the energy decrease satisfies the Armijo rule [START_REF] Nocedal | Numerical optimization[END_REF]. To preserve the structure of the Lie group SE(3), the update formula is a left matrix multiplications, and reads: (g 0 ) n+1 = A n (g 0 ) n with the matrix A n ∈ SE(3) given by

A n = exp t n ξn , (5)
The relative transformation A n is defined from an unitary direction ξ n = W/||W||, and a step size t n determined by the Armijo rule [START_REF] Nocedal | Numerical optimization[END_REF]. Then, the transformation A n is computed thank to the exponential map, mapping an element from the Lie algebra ξn ∈ se(3) to an element of the group. After some iterations of the descent methods approaching the minimum, we accelerate the code via a second-order method based on a Netwon-Raphson algorithm. The reported method provides the equilibrium state, and gives the starting point in the configuration space G to move on M g0 .

Then, we want to move on the manifold M g0 by modifying the body shape while preserving the head configuration g 0 ∈ G. To address this problem, we introduce the concept of retraction, which is used in differential geometry to construct a path from the elements of the tangent space of a smooth Riemannian manifold [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF]. The simplest way to move on a manifold is to follow a geodesic curve on a manifold endowed by a Riemannian structure. However, finding these geodesic curves for an arbitrary manifold can be computationally expensive providing this update exists. An alternative method consists in constructing an approximation of the geodesic curve, at least at order o(t). This method called a retraction is a geometrical construction connecting two infinitesimally closed points of the manifold, so that it shares the two leading order properties of the Taylor expansion of a geodesic curve. A mapping u → R(q, u) from a tangent space T q M g0 into M g0 , with q and R(q, u) ∈ M g0 , is a retraction if and only if the curve γ(q, u, t) = R(q, tu) with t ∈ R + , satisfies γ(q, u, 0) = q and γ(q, u, 0) = u.

(

) 6
The reader can find in [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF] a more formal definition of the retraction. Here, we use a projection-like retraction (also called orthographic retraction [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF]) illustrated in Fig. 4 that is defined by

R(q, tu) = q + tu + h [tu] , (7)
with u ∈ T q M g0 and h a orthographic retractor belonging to the orthogonal complement N q M g0 of T q M g0 , which are subspaces of S. It is worth noting that the update defined by Eq. 7 doesn't require a manifold endowed with a Riemannian structure so that no metric is necessary in our approach, only an inner product given by the dot product of R 2n .

To identify the local tangent space and its complement, we compute the exact differential d W of W. Then, we state that for any variation on the set of equilibrium configurations, the variations of body shape δq and head configuration δ ψ (defined below) have to satisfy d W = 0, which reads

D ψ W δ ψ + D q W δq = 0, (8)
with the jacobian matrices D ψ W ∈ R 3×3 and Dq W ∈ R 3×2n . Equation 7 plays the role of a connection on a principal fiber bundel. However, the jacobian matrices depend on the configuration g 0 , i.e. Eq. 8 is not left invariant, and doesn't define a connexion. The matrix D ψ W is called the hydrostatic stiffness matrix, and it is computed from the variation of the head configuration δg parametrized by δ ψ ∈ se(3) so that δg 0 = δ ψg 0 . We introduce the mapping δ ψT = (δz, δθ x , δθ y ) → δψ T = (0, 0, δz, δθ x , δθ y , 0), with δψ the coordinates of the twist. This transformation allows to map elements from R 3 to se(3). For a stable equilibrium configuration, the opposite of the eigenvalues λ i of D ψ W give the positive hydrodynamic stiffness characterizing the local stability of the equilibrium [START_REF] Fossen | Handbook of marine craft hydrodynamics and motion control[END_REF]. For a straight body (no bending), the eigenvalue associated with the eigenvector in the longitudinal direction e 1 characterizes the stability of the rolling mode, while the two others eigenvectors correspond to two modes of oscillations coupling the heaving and the pitching. If one of the eigenvalues becomes positive, which corresponds to a negative stiffness, the equilibrium is unstable. Note that this matrix D ψ W is also used in the Netwon-Raphson method reported in the optimization on the Lie group. The matrix Dq W being rectangular, we perform a singular value decomposition (SVD) to compute its kernell, which is generated by the basis vectors {v } spanning the local tangent space T q M g0 . Then, the jacobian matrix Dq W can be reduced to

D q W ⊥ = U ⊥ ΣV T ⊥ ( 9
)
with Σ the non-zero diagonal matrix obtained by the SVD, and V T ⊥ the matrix with the basis vectors {v ⊥ }, spanning the orthogonal complement N q M g0 of T q M g0 in the subspace of the joint coordinates. In our study, the rank of the matrix D q W ⊥ is always 3. A displacement tu on the tangent space produces a residual W (g 0 , q + tu) of order o(t). The orthographic retractor h is thus define to correct the small increase of wrench required to compensate any head variation, i.e. d ψ = 0. To compensate this drift from the manifold M g0 , the orthographic retractor reads

h [tu] = -α D q W + ⊥ W (g 0 , q + tu) (10)
with D q W + ⊥ the pseudo-inverse. The parameter α(u, t) is a free real parameter, which is in practice close to one for t 1, defining the step size along the admissible direction to come back to the manifold [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF]. By construction, the retractor h [tu] belongs to N q M g0 and is of order O(t 2 ), so that the combination of Eqs. 7 and 10 is a retraction satisfying Eqs. 6. Now, the path is constructed by sequential updates, such that the n-th update is characterized by its parallel step direction u n and step size t n . To reach the goal body shape q goal , the unitary step direction u n is constructed from the projection of the update goal direction q goal -q n onto the tangent space, so that

ũn = V V T q goal -q n , u n = ũn |ũ n | (11)
with V the rectangular matrix with the basis vectors {v }.

The step size t n is a slowly varying function of n, such that the system converges slowly to the goal.

IV. THE SIMULATIONS

Thanks to the reported algorithms, we will investigate to what extent the buoyancy and the body deformations can be used to stabilize the head of the snake-like robot. Initially, the model of our robot is composed of 7 bodies: one head plus 6 modules with 6 flexuraland rolling angles. First, we study the effect of the volumetric mass on the local stability of a straight robot (q = 0) at equilibrium. In Fig. 5, we have reported the evolution of the largest eigenvalue λ 3 of the hydrostatic stiffness matrix, which is associated with a pure net rolling motion, as a function of the averaged volumetric mass of the robot. As expected, the level of immersion, illustrated by the embedded figures, goes up when the density increases. The eigenvalue, which is initially negative, also increases with the density, until the sign of the eigenvalue changes at ρ c = 497kg.m -3 . At this critical point, the equilibrium becomes unstable: any infinitesimal rotation along the longitudinal axis will overturn the robot. This result confirms that the robot NATRIX has been designed to remain stable with a low density ρ = 312kg.m -3 significantly smaller than ρ c . This stability can also be improved by the anguilliform swimming gait. We compute a sequence of equilibrium states by moving on the manifold M g0 via intermediate goal body shapes (see left Fig. 6) defined by a locomotive gait. The robot's density is fixed to 0.3. The values of the six flexural angles for the intermediate goal shapes are reported on table I for an half period, the setpoints for the other half are given by the opposite values. The rolling angles are set to zero. These goal shapes are used in Eq. 11 to define the admissible direction on M g0 . In Fig. 6 (right), we have reported the evolution of the flexural (thick curve) and rolling (dashed curve) angles of the second module as a function of the number of iterations, which can vary to reach the intermediate goal body shape. The flexural angle displays a sine-like function while the rolling angle remains null, as the other rolling angles. Indeed, the algorithm selects a path with solely flexions to preserve the stability. Then, the evolution of the eigenvalues of the hydrodynamic stiffness are reported on Fig. 6 (left) as a function of number of iterations of the process. The first iterations (up to 300) correspond to the transition from the straight configuration to the first goal body shape. All the eigenvalues remain negative, which implies that the quasi-static motion is always stable. However, their evolution differs significantly. The eigenvalues λ 1 and λ 2 (dashed curves) correspond to modes combining heaving and pitching motions. Their values are much smaller than the third eigenvalues, between one and two order of magnitude, and vary relatively slightly during the motion. The third eigenvalue λ 3 , which corresponds mainly to a rolling motion, varies significantly during the gait. First, the stiffness is increased by a factor 8 when the body shape departs from the straight configuration, and reaches the first goal body shape. The smallest eigenvalue, which represents the most stable state, is reached at t = T /4 modulo T /2, and can be divided by two when the system is in the least stable state during the gait. This simulation demonstrates that the robot is particularly sensitive to the rolling motion, and its stability can significantly change during a gait.

t (q 1 ) f (q 2 ) f (q 3 ) f (q 4 ) f (q 5 ) f (q 6 ) f 0 30 o -
We have tried to implement other goal body shapes with different setpoints for the rolling angles. The results were initially surprising: the rolling angles remain null regardless of the setpoints. The explanations is given by the physics of the system. When one module rolls, its immersed volume decreases, such that the equilibrium is no more satisfied. To balance this effect, the altitude of the head must decrease to dive and recover the immersed volume. However, this motion is forbidden in our algorithm: g 0 can't be changed. Therefore, the distributed rolling system can't be used when the pitching and the heaving of the head is not allowed.

To overcome this problem, we have added to our robot a neck between the head and the first body module. This neck has two pitch angles to compensate the pitching and the heaving of the body. This new robot model is composed of 8 bodies: one head, one neck and 6 modules with 5 flexural joints (one is replaced by the pitch joint) and 6 rolling joints. In Fig. 7, the eigenvalue λ 3 associated with the rolling motion are reported for three different shape control strategies. The first one (black curve) corresponds to a pure flexural motion without rolling angle (q r = 0), which is the reference curve. Then, we have implemented a rolling motion inspired from aquatic snake, such that each shell rotates inward when the module departs from the straight configuration (see the small schema in the right corner of Fig. 7). This control law is surprisingly not effective since the eigenvalue (dot-dashed curves) remains above the reference curve (black curve). Finally, the last control law is an alternated rolling motion between each segment, such that the rolling angles changes its sign relatively to the previous rolling angle and also between each goal shape. Here, we observe a small improvement of the stability (the dotted curve having the smallest mean value) relatively to the pure flexural body motion. V. CONCLUSIONS The work presented in this paper represents the first steps in a program of research aiming to control the stability of snake-like robots swimming on water surface. First, we have introduced a new robot, called NATRIX, with an additional degree of freedom allowing to rotate the shell of each body module. To predict the stability of the robot, we have reported a method to deform the body shape while preserving the head configuration. Our study has demonstrated that the net rolling motion is the most dangerous degree of freedom for slender snake-like robot. Our simulations shows that the stability can be improved by the undulatory gait with a significant modulation of the hydrodynamic stiffness. A neck with two additional pitch angle is required to use the distributed rolling motion. Our simulations show that this simple additional degree of freedom can be easily controlled, and improve the stability of the robot. A step further would consist in seeking the optimal control of the body shape to provide a stable and controllable robot on water surface. Our future work will also consist in considering the dynamical aspects of the stability, and investigate to what extent the quasi-static stability can be used. Moreover, we are working on the instrumentation of the robot to compare these results with the reported results.

Fig. 1 .

 1 Fig. 1. The robot NATRIX with its distributed controlled rolling motions.

Fig. 2 .

 2 Fig. 2. Geometrical representation of the process of head stabilization in the principal fiber bundle SE(3) × S. A sequence of equilibrium states with different body shapes q n can be represented by any trajectory on the manifold Mg 0 , subset of the base state S, while the head configuration g 0 remains unchanged in G = SE(3)/SE(2).

Fig. 3 .

 3 Fig. 3. Parametrization of the robot configuration

Fig. 4 .

 4 Fig. 4. Illustration of the retraction

Fig. 5 .

 5 Fig. 5. Largest eigenvalue of the hydrodynamic stiffness matrix, corresponding to a rolling motion, as a function of the volumetric mass ρ for a straight body q = 0.

Fig. 6 .

 6 Fig. 6. (Left) Eigenvalues of the hydrodynamic stiffness matrix as a function of the iterations, which are required to undergo one period of a undulatory gait, illustrated by the intermediate body shapes. (Right) Evolution of the flexural and rolling joint angles of the second module as a function of the iterations.

Fig. 7 .

 7 Fig. 7. Eigenvalue λ 3 of the hydrodynamic stiffness matrix as a function of the iterations for three different cases (see text).

TABLE I INTERMEDIATE

 I FLEXURAL ANGLES FOR THE GOAL BODY SHAPE

IMT Atlantique, LS2N, UMR CNRS 6004, F-44307 Nantes, France, johann.herault@imt-atlantique.fr

ACKNOWLEDGMENT

This work was supported by the French National Research Agency ANR (grant no. ANR-2019-CE33-0004-01) and La Region Pays de la Loire (program " Étoile montante").


Kolade Augustine Oyekan

Maarten Van De Put

Sabyasachi Tiwari

Carole Rossi

Alain Esteve

William Vandenberghe
email: william.vandenberghe@utdallas.edu

Re-examining the Role of Subsurface Oxygen Vacancies in the Dissociation of H 2 O Molecules on Anatase TiO 2 Authors

Keywords: DFT, Anatase TiO 2, Oxygen vacancy, Water dissociation, Full coverage, Hydrogen

Titanium (IV) oxide is a reference photocatalytic material for a number of chemical processes of interest in the energy transition, including the reduction of water into hydrogen gas. Density functional theory calculations have previously revealed that the decomposition into hydrogen and hydroxyl ions could be catalyzed by a subsurface oxygen vacancy (V O ). However, these theoretical investigations have been carried out with models considering only isolated water molecules. We investigate the role of the subsurface V O at different levels of coverage up to full monolayer (ML) coverage. We observe that migration of the vacancy concomitantly to water decomposition is mandatory to obtain an exothermic pathway. Furthermore, the migration of the V O only occurs in the presence of an under-coordinated surface oxygen atom. As coverage and decomposition increase, surface oxygen atoms are saturated which inhibit further migration of the subsurface V O ; hence the V O shows no catalytic effect on the splitting of water into hydrogen and hydroxyl ions at full water coverage.

Introduction

Titanium (IV) oxide (TiO 2 ) is arguably one of the most investigated semiconductor materials due to its exquisite physical and chemical properties and consequent numerous applications [START_REF] Diebold | The surface science of titanium dioxide[END_REF]. Among other metal oxides, TiO 2 's chemical stability, low toxicity, environmentally friendly nature, and reusability, have led to its widely accepted consideration as a top leading photocatalyst [START_REF] Hoffmann | Environmental Applications of Semiconductor Photocatalysis[END_REF], particularly since the seminal work by Honda and Fukushima in 1973 [START_REF] Fujishima | Electrochemical Photolysis of Water at a Semiconductor Electrode[END_REF], which paved the way to photoelectrochemical materials and systems for hydrogen production through the splitting of water molecules. Since then, much effort has been focused on overcoming the TiO 2 weaknesses, such as the high electron/hole recombination rate and its optical absorption limited to the UV region, in both theoretical (mostly density functional theory, DFT) and experimental communities [START_REF] Mendoza-Diaz | On the UV-Visible Light Synergetic Mechanisms in Au/TiO2 Hybrid Model Nanostructures Achieving Photoreduction of Water[END_REF][START_REF] An | Hierarchical Nanotubular Anatase/Rutile/TiO2(B) Heterophase Junction with Oxygen Vacancies for Enhanced Photocatalytic H2 Production[END_REF][START_REF] Ismail | Photochemical splitting of water for hydrogen production by photocatalysis: A review[END_REF]. While several experimental research routes including doping [START_REF] Asahi | Visible-Light Photocatalysis in Nitrogen-Doped Titanium Oxides[END_REF][START_REF] Liu | One-step fabrication of N-doped mesoporous TiO2 nanoparticles by self-assembly for photocatalytic water splitting under visible light[END_REF][START_REF] Ohno | Photocatalytic Activity of S-doped TiO2 Photocatalyst under Visible Light[END_REF], heterojunctions with other semiconductors [START_REF] Coy | Exploring the effect of BN and B-N bridges on the photocatalytic performance of semiconductor heterojunctions: Enhancing carrier transfer mechanism[END_REF], and embedding of plasmonic nanoparticles [START_REF] Mendoza-Diaz | On the UV-Visible Light Synergetic Mechanisms in Au/TiO2 Hybrid Model Nanostructures Achieving Photoreduction of Water[END_REF][START_REF] Cure | A triptych photocatalyst based on the Co-Integration of Ag nanoparticles and carbo-benzene dye into a TiO2 thin film[END_REF] have been pursued, DFT calculations have concentrated on surface chemistry for a number of applications and associated molecular decomposition of probe molecules such as methanol, ethanol, and water [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Lang | Interplay between Methanol and Anatase TiO2(101) Surface: The Effect of Subsurface Oxygen Vacancy[END_REF][START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF][START_REF] Vittadini | Structure and Energetics of Water Adsorbed at TiO 2 Anatase \(101\) and \(001\) Surfaces[END_REF], which is the focus of the present article. TiO 2 has three stable polymorphs: rutile, anatase, and brookite. Rutile and anatase are widely studied for catalytic applications, but anatase has been shown to be more photo-reactive, and it is also the most stable form found in nanocrystals [START_REF] Vittadini | Chemistry of and on TiO2-anatase surfaces by DFT calculations: a partial review[END_REF][START_REF] Lazzeri | Structure and energetics of stoichiometric TiO 2 anatase surfaces[END_REF] . For these reasons, anatase is the most utilized form of TiO 2 in photocatalytic applications, and consequently in DFT-based investigations. The (101) surface of anatase was established to be the most stable one [START_REF] Vittadini | Chemistry of and on TiO2-anatase surfaces by DFT calculations: a partial review[END_REF][START_REF] Lazzeri | Structure and energetics of stoichiometric TiO 2 anatase surfaces[END_REF][START_REF] Arrouvel | Effects of morphology on surface hydroxyl concentration: a DFT comparison of anatase-TiO2 and γ-alumina catalytic supports[END_REF], and this has led to a number of DFT investigations dedicated to the water splitting process [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF][START_REF] Valdés | Water adsorption and dissociation on gold catalysts supported on anatase-TiO2(101)[END_REF].

Results from the literature indicate that the limiting step in the dissociation of a single water molecule on the anatase (101) surface lies on the ability of an adsorbed water molecule to decompose into an adsorbed hydroxyl group and a hydrogen ion in solution [START_REF] Li | Mechanism and Activity of Photocatalytic Oxygen Evolution on Titania Anatase in Aqueous Surroundings[END_REF]. In the calculations by Aschauer and co-workers, the activation barrier of this step for a single water molecule on a pristine anatase (101) surface is 0.52 eV, and the overall endothermic reaction energy is 0.38 eV, making it thermodynamically unfavorable. [START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF]. Interestingly, studies have shown that the dissociation process can be enhanced by the presence of an oxygen vacancy in the anatase slab [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF], and many studies have reported that a subsurface vacancy has a lower formation energy than a surface vacancy [START_REF] He | Evidence for the Predominance of Subsurface Defects on Reduced Anatase TiO 2 ( 101 )[END_REF][START_REF] Cheng | Energetics and diffusion of intrinsic surface and subsurface defects on anatase TiO2(101)[END_REF]. Rather than a surface vacancy, a defective slab with a subsurface oxygen vacancy led to a reduced kinetic barrier (by 0.26 eV), relative to a stoichiometric slab, but the water dissociation process was shown to remain thermodynamically unfavorable (by 0.32 eV). Experiments at low temperature and pressure on a defective slab with subsurface defects also indicate higher binding energies for water on a defective slab, compared to a pristine slab [START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF]. In contrast, some results from the literature, including the present study, show that not only is the kinetic barrier reduced, but the overall reaction becomes thermodynamically favorable in the presence of a subsurface oxygen vacancy [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF]. However, it is well documented that molecular water adsorption is exothermic [START_REF] Vittadini | Structure and Energetics of Water Adsorbed at TiO 2 Anatase \(101\) and \(001\) Surfaces[END_REF][START_REF] Valdés | Water adsorption and dissociation on gold catalysts supported on anatase-TiO2(101)[END_REF][START_REF] Tilocca | Vertical and Lateral Order in Adsorbed Water Layers on Anatase TiO2(101)[END_REF][START_REF] He | Local ordering and electronic signatures of submonolayer water on anatase TiO2(101)[END_REF] and that it is thermodynamically favorable for the surface to be completely covered with chemisorbed water molecules in their molecular form. The effect of defective surfaces has been explored to a lesser extent, while being of much technological importance. In particular, DFT studies examining the effect of subsurface V O have only dealt with a single or sparse concentration of water molecules on the anatase (101) surfaces. Importantly, it was reported that the enhanced dissociation in the presence of a subsurface oxygen vacancy is accompanied by the migration of the vacancy towards the surface [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF], leading to its annihilation. As a result, the exact role of the vacancy defect in the photocatalysis of water on the TiO 2 anatase (101) surface requires further investigation. It is also important to investigate the issue of the impact of water coverage on the dissociation chemistry when considering defective TiO 2 .

In this work, we study theoretically water dissociation on defective anatase (101) accounting for full water coverage. Using DFT, we investigate and compare the dissociation of water molecules on a perfect and a defective surface at different levels of coverage: from a minimum coverage of 17% up to 100% (1 ML coverage).

We examine the behavior of the subsurface oxygen vacancy during dissociation. Our results indicate that dissociation of adsorbed water molecules on a defective surface is exothermic up to 67% coverage (2/3 ML). At higher, more realistic, coverage the dissociation becomes endothermic, similar to the endothermic dissociation obtained for a perfect surface at all levels of coverage. Therefore, we thus find that the V O has no catalytic effect on the splitting of water into hydrogen and hydroxyl ions on a defective anatase (101) surface at high coverage.

Method

DFT Calculations

We performed spin-polarized DFT calculations using the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzehof (PBE) [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] exchange-correlation potential as implemented in the Vienna ab-Initio Simulation Package (VASP) [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Ab initio molecular dynamics for liquid metals[END_REF][START_REF] Kresse | Ab initio molecular-dynamics simulation of the liquid-metal--amorphous-semiconductor transition in germanium[END_REF]. Ion-electron interactions were described by projector augmented waves (PAW) pseudopotentials [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF][START_REF] Blöchl | Projector augmented-wave method[END_REF], and the electron wave functions were expanded in terms of plane wave basis sets, with a kinetic energy cutoff of 400 eV. Wavefunctions were optimized until the change in total energy was less than 10 -6 eV. We employed Grimme's DFT-D3 method to include dispersion contributions to the energy values [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF]. Potential energy surfaces to determine the kinetic barrier were calculated using the nudged elastic band (NEB) method [START_REF] Jonsson | Nudged elastic band method for finding minimum energy paths of transitions[END_REF].

The anatase (101) surface supercell is built from a three-layer geometry (see Fig. 1), with a vacuum of 10 Å between slabs along the direction normal to the surface. The 1 x 3 surface supercell measuring 10.376 x 11.256 Å 2 has a total of 108 atoms. The bottom layer of TiO 2 was fixed while all other atoms were allowed to relax in all directions during minimization of the total energy. Due to the large size of the supercell, the Brillouin zone was sampled at the gamma point only. The unit cell from which the surface slab is cut was relaxed using a kinetic energy cutoff of 520 eV, and the lattice parameters obtained are: a = b = 3.8073 Å, c = 9.7318 Å, which are in good agreement with values from other computational and experimental studies [START_REF] Liu | Electronic and elastic properties of the tetragonal anatase TiO2 structure from first principle calculation[END_REF][START_REF] Burdett | Structural-electronic relationships in inorganic solids: powder neutron diffraction studies of the rutile and anatase polymorphs of titanium dioxide at 15 and 295 K[END_REF]. The single oxygen vacancy defective surface slab was obtained by removing the oxygen atom labelled O sub in Fig. 1. Previous studies, both experimental and theoretical, have established that the subsurface vacancy created by removing O sub has a lower formation energy than the most stable surface vacancy created by removing O 2c , or other subsurface oxygen sites [START_REF] Cheng | Energetics and diffusion of intrinsic surface and subsurface defects on anatase TiO2(101)[END_REF][START_REF] Wang | Interplay between invasive single atom Pt and native oxygen vacancy in anatase TiO2(1 0 1) surface: A theoretical study[END_REF][START_REF] Setvin | Aggregation and electronically induced migration of oxygen vacancies in TiO 2 anatase[END_REF].

Adsorption Energy Calculations

An incremental energy needed for the molecular adsorption of the n th water molecule on each surface, is defined as: [START_REF] Diebold | The surface science of titanium dioxide[END_REF] where n = 6 means full coverage. and are the total energies of a slab supercell containing n and (n-1) adsorbed water molecules, respectively. is the total energy of an isolated water molecule, as being positioned at an infinite position from the surface, in vacuum. A negative value of indicates an exothermic process, while a positive value indicates an endothermic process.

The molecular adsorption energy per molecule is given by: where is the total energy of a surface slab, and n is the number of water molecules adsorbed on the surface. The dissociative adsorption energy per molecule is:

where is the total energy of a slab with n dissociated water molecules on the (1 x 3) surface.

Results and Discussion

Single water molecule interaction

In table 1, we list the molecular adsorption and dissociative adsorption energies we calculate for a perfect and a defective surface. The interaction of a single water molecule with a perfect and a defective surface is a reference system, which has been discussed in several DFT-based papers [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF]. While the energy profile of the water molecule decomposition is largely endothermic on the perfect surface, the presence of a subsurface oxygen vacancy (as defined in Figure 1) has led to contrasting results, showing slight endothermic [START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF] or clear exothermic profiles [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF], as seen in Table 1. The issue behind the differences between the four results presented in table 1 lies in the stability of the subsurface oxygen vacancy. In Ref. [START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF], the oxygen vacancy remains stable, leading to an endothermic reaction pathway for dissociation of H 2 O. In Fig. 3, we present the energy pathway for the entire dissociation process of an adsorbed water molecule, accompanied with the migration of the V O to the surface from the structure in Fig. 2a to the structure in Fig. The calculated values of incremental adsorption energy on both surfaces, up to 1 ML coverage, are shown in table 2. The results show that 1 ML coverage with adsorbed water molecules is possible on both stoichiometric and defective surfaces since the successive adsorption energies at all levels of coverage are negative. The adsorption energy on the stoichiometric surface ranges from -0.64 eV to -0.91 eV, while on the defective surface, the values vary between -0.80 eV and -0.99 eV. Overall, molecular adsorption is more exothermic on the defective surface than on the perfect surface. Our results for the stoichiometric surface agree with previously published results for a similar (1 x 3) surface slab showing successive adsorption energies between -0.70 eV and -0.85 eV [START_REF] Valdés | Water adsorption and dissociation on gold catalysts supported on anatase-TiO2(101)[END_REF], and also agree with experimental results [START_REF] He | Local ordering and electronic signatures of submonolayer water on anatase TiO2(101)[END_REF].

In Table 3, we show the effect of the subsurface V O on the energy required for dissociation of adsorbed water molecules on the slab surface at different levels of coverage, with a comparison to the adsorption on the perfect surface. To compare, we calculate the adsorption energy per molecule for each type of adsorption, using equation [START_REF] Fujishima | Electrochemical Photolysis of Water at a Semiconductor Electrode[END_REF]. For a stoichiometric surface, molecular adsorption is always energetically more favorable than dissociative adsorption at all levels of coverage with water molecules. While non-dissociative adsorption energy is stable with increasing coverage, dissociative adsorption energy slightly increases (more than 0.1 eV gain). For a defective surface, the energy released on dissociative adsorption of the first water molecule on the surface is 2.40 eV, relative to the infinitely separated water molecule and slab, compared to the 0.99 eV released on molecular adsorption. In other words, 1.41 eV is released when an isolated water molecule adsorbed on the surface dissociates. A similar exothermic profile persists on the defective surface for the first four adsorbed water molecules, which corresponds to 2/3 ML coverage. At 2/3 ML coverage (coverage n = 4), the energy released per adsorbed water molecule drops down to 0.09 eV. Beyond this level of coverage, nondissociative molecular adsorption becomes more favorable, similar to the dissociation on perfect surfaces. When coverage increases to 5/6 and 6/6 water molecules, the energy required for dissociation of each adsorbed water molecule on both surfaces is more than 0.2 eV, clearly illustrating the switch from an exothermic to an endothermic profile as coverage increases beyond 2/3 coverage. -

No of water molecules stoichiometric defective Fig. 6 shows the results of charge density difference analysis, which illustrates the distribution of the excess charge in the system due to the oxygen vacancy. The position of the vacancy is colored green, indicating a loss of charge while the regions gaining charge, colored yellow, are around Ti atoms. The significant reduction observed in the dissociative adsorption energy can be attributed to the charge distribution pattern seen in Fig. 6. Out of the six surface Ti 5c atoms to which water molecules are attached, the only Ti 5c with a significant concentration of excess charge is the one to which the first water molecule is adsorbed. The We further introduced a second subsurface V O to investigate the effect of vacancy concentration or the formation of a subsurface vacancy cluster, which has been observed experimentally [START_REF] Setvin | Aggregation and electronically induced migration of oxygen vacancies in TiO 2 anatase[END_REF], on the dissociation of adsorbed water at 5/6 and 1ML coverage. We considered two possible cases: separated vacancies and a cluster of vacancies. The results, as presented in Table S1 and Figs S2-S4 show that the dissociative adsorption remains endothermic relative to molecular adsorption at these levels of coverage. The energy required to dissociate a molecule at 5/6 ML coverage reduces for both slabs with two vacancies. For the slab with separated vacancies, the final structure is further stabilized because a hydrogen bond forms between the hydrogen atom of a dissociated water molecule and a neighboring oxygen atom as seen in Fig S3(c), and the energy required for dissociation reduces from 0.24 eV in the single V O structure to 0.09 eV. In the slab with a cluster, the required energy for dissociation is smaller: 0.05 eV, due to the relaxation of the two-fold coordinated (O 2c ) oxygen atom on the surface, as illustrated in Fig S4 (c). A recent study on the more reactive but less stable (001) surface of anatase also shows that water splitting remains endothermic in the presence of two surface oxygen vacancies [START_REF] Mete | Elucidating the Barriers on Direct Water Splitting: Key Role of Oxygen Vacancy Density and Coordination over PbTiO3 and TiO2[END_REF], but the comparison between that investigation and this work is limited because of the differences between the facets of anatase studied and the positions of the oxygen vacancies considered. The diagrams of the relaxed structures in Figs S3 and S4 are similar to those obtained for a single V O . There is a partial migration of one V O towards the surface at 5/6 ML coverage, while no significant V O migration is seen at 1 ML coverage for both types of double V O .

Previous investigations with only one adsorbed water molecule have concluded that the presence of a subsurface vacancy enhances dissociation by either reducing the energy barrier required for dissociation while the dissociation becomes less endothermic [START_REF] Aschauer | Influence of Subsurface Defects on the Surface Reactivity of TiO2: Water on Anatase (101)[END_REF], or changing the dissociation step to an exothermic process accompanied by a reduced dissociation barrier [START_REF] Fisicaro | Wet Environment Effects for Ethanol and Water Adsorption on Anatase TiO2 (101) Surfaces[END_REF][START_REF] Li | Interplay between Water and TiO 2 Anatase (101) Surface with Subsurface Oxygen Vacancy[END_REF]. But the described change in the energy profile during dissociation is rather local, and coverage dependent. Our results indicate that at 1 ML full coverage and in the presence of either one or two subsurface V O s, dissociation of adsorbed water molecules is still endothermic, and the vacancy has no catalytic effect on the splitting of water molecules adsorbed on the defective surface. When O C is fully coordinated, as seen in Figs. 7(e) and (f), the vacancy remains in a subsurface position and the splitting of the adsorbed water molecule is endothermic, making it unlikely to take place.

We have investigated the role of a subsurface oxygen vacancy in the splitting of chemisorbed water molecules on the anatase (101) surface at different levels of coverage, up to 1 ML coverage. Up to 2/3 ML coverage, the subsurface V O enhanced the splitting chemistry, making dissociative adsorption more thermodynamically favorable than molecular adsorption. In all instances, this enhancement was accompanied by a migration of the V O towards the surface, and the migration was made possible by the presence of an under-coordinated surface oxygen atom near the vacancy. In the absence of such a surface atom, and especially at full ML coverage, molecular adsorption is more thermodynamically favorable than dissociative adsorption on the defective surface with a subsurface V O , similar to the results obtained for a stoichiometric surface.

Our findings indicate that contrary to previous findings, subsurface V O are unlikely to catalyze water splitting.

While a subsurface V O enhances the photo-splitting on a surface with a defect, it only does this when not accounting for water coverage.

2. 2

 2 Anatase TiO 2 model surfaces The stoichiometric anatase (101) surface slab is shown in Fig. 1. It has a saw tooth-like morphology, with rows of fully coordinated (Ti 6c ) and under coordinated (Ti 5c ) Ti atoms along the [010] direction. These rows of Ti atoms are joined by three-fold coordinated (O 3c ) and two-fold coordinated (O 2c ) oxygen atoms. Water molecules are most easily adsorbed on the Ti 5c atoms, and for the 1 x 3 supercell used in the calculations, there are six such Ti 5c sites (two rows of 3 Ti 5c ).

Figure 1 :

 1 Figure 1: Side view of the 1 x 3 anatase (101) slab. Titanium atoms are cyan while oxygen atoms are red. O sub is the oxygen atom removed to create the subsurface vacancy in the defective surface. Atoms underneath the dashed line are held fixed during minimizations of the total energy.

Figs 2 (

 2 Figs 2(a) -(c) show the molecular adsorption, metastable state of dissociative adsorption, and final state of dissociative adsorption of a water molecule on a defective surface, respectively. The oxygen atoms labeled 1, 2, and 3 in Fig 2(a) are involved the migration of the subsurface V O towards the surface, and the arrows in Fig 2(b) show the movement of each atom during the migration. In agreement with Refs. [12 and 15], these figures show that the subsurface V O migrates towards the surface as the water molecule dissociates,according to the following process: (i) after adsorption, reconstruction around the vacancy (V O ) is operating with the oxygen atom labeled 1 moving towards the vacancy site, in a bridge site situation between two Ti atoms (Fig.2(a)); (ii) as the water molecule dissociates, 1 relaxes back downwards to its stoichiometric crystal position, while oxygen 2 migrates to occupy the initial V O site. At this stage the vacancy is located immediately below the surface, i.e. in the third oxygen layer; 3 then moves downwards, towards the initial position occupied by oxygen 2 before its migration, oxygen 3 stabilizes as a bridge site, extending the defect site to two oxygen vacancy sites on the surface and immediate subsurface from either side of the oxygen 3 atom, as seen in Fig.2(b). The structure in Fig.2(b) is a metastable state and after overcoming another activation barrier the hydroxyl group attached to the surface migrates to annihilate the vacancy defect site as described in (ii) results in a final stable state as shown in Fig.2c.

2c.

  The limiting step is the splitting of the adsorbed water molecule into a hydrogen atom, which binds to an oxygen atom on the surface, and a hydroxyl group which remains attached to the Ti atom. The energy barrier of 0.33 eV between the metastable and final states in Figs 2(b) and (c) respectively is the energy needed to move the hydroxyl group towards the surface V O shown in Fig. 2(b). The entire process is exothermic, releasing a total of 1.41 eV, and has a low barrier of 0.35 eV.

Figure 2 :Fig. 3

 23 Figure 2: Diagrams of (a) molecular adsorption (b) metastable state of dissociative adsorption, and (c) final state of dissociative adsorption of a single water molecule on the defective surface. Oxygen atoms associated with vacancy migration are colored green and labeled 1, 2, and 3, while position of V O is shown in yellow. Cyan: Ti atom; red: O atom; white: H atom

Figure 4 :

 4 Figure 4: Ball and stick model of the stoichiometric TiO 2 anatase (101) slab, covered with 1 ML of adsorbed water molecules. Titanium atoms are cyan, oxygen atoms are red, and hydrogen atoms are white.

Figure 5 :

 5 Figure 5: Energy required (per molecule) to dissociate an adsorbed water molecule on both stoichiometric and defective surfaces, at different levels of coverage.

  excess charge facilitates a stronger bond between the Ti 5c and the oxygen atom of the water molecule. The resulting weakening of O-H bonds thereby enhances the dissociation of the water molecule into H and OH species. The system is stabilized by the dissociation of the water molecule, leading to the relatively large energy gain of 1.41 eV upon dissociation.

Figure 6 :

 6 Figure 6: Charge density difference analysis showing the net flow of charge due to the creation of a subsurface vacancy in the slab. Charge flows from regions highlighted in green to regions highlighted in yellow.

Figure 7 :

 7 Figure 7: Diagrams of molecular (a) -(c) and dissociative (d) -(f) adsorption of water molecules on the defective surface for 2/6, 5/6, and 6/6 ML coverages. Oxygen atoms of TiO 2 slab associated with vacancy migration (O A , O B , and O C ) are colored green, while position of an oxygen vacancy is shown in yellow. O w is an oxygen atom of a water molecule. The arrows in Fig. 7(a) show the respective movement of the O atoms during reconstruction and vacancy migration. Cyan: Ti atom; red: O atom; white: H atom

Table 1 : Comparison of molecular and dissociative adsorption energies of a single water molecule on perfect and defective anatase (101) surfaces

 1

	Source	perfect surface	defective surface
		(eV)	(eV)	(eV)	(eV)
	This work	-0.81	-0.42	-0.99	-2.40
	Ref 12	-0.70	-0.32	-0.85	-1.83
	Ref 14	-0.72	-0.34	-1.01	-0.70
	Ref 15	--	--	-0.84	-1.59

Table 2 : Incremental non-dissociative adsorption energies of water molecules on both stoichiometric and defective TiO 2 anatase (101) surfaces

 2

	Coverage (n)	∆E ads,n (eV)
		stoichiometric	defective

Table 3 : Values of average molecular and dissociative adsorption energies on stoichiometric and defective surfaces perfect surface defective surface

 3 shows how the energy required to dissociate each adsorbed water molecule (∆E diss /molecule) changes with the number of adsorbed molecules for both the defective and the stoichiometric surface. The presence of a subsurface V O enhances dissociation of adsorbed water molecules up to 2/3 ML (four water molecules).

	Coverage (n)	E ads (eV)	E diss (eV)	E ads (eV)	E diss (eV)

Acknowledgements

The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu