E Duchêne

A S Fraenkel

S Gravier

R J Nowakowski

Another bridge between Nim and Wythoff

Keywords: Nim, Wythoff's game, Sprague-Grundy function, impartial combinatorial game AMS subject classification: 91A46 Another Bridge between Nim and Wythoff

The P positions of both two-heap Nim and Wythoff's game are easy to describe, more so in the former than in the latter. Calculating the actual G values is easy for Nim but seemingly hard for Wythoff's game. We consider what happens when the rules for removing from both heaps are modfied in various ways.

Introduction

In the game of Nim, played on two heaps of tokens, the two players alternate in choosing a heap and removing any positive number of tokens from that heap. Wythoff's game is also played with two heaps, as in Nim, a player may remove any positive number from a single heap or the same positive number from both heaps, subject to the proviso that every heap remains of nonnegative size at all times. For both games, and in all other games considered in the sequel, the player first unable to move loses.

Definition 1 For a position v of an impartial game, let Op(v) be the set of the options of v. That is Op(v) is the set of all the positions that can be reached from v in one move.

Definition 2 A P position is one in which the next player has no winning move and in an N position, the next player does have a winning move.

In Wythoff and in nim, the position in which both heaps are empty is a P position and any position that has a move to a P position is an N position. Every move from a P position leads to an N position.

Let U ⊂ Z ≥0 , U = Z ≥0 . The Minimum EXcluded value of U , denotd by M ex(U ), is the smallest nonnegative integer not in U . In particular, M ex(∅) = 0.

The P and N classsification of positons can be refined. To each game position v of an impartial game we associate a nonnegative integer value G(v), called the G-value of v. This function G is called the Sprague-Grundy function. It can be defined recursively as follows:

G(v) = M ex({G(u) : u ∈ Op(v)}).
It is well-known that the 0s of the Grundy function constitute the P positions of a game. (See e.g. [START_REF] Berlekamp | Winning Ways for Your Mathematical Plays (in 4 volumes)[END_REF][START_REF] Conway | On Numbers and Games[END_REF][START_REF] Fraenkel | The generalized Sprague-Grundy function and its invariance under certain mappings[END_REF] for more information on the Gfunction.) Note that this function exists uniquely for any finite impartial acyclic game. In an acyclic game, each game position is reached at most once.

Definition 3 Given two positive integers a and b, a mod b denotes the smallest non-negative remainder of the division of a by b. Definition 4 We denote the nim-sum of a and b by a ⊕ b, that is, addition in binary without carries. (Also known as XOR or addition over GF [START_REF] Blass | The Sprague-Grundy function for Wythoff 's game[END_REF] of a and b.)

In two heap Nim, the G-value of (a, b) is a⊕b and is a P position precisely when a = b. In Wythoff's game, the P positions are ( nτ , nτ 2 ), where τ = (1 + √ 5)/2 is the golden ratio. The non-zero G-values appear to be difficult to calculate (cf. [START_REF] Blass | The Sprague-Grundy function for Wythoff 's game[END_REF][START_REF] Nivasch | More on the Sprague Grundy function for Wythoff 's game[END_REF]), however they exhibit additive periodicity. See [START_REF] Dress | Additive periodicity of the Sprague-Grundy function of certain Nim games[END_REF]. A much simplified proof is given in [START_REF] Landman | A simple FSM-based proof of the additive periodicity of the Sprague-Grundy function of Wythoff 's game[END_REF].

In the literature, several variations of Wythoff's game were investigated, some concerning its P positions, others its Sprague-Grundy function. The variations can be subdivided mainly into two categories: (i) extensions, i.e., adjoining new rules to those of Wythoff, and (ii), restrictions, where only certain subsets of Wythoff's moves are permitted. Most investigations concern (i). Examples are [START_REF] Fraenkel | How to beat your Wythoff games' opponent on three fronts[END_REF], where the "diagonal" move (taking from both piles) is relaxed to taking k > 0 from one pile, > 0 from the other, subject to |k -| < a, where a is a fixed positive integer parameter. In [START_REF] Fraenkel | Heap games, numeration systems and sequences[END_REF], this rule is further extended to permit diagonal moves of the form 0 ≤ k -< (s -1)k + t, k ∈ Z ≥1 , where s, t ∈ Z ≥1 are fixed parameters. See also [START_REF] Connell | A generalization of Wythoff 's game[END_REF][START_REF] Fraenkel | Euclid and Wythoff games[END_REF][START_REF] Fraenkel | The game of End-Wythoff[END_REF]. Still other extensions are generalizations to more than 2 piles [START_REF] Fraenkel | The structure of complementary sets of integers: a 3-shift theorem[END_REF], [START_REF] Sun | On Fraenkel's N -heap Wythoff's conjectures[END_REF], [START_REF] Sun | Wythoff 's sequence and N -heap Wythoff 's conjectures[END_REF].

Examples of (ii) are [START_REF] Fraenkel | The game of End-Wythoff[END_REF], where m ≥ 2 piles are considered to be components of a vector, and removals can be made only from the first and the last end piles of the vector. In [START_REF] Fraenkel | Nimhoff games[END_REF], the diagonal moves of Wythoff are restricted in certain ways. In [START_REF] Fraenkel | A generalization of Wythoff 's game[END_REF] the moves from a single pile are restricted, and the diagonal moves are both restricted and extended!

In this paper, we define a new variation of Wythoff's game called W yt K . The rules are more restrictive than in classical Wythoff's game but are in the same spirit as [START_REF] Fraenkel | Nimhoff games[END_REF], which appears to constitute the first bridge between Nim and Wythoff's game. In [START_REF] Fraenkel | Nimhoff games[END_REF], the authors deal with games where the "diagonal move" (i.e., taking k > 0 from one pile and > 0 from the other) is subject to a relation R(k, ). Such games are called Nimhoff games. Wythoff's game is Nimhoff's game where R(k, k) for all k > 0, whereas Nim is the game where no pair (k, ) satisfies R. The main objective there is to a find a closed formula for the G-values of Nimhoff games, for some particular relations R. The family of cyclic Nimhoff games is widely studied in [START_REF] Fraenkel | Nimhoff games[END_REF]. This family contains games of the type R(k, ) if 0 < k +l < h, where h is a fixed positive integer. In addition to these results, the cases R(1, 1) and R(k, k) for all k being a power of 2 are investigated. A generalized Nim sum is provided to ensure the polynomiality of the Gfunction of such games.

In the present paper, the games W yt K that we investigate are exactly the subset of Nimhoff games corresponding to restrictions of Wythoff's game. Unlike the previous paper, we here focus on the regularity of Gfunctions (defined in Sect. 2). The set W yt K is an illustration of games having a certain regular G-function that we call p-Nim regular. For that purpose, the games R(1, 1) and R(2 q , 2 q ) of [START_REF] Fraenkel | Nimhoff games[END_REF] are also considered here as instances of W yt K having such a regular G-function. Besides, we deal with the conditions for which a game does or does not have a p-Nim regular G-function.

Definition 5 Let K be a subset of the positive integers. The game Wyt K is played with two heaps of tokens and

Op(a, b) = {(a -i, b) : 0 < i ≤ a} ∪ {(a, b -j) : 0 < j ≤ b} ∪{(a -k, b -k) : 0 < k ≤ min{a, b}, k ∈ K}.
That is, for a given K, Wyt K is Wythoff's Nim but with a restricted set K of moves along the diagonal. Specifically, we focus on the following questions:

1. What are the P positions for Wyt K ? 2. For any non-negative integer j, is there an a j such that (a j , a j + j) is a P position?

3. What are the G-values and do they exhibit any regularity?

The interest in the first question is clear. The second is an indication of how close the game is to Wythoff's game. The third is clear but needs a little explanation. Subtraction games have periodic Sprague-Grundy functions; many infinite octal games (including one-heap nim) have arithmetic periodic Sprague-Grundy functions. As noted, the rows of the Sprague-Grundy function of Wythoff's game are ultimately additive periodic [START_REF] Dress | Additive periodicity of the Sprague-Grundy function of certain Nim games[END_REF], [START_REF] Landman | A simple FSM-based proof of the additive periodicity of the Sprague-Grundy function of Wythoff 's game[END_REF]. For many other games, when a player is analyzing a new game, hand calculations are usually tried first, varying the value of just one or two heaps, say of size k, and calculating the corresponding Sprague-Grundy function, call it G(k), Even though, initially, this sequence can hold the promise of regularity, the appearance of values k where G(k ) ≥ k is a typical indicator of impending chaos. This is the motivation behind our definition of nim-regularity in section 2, where we give an automatic test (one suitable for computers) for checking for this regularity. This test forms the basis for our positive results, but it also leads, later on, to conditions where games do not have any of the aforementioned periodicities, though it may have other regularities. This negative result does not appear (explicitly) in [START_REF] Fraenkel | Nimhoff games[END_REF][START_REF] Blass | The Sprague-Grundy function for Wythoff 's game[END_REF][START_REF] Blass | How far can Nim in diguise be stretched?[END_REF].

In section 3 we show that when K = {k}, the P positions of Wyt K are nim-regular. For k even, this is (essentially) Lemma 10 of [START_REF] Fraenkel | Nimhoff games[END_REF]. We complete the picture for k odd in this section.

In section 4 we consider the case where K contains only powers of 2. In 4.1 we handle the case K = {1}, followed by stating, in the present language, the case K = {2 k } for fixed k > 0, already given in [START_REF] Fraenkel | Nimhoff games[END_REF]. We then state and prove the negative result alluded to earlier. In section 4.2 we deal with the case where |K| is an infinite set of powers of 2. The case K = {1, 2 i , i ∈ I ⊆ Z ≥1 }, turns out to be equivalent to the game where K = {1}. In that case, we show a surprising regularity of the G-values modulo 3. We wrap up with a brief final section 5.

Closed p-Nim Regularity Check

Our basic definitions are the following.

Definition 6 Let A be a doubly, semi-infinite matrix and A p be the finite matrix consisting of the first p rows and first p columns of A. The matrix

A is called p-nim-regular if A(a, b) = p a p ⊕ b p + A p (a mod p, b mod p) for a, b ∈ Z ≥0 ;
if, in addition, each row and each column of A p contains all the integers 0 through p -1 then A is called closed p-nim-regular. A game whose G-values constitute a (closed) p-nim-regular matrix, is said to be a (closed) p-nim-regular game.

Figure 1 illustrates this definition with p = 4. Roughly speaking, one can say that a p-nim-regular matrix is obtained by tiling the quarter-plane with copies of A p scaled according to the Nim matrix. Notation 1 For a given Wyt K and p ∈ Z ≥1 , denote by A K (p) a p×p matrix where A K (p)(i, j) = G(i, j), i.e. the (i, j) entry of A K (p) is the G-value of the position with heaps of size i and j, for all 0 ≤ i, j < p.

Note that if A is a closed p-nim-regular matrix, then in the matrix A 2 k p , every row and column contains all the values 0 through 2 k p -1; i.e., this matrix is a latin square.

Theorem 1 of [START_REF] Fraenkel | Nimhoff games[END_REF] shows that the following game is p-nim-regular: given 2 heaps and allowing the subtraction of a i from heap i, i ≤ 2 where a 1 +a 2 < p. Here, Theorem 5 notes that not all Wyt K with |K| = 1 are closed p-nimregular.

We present an automatic check for closed p-nim-regularity for any finite set K.

Lemma 1 Let K ⊆ Z ≥1 be a finite set and A be the matrix with entries G(a, b). If there is a positive integer p > max K such that (i) each row and each column of A K (p) contains all the integers 0 through p -1; and

(ii) G(a + p -k, b + p -k) = G(a, b) for 0 ≤ a + p -k, b + p -k < p, 0 ≤ a, b < p and k ∈ K,
then A is closed p-nim regular.

Figure 1 above shows the first few G-values of Wyt K for K = {2}. In that case, the value p = 4 makes A K (p) satisfy both conditions of Lemma 1.

Proof: Assume there exists some p > max K satisfying both conditions (i) and (ii). Now, let (M n ) be the following sequence of matrices:

M n = M n-1 M n-1 + 2 n-1 p M n-1 + 2 n-1 p M n-1 for all n ≥ 1. Set M 0 = A K (p).
We will now prove four properties about the sequence (M n ):

1. M n = A K (2 n p). 2. M n (a, b) = p( a p ⊕ b p ) + A K (p)(a mod p, b mod p) for 0 ≤ a, b < 2 n p.
3. Each row and each column of M n contains all the integers 0 through

2 n p -1. 4. M n (a + 2 n p -k, b + 2 n p -k) = M n (a, b) for 0 ≤ a + 2 n p -k, b + 2 n p -k < 2 n p, 0 ≤ a, b < 2 n p and k ∈ K.
One can check that these properties are true for M 0 = A K (p). Now suppose that they are true for some matrix M n-1 with n ≥ 1, and consider the matrix M n .

1. We will prove that

M n = A K (2 n p). Let 0 ≤ a, b < 2 n-1 p. By the induction hypothesis, we have G(a, b) = M n (a, b) = M n-1 (a, b).
Now consider the position (a, b + 2 n-1 p). According to the rules of the game, we have

G(a, b + 2 n-1 p) = mex({G(a -i, b + 2 n-1 p) : 0 < i ≤ a} ∪{G(a, b + 2 n-1 p -j) : 0 < j ≤ b + 2 n-1 p} ∪{G(a -k, b + 2 n-1 p -k) : 0 < k ≤ a, k ∈ K}).
One may assume inductively that G(s, t) = M n (s, t) for all the pairs (s, t) = (a, b + 2 n-1 p) satisfying 0 ≤ s ≤ a and 0 ≤ t ≤ b + 2 n-1 p. With this hypothesis and by construction of M n , we have

{G(a -i, b + 2 n-1 p) : 0 < i ≤ a} = {G(a -i, b) + 2 n-1 p : 0 < i ≤ a}. Similarly {G(a, b + 2 n-1 p -j) : 0 < j ≤ b + 2 n-1 p} = {G(a, b + 2 n-1 p -j) : 0 < j ≤ b} ∪{G(a, b + 2 n-1 p -j) : b < j ≤ b + 2 n-1 p} = {G(a, b -j) + 2 n-1 p : 0 < j ≤ b} ∪{0, 1, 2, . . . , 2 n-1 p -1} and {G(a -k, b + 2 n-1 p -k) : 0 < k ≤ a, k ∈ K} = {G(a -k, b + 2 n-1 p -k) : 0 < k ≤ a, b, k ∈ K} ∪{G(a -k, b + 2 n-1 p -k) : b < k ≤ a, k ∈ K} = {G(a -k, b -k) + 2 n-1 p : 0 < k ≤ a, b, k ∈ K} ∪{G(a -k, b + 2 n-1 p -k) : b < k ≤ a, k ∈ K}. Since {G(a -k, b + 2 n-1 p -k) : b < k ≤ a, k ∈ K} ⊆ {0, 1, 2, . . . , 2 n-1 p -1}, we have that G(a, b + 2 n-1 p) = mex({G(Op(a, b)) + 2 n-1 p} ∪ {0, 1, 2, . . . , 2 n-1 p -1}) = G(a, b) + 2 n-1 p. Now consider the position (a + 2 n-1 p, b + 2 n-1 p). Then G(a + 2 n-1 p, b + 2 n-1 p) = mex({G(a + 2 n-1 p -i, b + 2 n-1 p) : 0 < i ≤ a + 2 n-1 p} ∪{G(a + 2 n-1 p, b + 2 n-1 p -j) : 0 < j ≤ b + 2 n-1 p} ∪{G(a + 2 n-1 p -k, b + 2 n-1 p -k) : 0 < k ≤ min(a + 2 n-1 p, b + 2 n-1 p), k ∈ K}).
As previously, suppose that G(s, t) = M n (s, t) for all the pairs (s, t) = (a+2 n-1 p, b+2 n-1 p) satisfying 0 ≤ s ≤ a+2 n-1 p and 0 ≤ t ≤ b+2 n-1 p. Then we have:

{G(a + 2 n-1 p -i, b + 2 n-1 p) : 0 < i ≤ a + 2 n-1 p} = {G(a -i, b) : 0 < i ≤ a} ∪ {2 n-1 p, . . . , 2 n p -1} {G(a + 2 n-1 p, b + 2 n-1 p -j) : 0 < j ≤ b + 2 n-1 p} = {G(a, b -j) : 0 < j ≤ b} ∪ {2 n-1 p, . . . , 2 n p -1} {G(a + 2 n-1 p -k, b + 2 n-1 p -k) : 0 < k ≤ min(a + 2 n-1 p, b + 2 n-1 p), k ∈ K} = {G(a -k, b -k) : 0 < k ≤ min(a, b), k ∈ K} ∪{G(a + 2 n-1 p -k, b + 2 n-1 p -k) : k > min(a, b), k ∈ K}. Hence G(a + 2 n-1 p, b + 2 n-1 p) = mex({G(Op(a, b)} ∪{2 n-1 p, . . . , 2 n p -1} ∪{G(a + 2 n-1 p -k, b + 2 n-1 p -k) : k > min(a, b), k ∈ K}). We already know that G(a, b) / ∈ {2 n-1 p, . . . , 2 n p -1}, since G(a, b) ∈ M n-1 . Moreover, G(a, b) / ∈ {G(a+2 n-1 p-k, b+2 n-1 p-k) : k > min(a, b), k ∈ K}. Indeed, if a + 2 n-1 p -k ≥ 2 n-1 p or b + 2 n-1 p -k ≥ 2 n-1 p, then G(a + 2 n-1 p -k, b + 2 n-1 p -k) ≥ 2 n-1 p > G(a, b). Otherwise, if a + 2 n-1 p -k < 2 n-1 p and b + 2 n-1 p -k < 2 n-1 p, it is true since M n-1 satisfies Property (4). Therefore, we conclude that G(a + 2 n-1 p, b + 2 n-1 p) = G(a, b).
2. Let 0 ≤ a, b < 2 n p. We will prove that the formula

M n (a, b) = p( a/p ⊕ b/p ) + A K (p)(a mod p, b mod p) is satisfied. • If 0 ≤ a, b < 2 n-1 p, then it is true by induction. • If 2 n-1 p ≤ a, b < 2 n p, then by construction of M n , we have M n (a, b) = M n (a -2 n-1 p, b -2 n-1 p). ( 2A
)
By the induction hypothesis, we expand (2A):

M n (a, b) = M n (a -2 n-1 p, b -2 n-1 p) = p( a/p -2 n-1 ⊕ b/p -2 n-1 ) + A K (a mod p, b mod p) = p( a/p ⊕ b/p ) + A K (a mod p, b mod p).
• If a ≥ 2 n-1 p and b < 2 n-1 p, then by construction:

M n (a, b) = M n (a -2 n-1 p, b) + 2 n-1 p. (2B)
By the induction hypothesis, we know that

M n (a -2 n-1 p, b) = p( a/p -2 n-1 ⊕ b/p ) + A K (a mod p, b mod p).
Finally, since a ≥ 2 n-1 p and b < 2 n-1 p, and according to the properties of the nim-sum, we get from (2B)

M n (a, b) = 2 n-1 p + p( a/p -2 n-1 ⊕ b/p ) + A K (a mod p, b mod p) = p( a/p ⊕ b/p ) + A K (a mod p, b mod p).
• If b ≥ 2 n-1 p and a < 2 n-1 p, then we reduce to the previous case by symmetry.

3. By construction of M n from M n-1 , and since M n-1 satisfies property (3) by the induction hypothesis, one can easily check that each row and each column of M n contains all the integers 0 through 2 n p -1.

We will show that

M n (a + 2 n p -k, b + 2 n p -k) = M n (a, b) for 0 ≤ a + 2 n p -k, b + 2 n p -k < 2 n p, 0 ≤ a, b < 2 n p and k ∈ K. The condition 0 ≤ a + 2 n p -k, b + 2 n p -k < 2 n p implies a, b < k.
And since k < p, we now consider that 0 ≤ a, b < p.

We now define two integers A and B such that:

A = a + (2 n -1)p B = b + (2 n -1)p.
By the formula proved in (2), we have

M n (a, b) = M n (A, B).
Let (X, Y ) be a position defined as follows:

X = a + 2 n p -k Y = b + 2 n p -k. Thus M n (X, Y ) = M n (A + p -k, B + p -k) = M n (a + p -k, b + p -k)
according to the formula. It now suffices to prove that

M n (a + p -k, b + p -k) = M n (a, b). Since a + 2 n p -k, b + 2 n p -k < 2 n p, we have a -k + p, b -k + p < p,

P positions

Theorem 2 Let K = {k}. Then in Wyt K , 1. For k = 2j, the P positions are

(i + 2kp, i + 2kp), i = 0, 1, . . . , k -1, p ∈ Z ≥0 and (k +2i+2kp, k +2i+2kp+1), i = 0, 1, . . . , j -1, p ∈ Z ≥0 .
2. For k = 2j + 1, the P positions are

(i + (2k + 1)p, i + (2k + 1)p), i = 0, 1, . . . , k-1, p ∈ Z ≥0 and (k+2i+(2k+1)p, k+2i+(2k+1)p+1), i = 0, 1, . . . , j, p ∈ Z ≥0 .
Proof: Denote by S ⊂ Z 2 the set of positions described by the theorem. Denote by A the set Z 2 \ S. We must prove that any move from S lands in a position of A, and that from any position of A, there exists a move leading to a position of S.

For any subset T of Z 2 , denote by T | x the set {(i, j) ∈ T : 0 ≤ i, j < x}. Denote by B the value 2k (resp. 2k + 1) if k is even (resp. odd). We first remark that there is exactly one position of S in each row and in each column. It is then straightforward to see that it is not possible to move from a position of S| B to another one. Besides, from each position of A| B , one can move to a position of S| B . Therefore, from any position of the sets {(i, j) : i ≥ B, j < B} and {(i, j) : j ≥ B, i < B}, one can reach a position of S| B . Now consider a position of S ∩ {(i, j) : B ≤ i, j < 2B}. Moves in a single heap clearly land in A. A move of length k in both heaps may land in Z 2 | B , but not in S| B (see figure 2). Also, from any position of

A ∩ {(i, j) : B ≤ i, B ≤ j < 2B ∨ B ≤ j, B ≤ i < 2B}, one can move to a position of S ∩ {(i, j) : B ≤ i, j < 2B}.
It now suffices to iterate the result for the sets S ∩ {(i, j) : pB ≤ i, j < (p + 1)B} and A ∩ {(i, j) : (pB ≤ i, pB ≤ j < (p + 1)B) or (pB ≤ j, pB ≤ i < (p + 1)B)}, with p > 1.

Remark 1 When K = {k}, this result ensures that the only integers j for which there exists an a j such that (a j , a j + j) is a P position are 0 and 1.

This remark leads us to the following conjecture:

Conjecture 1 Let K be a finite set, then there exists an integer

J K > 0 such that if (a, b) is a P position for Wyt K then |a -b| < J K . 4 G-Values 4.1 K is finite Theorem 3 Let K = {1}. Then in Wyt K G(3m + i, 3n + j) = 3(m ⊕ n) + A K (3)(i, j), 0 ≤ i, j < 3, ∀m, n ∈ Z ≥0
where,

A K (3) =   0 1 2 1 2 0 2 0 1   .
Figure 3 illustrates this result by depicting the table of the first G-values. Proof: One can check that A K (3) contains the first G-values of the game Wyt K with K = {1}. With p = 3 > max K, it is now straightforward to see that the conditions of Lemma 1 are satisfied.

For completeness, we report the cases where K = {2 j } since the language used in [START_REF] Fraenkel | Nimhoff games[END_REF] does not immediately lend itself to this interpretation.

Theorem 4 Let K = {2 k } for fixed k > 0. Then in Wyt K G(2 k+1 m + i, 2 k+1 n + j) = 2 k+1 (m ⊕ n) + A K (2 k+1 )(i, j), 0 ≤ i, j < 2 k+1 , ∀m, n ∈ Z ≥0
and, as a 2 × 2 array of matrices,

A K (2 k+1 ) = A K (2 k ) A K (2 k ) ⊕ 2 k A K (2 k ) ⊕ 2 k A K (2 k ) ⊕ 1 , where A K (2 k ) is the 2 k × 2 k matrix of the G-values of Nim with two heaps of sizes lower than 2 k .
As an illustration of that case with k = 1, one can refer to Figure 1.

Theorem 5 Let K = {k : k = 2 q = k + 1 ∀q ∈ Z ≥1 } (for any integer q > 0). Then there is no p such that A K (p) satisfies Lemma 1, condition (i), i.e. Wyt K is not closed p-nim-regular for any p.

Proof: We suppose that such a p exists. First, let K = {2j}, j = 2 q . From Theorem 2, the P-positions repeat with period 4j, consequently, then p would be also be a multiple of 4j. There exists an i such that 2 i < 2j and with the property that 2 i ⊕ 2j = 2 i + 2j. Moreover, since 2j is not a power of 2 then 2 i+1 < 2j also holds and there is no diagonal move available from (2 i+1 , 4j -1). Hence the G-value of the position is 2 i+1 ⊕ 4j -1-this is 2-heap nim. Then we have 2 i+1 ⊕ 4j -1 = 2 i+1 + 4j -1 ≥ 4j. Therefore, since A K (4j) contains a number greater than 4j -1 not every row and column can contain all the numbers 0 through 4j -1.

For K = {2j + 1}, the argument is similar. ¿From Theorem 2, p would be a multiple of 4j + 3. Since 2j + 2 is not a power of 2 then there exists 2 i ∈ 2j + 1, 1 < 2 i < 2j + 1 so that 2 i ⊕ 2j + 1 = 2 i + 2j + 1. Therefore, 13 2 i+1 ⊕ 4j + 2 = 2 i+1 + 4j + 2 > 4j + 3. The position (2 i+1 , 4j + 2) has a G-value equal to 2 i+1 ⊕ 4j + 2. Therefore, since A K (4j + 3) contains a number greater than 4j + 3 not every row and column can contain all the numbers 0 through 4j + 2.

|K| is an infinite set of powers of 2

When K is a subset of the powers of two including 1, we show that the G-values of Wyt K are those described by Theorem 3 (see Figure 3).

Theorem 6 Let K = {1, 2 i : i ∈ I ⊆ Z ≥1 }. Then in Wyt K , G(3m + i, 3n + j) = 3(m ⊕ n) + A K (3)(i, j), 0 ≤ i, j < 3, ∀m, n ∈ Z ≥0 , where A K (3) is as defined in Theorem 3.
Proof: Denote by G 1 (a, b) the G-value of the position (a, b) for Wyt K with K = {1}. The function G 1 is described in Theorem 3. We aim at proving that G 1 and G for Wyt K where K = {1, 2 i : i ∈ I ⊆ Z ≥1 } are identical.

Since the moves of Wyt According to Theorem 3, we can assert that G 1 (x, y) mod 3 = A K (3)(x mod 3, y mod 3) for any position (x, y). This implies that A K (3)(a mod 3, b mod 3) = A K (3)((a -k) mod 3, (b -k) mod 3). When looking at the matrix A K (3), we notice that each value 0, 1, 2 appears exactly once in each diagonal (fig. 4).

K with K = {1} are included in those of Wyt K with K = {1, 2 i : i ∈ I ⊆ Z ≥1 }, it suffices to check that G 1 (a, b) = G 1 (a -k, b -k) for all k in {2 i : i ∈ I ⊆ Z ≥1 }.
Therefore, A K (3)(a mod 3, b mod 3) = A K (3)((a-k) mod 3, (b-k) mod 3) if and only if a mod 3 = (a-k) mod 3 and b mod 3 = (b-k) mod 3. These equalities now imply k mod 3 = 0, which is impossible since k is a power of 2.

Remark 2 From Theorem 6, one can see that in Wyt K with K = {2 i : i ≥ 0}, we have G(3m+i, 3n+j) = 3(m⊕n)+A K (3)(i, j), 0 ≤ i, j < 3, ∀m, n ∈ Z ≥0 .

Conjecture 2 Let K = {2 i : i > 0}. Then for all non-negative integers j, there is an a j such that (a j , a j + j) is a P position. Other games that are both an extension and a restriction of Wythoff are suggested by the present paper. For example, let a ∈ Z ≥1 , K ⊂ Z ≥1 . The diagonal move is extended as follows: take k > 0 from one pile and > 0 from the other subject to |k -| < a (extension -see [START_REF] Fraenkel | How to beat your Wythoff games' opponent on three fronts[END_REF]) and k ∈ K (restriction). The extension [START_REF] Fraenkel | Heap games, numeration systems and sequences[END_REF] can be restricted similarly. (Although in [START_REF] Fraenkel | The game of End-Wythoff[END_REF] there are m ≥ 2 piles, this is not a genuine extension, since all moves are restricted to taking from at most 2 piles.)

Figure 1 :

 1 Figure 1: The first G-values of Wyt K when K = {2}: an example of a closed 4-nim-regular matrix

  which means that the positions (a, b) and (a + p -k, b + p -k) both belong to the square [0, . . . , p -1] × [0, . . . , p -1]. As A K (p) satisifes condition (ii), we deduce that M n (a + p -k, b + p -k) = M n (a, b).

Figure 2 Figure 2 :

 22 Figure 2: A view of the structure S.

12

 12

Figure 3 :

 3 Figure 3: The first G-values of Wyt K for K = {1}.

  Suppose that there exists a position (a, b) and an integer k= 2 i , i ≥ 1 such that G 1 (a, b) = G 1 (a -k, b -k). Then we also have G 1 (a, b) mod 3 = G 1 (a -k, b -k) mod 3.

Figure 4 :

 4 Figure 4: The three diagonals of A K (3) modulo 3


Hugo Girard

Behrad Koohbor
email: koohbor@rowan.edu

Aurélien Doitrand

Robert Livingston

Experimental Characterization of In-Plane Debonding at Fiber-Matrix Interface Using Single Glass Macro Fiber Samples

Keywords: A. Glass fibres, B. Debonding, C. Damage mechanics, D. Mechanical Testing

Fiber-matrix interface debonding is investigated using a hybrid experimental-computational approach. Model composite specimens are prepared using a single glass macro fiber embedded in an epoxy test piece and then subjected to uniaxial tension. Displacement and strain fields developed in the vicinity of the fiber-matrix interface are measured experimentally using optical digital image correlation. Debonding initiation and propagation stages are characterized using the transient strain evolutions at the fiber-matrix interface.

Correlations between local and global stresses are identified by combining full-field measurement data and analytical solutions. Local deformation data obtained experimentally is then used to identify a stress-based traction-separation model. Experimental measurements are verified and supplemented by finite element analyses. Sources of uncertainty and the limitations associated with the proposed hybrid approach are discussed. The approach presented herein provides a systematic methodology for the quantification of the mechanics and failure at the fiber-matrix interfaces in unidirectional composites under transverse loading conditions.

Introduction

Fiber-matrix interfacial debonding is a type of damage that occurs at relatively low stresses in fiber composites subjected to transverse tension [1]. Upon the application of an increasing tensile load, the debonded zone propagates circumferentially along the fiber-matrix interface, eventually kinking out toward the matrix [2]. The coalescence of the kinked-out microcracks among neighboring fibers leads to the development of larger cracks that grow rapidly across the entire thickness of the lamina, eventually causing through-thickness cracks in transverse plies. After reaching interlaminar regions, these transverse cracks can lead to delamination and/or fiber breakage [3]. The fundamental mechanisms governing the mechanical and failure behaviors at the fiber-matrix interface can be characterized in a number of ways. Most common experimental techniques developed to measure the interface strength involve pulling (or pushing) the fiber out of the matrix while measuring the force required to overcome the interface shear strength. Among these techniques are the well-established pull/push-out [4], microbond (or microdroplet) [START_REF] Tamrakar | Rate dependent mode II traction separation laws for S-2 glass/epoxy interface using a microdroplet test method[END_REF], fragmentation [START_REF] Kim | Observation of fiber fracture and interfacial debonding phenomena using the fragmentation test in single fiber composites[END_REF], and Broutman tests [START_REF] Ageorges | Single-fibre Broutman test: Fibre-matrix interface transverse debonding[END_REF]. While these well-established techniques have been utilized for decades and are capable of providing useful information on the interfacial shear strength (i.e., where debonding growth is parallel to the fiber axis), they fail to provide realistic measures for the tensile strength at the interfaces; needless to say, the latter is the dominant failure mode developed at the fiber-matrix interface in transverse loading conditions.

Investigations reporting the multiscale mechanics and failure behaviors at the fiber-matrix interface and the transverse cracking in composites overwhelmingly use modeling approaches [START_REF] Tan | The role of interfacial properties on the intralaminar and interlaminar damage behaviour of unidirectional composite laminates: Experimental characterization and multiscale modelling[END_REF][START_REF] Velasco | BEM multiscale modelling involving micromechanical damage in fibrous composites[END_REF][START_REF] Velasco | Interaction between fibers in the transverse damage in composites[END_REF][START_REF] Zhang | IGFEM-based shape sensitivity analysis of the transverse failure of a composite laminate[END_REF][START_REF] Shakiba | Detecting transverse cracks initiation in composite laminates via statistical analysis of sensitivity data[END_REF][START_REF] Sepasdar | Micromechanical study of multiple transverse cracking in cross-ply fiberreinforced composite laminates[END_REF]. However, the development of high magnification full-field techniques, especially digital image correlation (DIC), has provided new opportunities for accurate experimental measurement of material behavior at dissimilar interfaces and over a broad spectrum of length scales. A perfect example of such measurements made possible by high magnification optical and scanning electron microscopy (SEM) DIC is the characterization of the debonding behavior at the fiber-matrix interface in transversely-loaded composites. DIC has been applied to a variety of composite systems, including ceramic and polymer-matrix composites, to identify the underlying mechanisms and quantify the local deformation behaviors at the fiber-matrix interface [START_REF] Tracy | Multiscale damage characterization in continuous fiber ceramic matrix composites using digital image correlation[END_REF][START_REF] Mehdikhani | Full-field strain measurements at the micro-scale in fiber-reinforced composites using digital image correlation[END_REF][START_REF] Montgomery | A robust patterning technique for electron microscopybased digital image correlation at sub-micron resolution[END_REF][START_REF] Koohbor | Identification of RVE length scale in fiber composites via combined optical and SEM digital image correlation[END_REF]. Recent studies have proved the unique advantages of SEM-DIC techniques in the direct measurement of local deformation and strain fields developed in the vicinity of glass and carbon fibers in polymer matrix composites [START_REF] Montgomery | A robust patterning technique for electron microscopybased digital image correlation at sub-micron resolution[END_REF][START_REF] Koohbor | Identification of RVE length scale in fiber composites via combined optical and SEM digital image correlation[END_REF]. In addition to SEM-DIC, optical DIC approaches have also been utilized to study the fiber-matrix mechanics in various composite systems. For example, it was shown by Tabiai et al. [START_REF] Tabiai | In-situ full field measurements during inter-facial debonding in single fiber composites under transverse load[END_REF] that high magnification optical DIC measurements carried out on single PTFE macro fiber (fiber diameter ~1 mm) model composite samples can provide useful information about the displacement and strain fields in the fiber vicinity. This data was then used to identify the three stages of the interface damage development, i.e., an initial symmetric debonding formed at diagonally opposite locations at the interface, followed by crack growth along the interface, and finally, matrix failure caused by a crack kinking process [START_REF] Tabiai | In-situ full field measurements during inter-facial debonding in single fiber composites under transverse load[END_REF]. This study was later expanded to include several different fiber-matrix systems. It was documented that as long as the image correlation parameters are selected such that the strong strain gradients at the interface are captured, the fiber size effects are insignificant, at least for the study of fundamental mechanics and failure characteristics at the interface.

Moreover, independent investigations performed on macro fibers with diameters ranging from a few hundred micrometers to 1 mm were compared, validated, and shown to provide consistent results as long as the fiber edge and the local strain extraction points are far enough from one another (e.g., one subset distance between the fiber edge and the local strain extraction point) [START_REF] Tabiai | In-situ full field measurements during inter-facial debonding in single fiber composites under transverse load[END_REF][START_REF] Tabiai | In-situ full field out of plane displacement and strain measurements at the micro-scale in single reinforcement composites under transverse load[END_REF]. A combined DIC-image processing technique was then proposed to facilitate autonomous and rapid detection of a debonded interface and subsequent matrix cracking phenomena [START_REF] Tabiai | Hybrid image processing approach for autonomous crack area detection and tracking using local digital image correlation results applied to single-fiber interfacial debonding[END_REF], further highlighting the significance of image-based analyses in the multiscale characterization of fiber-matrix debonding. Besides pure experimental or modeling studies, reports that combine the two approaches also exist. For instance, Naylor et al. [START_REF] Naylor | Mechanically regularized FE DIC for heterogeneous materials[END_REF] extended a regularized finite element (FE) DIC approach to detect the damage initiation from SEM observations of heterogeneous carbon fiberepoxy matrix composite by using a microstructure-based consistent mesh and considering dissimilar material properties in the mechanical regularization. These hybrid approaches can potentially reduce the uncertainty in the determination of debonding angle by FE discretization that avoid elements located on the fiber-matrix interface. Furthermore, the hybrid approach proposed by Naylor et al. [START_REF] Naylor | Mechanically regularized FE DIC for heterogeneous materials[END_REF] was capable of measuring interfiber displacement fields with errors as small as a few nanometers.

The potential applications of image-based characterizations in the analysis of fiber-matrix interface mechanics have become evident. However, certain gaps still exist in the practical and quantitative protocols that can be useful in correlating global and local stress and deformation fields in this problem. The present work addresses one of these gaps by providing a methodology that enables a quantitative assessment of the fiber-matrix interface debonding process. To this end, an experimental protocol is designed to enable the measurement of local deformation and strain fields in the vicinity of a single glass macro fiber in an epoxy matrix. Local strain data measured by high magnification optical DIC are then used in conjunction with the bi-material interfacial crack theory to establish a correlation between local and global stress histories. Finally, the fiber/matrix interface strength under normal loading conditions is evaluated by proposing simple tractionseparation laws. The experimental measurements herein are accompanied by FE analyses to confirm the accuracy and reliability of the measurements while highlighting the potential sources of uncertainty associated with experimental characterizations. A unique contribution of this study compared with the previously published reports (e.g., Tabiai et al. [START_REF] Tabiai | In-situ full field measurements during inter-facial debonding in single fiber composites under transverse load[END_REF][START_REF] Tabiai | In-situ full field out of plane displacement and strain measurements at the micro-scale in single reinforcement composites under transverse load[END_REF]) is the establishment of a theoretically supported and modelingsupplemented experimental approach that enables the quantification of critical stages of the fiber-matrix debonding process. The detailed quantitative investigations facilitate the determination of traction-separation laws that describe the mechanics of the debonding process.

Materials and Methods

Modeling-Informed Specimen Size Determination

The first step in the design of experiments is to determine the dimensions and geometry of the specimens. The most significant consideration in the preparation of single fiber specimens is to ensure that the deformation fields developed in the vicinity of the macro fiber remain unaffected by the free boundaries (edges) of the tensile specimens. Ideally, the smallest possible fiber diameter is favored, as it ensures the absence of any possible boundary effect. However, for practical image correlation purposes, one must also consider the lower limits of the imaging system, as it directly affects the resolution and accuracy of the fullfield optical DIC measurements, as discussed in detail in Sec. 2.2.

The specimen size determination in this work was based on a FE analysis performed to characterize the sensitivity of the results to the specimen width for a fixed glass macro fiber diameter. The optimal tensile specimen width is considered to be one in which the deformation fields around the embedded fiber remain the least affected by the free boundaries of the test piece. However, for practical and quantitative image-based characterizations, the measurement resolution must also be carefully taken into consideration. In the present study, the optimal size determination was performed by FE modeling of virtual tensile specimens with widths ranging from 4 to 32 mm, incorporating a single fiber at their center. Due to symmetry, only a quarter of the entire geometry was modeled (see Figure 1a). The ratio between the tensile specimen width and the macro fiber radius was varied at a constant global tensile stress. Variations of the local longitudinal stress developed at various locations and along various paths were considered the criteria for selecting the proper specimen width to fiber radius ratios. Figure 1b-d show the variation of local longitudinal stress at various locations and paths across and/or near the fiber-matrix interface obtained for different w/r ratios. The data in this figure (specifically those shown in Figure 1b) suggest that the boundary effect is more significant at w/r<5, while convergence is observed at larger w/r ratios. For example, considering the local stress values plotted with respect to (w/r) ratio in Figure 1b, it is observed that the local stress varies insignificantly (i.e., differences <5%) at w/r ratios greater than 5. Furthermore, since both initiation and propagation stages are assessed in the following, it is also necessary to consider not only the stress behavior but also energy metrics. Similar sensitivity analyses, i.e., the variation of the energy release rate (ERR) as a function of w/r, were obtained, highlighting the negligible effect of boundary conditions for w/r>5 (see Supplementary Information for more details).

Based on the results obtained herein and the availability of macro fibers, 2 mm glass macro fibers and the dog bone specimen width of 13 mm were selected. The utilized fiber diameter and tensile specimen width lead to w/r=5.5 which is large enough to ensure a negligible boundary effect. Note that the 13 mm specimen width and other dimensions were also consistent with the ASTM D638 standard, ensuring that the global stress state is uniaxial and independent of the gauge length.

Preparation of Single-Fiber Specimens

Tensile specimens embedding a single glass macro fiber were prepared by incorporating a 2 mm diameter quartz glass macro fiber (Momentive Performance Materials, OH, USA) in a clear epoxy matrix.

The two-part epoxy (ProMarine Supplies Co., MI, USA) was mixed and degassed first and then poured inside the 3D printed dog bone mold. The glass macro fiber was then inserted inside the uncured resin and held in place with the aid of a 3D printed jig until the epoxy matrix was fully cured. The utilized epoxy material was cured at room temperature. The room temperature curing effectively mitigates the possibility of thermal residual stress generation at the fiber-matrix interface. grit sandpapers to surface imperfections and ensure a level surface flush with the cross-section of the glass macro fiber. The polished surface was then coated with a thin layer of matte white paint. A dispersed speckle pattern of fine black paint particles was then applied on the white substrate to enable high magnification digital image correlation (DIC) measurements. The coating of the black speckle pattern was achieved by suspending the specimen horizontally with the white-painted surface facing downward. The specimen was then speckled by spraying the black paint from a perpendicular direction to the white face for two to three passes, then rotating the specimen 180 degrees with the white face still facing down. This process was repeated until the deposition of an even speckle pattern was confirmed by visual inspection. This method produced an average speckle particle of ca. 30 µm. The painted surface was allowed to dry for several days before the onset of experimental characterizations. Complete drying of the paint is essential to ensure its adherence to the substrate and that the paint has enough brittleness to crack instantaneously upon the initiation of fiber-matrix debonding. The specimen geometry including the location of the single fiber, the specklepatterned surface, and its grayscale characteristics are shown in Figure 2a-c. Dog bone specimens were subjected to tensile loads applied uniaxially at a constant rate of 1 mm/min until complete failure. A 10 kN load cell with a data acquisition rate of 1 Hz was used. A 2D DIC camera system (5 megapixels, 9.33 µm/pixel) was utilized to capture images from the patterned surface of the specimen during mechanical tests. The imaging rate was synchronized with the load data acquisition. Images captured from the deforming specimen were analyzed using the commercial DIC software Vic-2D (Correlated Solutions, Inc., SC, USA) with the subset and step sizes of 29 pixels (272 µm) and 7 pixels (65 µm), using the local DIC method. Full-field strain maps were calculated using a Gaussian filtering approach with a filter size of 5 data points. The utilized DIC parameters were selected after a sensitivity analysis that identified optimized parameters for a minimal noise-to-strain ratio and ~100% correlation in the images [START_REF] Montgomery | A robust patterning technique for electron microscopybased digital image correlation at sub-micron resolution[END_REF][START_REF] Koohbor | Identification of RVE length scale in fiber composites via combined optical and SEM digital image correlation[END_REF][START_REF] Koohbor | Experimental determination of Representative Volume Element (RVE) size in woven composites[END_REF].

Following the procedure detailed in Koohbor et al. [START_REF] Koohbor | Experimental determination of Representative Volume Element (RVE) size in woven composites[END_REF], the strain measurement bias and noise floor associated with the above DIC parameters were determined as 7.26×10 -5 and 6×10 -4 , respectively.

All the above steps were also applied to neat epoxy specimens (i.e., dog bone specimens without an embedded glass macro fiber) to characterize their mechanical response. Figure 2d shows the stress-strain response of the neat epoxy matrix. Table 1 lists the mechanical properties of the epoxy and glass macro fibers used in this work. Further details about the mechanical characterization of the constituents can be found in Livingston et al. [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by fullfield measurements[END_REF].

Results and Discussion

Full-Field Strain Maps

Figure 3 shows the development of in-plane strain components at different global stress increments.

The longitudinal (εyy) strain maps show an initial development of narrow and localized strain zones formed symmetrically at the top and bottom poles of the glass macro fiber at global stresses as low as 7 MPa. The high strain zones grow into strongly heterogenous strain fields that propagate along the fiber-matrix interface and also through the matrix. The local strains evolve into highly localized strain bands prior to the complete failure of the specimens at ~41 MPa. The formation of these large strain bands is also evident in the other two components of strain, i.e., εxx and εxy. A noteworthy observation is the numerical values of the local strains that far exceed the failure strain of the matrix (see Table 1). As discussed in detail in a previous study, the seemingly high local strains can be associated with the generation of matrix microcracks that remain invisible to the optical measurement system [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by fullfield measurements[END_REF]. The local strain fields associated with such microcracks give rise to values that appear to be higher than the failure strain of the matrix. Evidence of matrix microcracks was presented in a previous study by Livingston and Koohbor [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by fullfield measurements[END_REF]. Another possible source of discrepancy can be due to the fact that the matrix failure strains reported in Table 1 are obtained from the application of a uniaxial one-dimensional stress state (i.e., a uniaxial tensile test). In the vicinity of a bi-material interface with strong property gradients, the stress state can be significantly different from those present in uniaxial stress conditions. As such, more complex plasticity criteria must be used to correlate the local strains with the local stress states. While beyond the scope of this work, we speculate that the development of complex stress states might have been a significant contributing factor to the large local strains developed in the matrix prior to the final failure.

Correlating Local and Global Strains

The full-field nature of the strain measurements in this work allows for the evaluation of local strain values at certain locations of interest. Recent studies suggest that the evolution of local longitudinal strain fields at certain regions around the fiber/matrix interface can be used as an effective metric for characterizing the interfacial debonding onset [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by fullfield measurements[END_REF]. These critical regions of interest include those that are located in diagonally opposite spots near the top and bottom poles of the interface and parallel to the direction of the farfield tensile load. As shown schematically in Figure 4, the matrix regions located above (and below) the interface experience partial elongation during loading and before the initiation of interface debonding. Upon the initiation of the fiber/matrix debonding process, the strained region will experience partial unloading due to the creation of separated traction-free surfaces [START_REF] Mallon | Fracture behavior of prestressed composites subjected to shock loading: A DIC-based study[END_REF]. By further increase of the tensile load, the local strain in the region of interest rises again until complete failure occurs in the specimen.

The above mechanism is investigated in this work by extracting the local strain values from two representative regions of interest located 250 µm away from the fiber-matrix interface, as shown in Figure 5.

Note that the 250 µm offset from the interface was selected to ensure that the local strain analyses remain unaffected by the strain gradients formed at the interface vicinity (see Figure 3). Figure 5 shows the evolution of local longitudinal strains (extracted from top 'T' and bottom 'B' regions of interest) with respect to global stress. For reference, the variation of global (macroscale) strain is also included in the figure . The global strain was determined by a significantly larger linear optical extensometer spanning vertically over the entire area of interest. A careful examination of the local curves (shown in the middle panel) reveals that the local strains increase at a high rate until reaching a peak. The local strain peaks, marked by a vertical dashed line, correspond to the initiation of the debonding at the fiber/matrix interface. The far-field stress value associated with these local strain peaks can be pinpointed by finding the extrema of the local strain curves, as shown in the top and middle panels in Figure 5. The far-field (global) tensile stress associated with the onset of the debonding process in this work is determined as σ ∞ =5.35 MPa. This stress value can be converted to the local normal stress acting on the interface using Goodier's analytical solution and appropriate bi-material constants, k=1.43 and m=1.67, under the plane-strain hypothesis [START_REF] Goodier | Concentration of stress around spherical and cylindrical inclusions and flaws[END_REF][START_REF] García | Debonding at the fibre-matrix interface under remote transverse tension. One debond or two symmetric debonds?[END_REF]. The conversion between local and global stress is given by, ( ) sin ( )

k m σ θ σ θ ∞   = -   , (1)
where, ( ) σ θ is the local radial stress as a function of the angular location along the interface and σ ∞ is the global stress. k and m denote the dimensionless bi-material constants discussed above. Eq. 1 is obtained under the assumption of small deformations and linear elasticity. The matrix material examined in this work exhibits moderate nonlinearity at strains larger than 2%. As a consequence, the linear elastic assumption seems reasonable when focusing on debonding initiation that occurs at strain levels around 0.2% [START_REF] Doitrand | Nonlinear implementation of Finite Fracture Mechanics: A case study on notched Brazilian disk samples[END_REF]. Considering

Eq. 1 and setting θ=90° (i.e., the top pole at the fiber-matrix interface), the local normal stress at the onset of debonding initiation can be determined as 7.65 MPa. Note that due to the relatively weak interface strength in the examined system (due to unsized fibers) and the nearly invisible debonded area, visual verification of the onset of debonding is cumbersome. However, the above method allows for a systematic overestimation of the normal in-plane debonding initiation stress without the need for visual and potentially subjective methods.

Nevertheless, the reliability of this method was further investigated by comparison with visual inspections, FE predictions, and a correlation residual method, discussed in the forthcoming sections.

Additional FE analyses were performed to evaluate the validity of the experimental measurements discussed above. A 2D plane-strain model under linear elasticity and small deformation assumptions was created (see Supplementary Information for details). The local longitudinal strain was extracted from a similar region of interest as in experiments for undamaged and damaged interfaces with different debonding angles. Figure 6a shows the variation of local strain derivative as a function of the interface debonding angle.

The local strain derivative decreases with increasing debonding angle, a behavior consistent with the experimental observations (Figure 5). However, the presence of experimental noise does not allow for a direct comparison between the model-predicted and measured strain derivatives. Instead, the comparison between simulations and experimental measurements can be made by tracking the variation of local axial strain as a function of the applied stress. Figure 6b shows this comparison for a variety of cases, from undamaged (fully bonded interface) to damaged interface with a 105° debonding angle. For comparison, the experimentally measured strains (shown earlier in Figure 5) are also included in this figure. Note that the model-predicted values in this figure are extracted from the same area of interest used for experimental analyses. Before the initiation of interface debonding, there is a good agreement between the local strain variations obtained numerically and experimentally. The larger the debonding angle, the lower the slope of the local strain variation. A comparison between experimental and numerical variations of the local strain slope after the debonding initiation indicates that the debonding angle just after initiation lies somewhere between ~55° and 90°. The latter observation is in agreement with previous reports [START_REF] García | Debonding at the fibre-matrix interface under remote transverse tension. One debond or two symmetric debonds?[END_REF][START_REF] Sandino | Numerical analysis of the influence of a nearby fibre on the interface crack growth in composites under transverse tensile load[END_REF] and is further elaborated in the forthcoming sections.

Characterization of Debonding Opening Displacement

The next step after the determination of debonding initiation stress is the characterization of opening displacements. The normal in-plane opening displacement was characterized through the use of an optical extensometer placed vertically at the top (and bottom) regions of the interface spanning across the interface (see Figure 7a). The vertical gauge displacement caused by the concurrent debonding opening and the stretching of the two materials can be extracted from this optical extensometer as the loading progresses.

Figure 7b shows the evolution of the total vertical gauge displacement, Δ t , as a function of global stress extracted from the top and bottom poles of the fiber-matrix interface. The initial length of the optical extensometers was ~440 µm, which was about 35% larger than the strain filter length (step size × strain filter size) in the DIC analyses. The total vertical displacement shows a slow rate of growth at stresses below the nominal debonding stress (i.e., 5.35 MPa). Upon the initiation of interface debonding, the rate of increase of the vertical gauge displacement increases rapidly, growing at a near steady pace until global stresses of ca. 40 MPa, i.e., the ultimate strength of the single fiber specimen. The stress-strain curves for single-fiber specimens are presented as Supplementary Information.

The evolution of vertical gauge displacement herein shows two notable characteristics. First, the insignificant variability between the data obtained from the top and bottom poles (characterized by the scatter bars) indicates a highly symmetric response, at least over a stress window with upper limits far exceeding the debonding initiation stress. Note that the aforementioned variability grows non-monotonically as the global stress reaches beyond 30 MPa and approaches the ultimate tensile strength of the specimen. The abrupt increase of variation in the vertical displacement is likely due to the unstable growth of the debonded area and/or the development of matrix cracks that impact the local deformation fields in the fiber-matrix interface vicinity. The other noteworthy remark in the interpretation of the vertical displacement data (Figure 7b) is the initial increase in the displacement at stresses below the debonding initiation, i.e., σ ∞ <5.35 MPa. As stated earlier, the displacement values measured across the fiber-matrix interface include the debonding opening displacement due to debonding as well as the contribution of the local strains developed in the parts of the fiber and the matrix that lie within the optical extensometer. However, once debonding initiates, strain relaxation around the debonded area occurs, mitigating the contributions of the fiber and matrix strains at higher global stresses. As a consequence, the experimentally measured Δ t provides an estimate of the in-plane vertical displacement at the uppermost location of the interface.

A fundamental requirement for the interface debonding opening displacement is that the measurements must be independent of the initial gauge length of the optical extensometer. The size independence of the measurement approach described here was verified by FE simulations. The numerical values of the total vertical gauge displacement, Δ t , were obtained from the FE simulations for various gauge length sizes, 𝐿𝐿 0 t . Figure 8 shows the variation of vertical gauge displacement, with respect to debonding angle for three different gauge lengths, 320, 440, and 660 µm. Although the debonding opening displacements show a large variation at small angles, for debonding angles larger than 40°, Δ t provides an accurate and converging estimate (with relative differences smaller than 3%) for the debonding opening displacement. Differences obtained for smaller angles are due to the interface debonding angle being too small such that Δ t is representative of the displacement due to the fiber and matrix local stretches rather than the actual opening of the debonded interface. Nevertheless, experimental evaluation of the interface displacement provides sufficiently accurate measures for the interface opening due to debonding after the initiation phase (i.e., at debonding angles greater than 55°, as discussed in Sec. 3.2).

Stress-Based Normal Traction-Separation Laws

The practical significance of the results discussed in the previous sections is that the output data can be used to identify simple stress-based traction-separation laws that describe the normal debonding behavior in a composite system. Such constitutive equations can be extremely useful in the prediction of transverse cracking behavior and fiber-to-fiber interactions in actual composite components. Before embarking on further details, it should be clarified that the identification of traction-separation models by a pure stress-based approach may lack essential information regarding the energy release rate as another factor influencing the modeling accuracy. Nonetheless, the direct identification of traction-separation models from experimental data and by considering both stress and energy metrics is not practically possible. Instead, there exist inverse techniques that facilitate the identification of high-accuracy traction-separation laws from experimental data by considering the mixed-mode nature of the debonding along the curved fiber-matrix interface. While the authors' current effort includes such inverse identification approaches, the present study is only focused on a 'stress-based' characterization of normal interface debonding behavior for brevity and consistency with the scope.

The evolution of local stress and interface opening displacement can be used as input to construct a simple bi-linear cohesive traction-separation law that describes the normal in-plane debonding behavior at the interface. Note that a traction-separation model identified in this way will be based on the assumption of linear elasticity, thus, neglects the presence of a process zone. The two critical data points required for this purpose are (1) the normal interfacial strength (i.e., the local stress corresponding to the debonding initiation)

and its corresponding interface opening displacement, and (2) the opening displacement that corresponds to complete separation, i.e., zero traction. Data sets required for the first part were presented in previous sections. Data needed for the second part practically includes one value representing the maximum displacement at complete separation. The latter can be identified indirectly from the graphs shown earlier in Following the procedure described above, a bilinear traction-separation model is constructed and shown in Figure 9. Note that the scatter bars in this figure represent the variation in displacement data extracted from the top and bottom sides of the interface, again suggesting a symmetric debonding behavior. In addition to the applicability of the identified traction-separation model in FE analyses, the model also enables the estimation of critical ERR, determined as the area under the traction-separation profile. As such, the nominal energy release rate associated with the in-plane fiber-matrix interface debonding in this work is determined as 25.45 J/m 2 . We should emphasize that the traction-separation parameters quantified herein are based on a linear elastic stress solution. Therefore, the model determined in this work is an upper bound of the actual traction-separation profile. Consequently, an upper bound for the interface critical ERR has also been quantified. In practice, due to the nonlinear constitutive response of the constituents and the interface, collectively leading to the formation of a process zone, the interfacial stresses may be bounded by values smaller than those obtained based on a linear elastic solution. It is worth noting that the tensile bond strengths (max. traction in Figure 9) determined here are lower than those reported in previous studies [START_REF] Ageorges | Single-fibre Broutman test: Fibre-matrix interface transverse debonding[END_REF][START_REF] Koyanagi | Mixed-mode interfacial debonding simulation in single-fiber composite under a transverse load[END_REF][START_REF] Ogihara | Evaluation of interfacial tensile strength in glass fiber/epoxy resin interface using the cruciform specimen method[END_REF][START_REF] Totten | In situ determination of the fiber-matrix interface tensile strength[END_REF]. The main sources of difference are identified as (1) the absence of fiber sizing in the system examined in this work, as further discussed in Sec. 5; and (2) the method of extraction of the bond strength in previous reports, i.e., computational and/or shear-dominated states of stress versus the normal stress state discussed here. In contrast, the critical ERR determined by the present approach is higher than those reported elsewhere [START_REF] Koyanagi | Mixed-mode interfacial debonding simulation in single-fiber composite under a transverse load[END_REF][START_REF] Rädel | Peridynamic analysis of fibre-matrix debond and matrix failure mechanisms in composites under transverse tensile load by an energy-based damage criterion[END_REF][START_REF] Varna | Transverse single-fibre test for interfacial debonding in composites: 2[END_REF], again due to the interfacial characteristics of the present system and the different approaches utilized in those previous reports.

Experimental Characterization of Interface Debonding Growth

While not essential to the characterization of debonding behavior discussed in previous sections, experimental results obtained here can be further utilized to characterize the interface debonding growth as a function of global deformation/stress. This information can be valuable in the validation of modeling analyses that are abundantly available in the literature without such in-depth experimental evidence.

Three separate approaches were utilized to characterize the propagation of the debonded interface in the single fiber specimens in this work. The first method used is the visual determination of the debonding angle, performed independently by two operators. The determination of the debonding angle is possible thanks to the in-situ optical measurement and the raw images obtained. Due to the loading configuration, debonding initiation occurs simultaneously at both poles. Each picture is imported into an image processing program (ImageJ, National Institutes of Health, USA) and zoomed in to observe the interface at the upper pole. Upon the debonding initiation, the formation of a debonded black area allows the identification of the two tips of the crack, enabling the measurement of debonding angle. Nevertheless, precise crack tip identification remains a challenge since the debonding opening becomes smaller when approaching the crack tip.

Another strategy to estimate the debonding angle consists of using DIC analyses. The image correlation can no longer be achieved locally in the presence of a black, previously nonexistent debonded zone. This results in the development of high correlation residuals. Therefore, the debonding angle can be determined by applying a threshold on the correlation residuals. The threshold is required so that no debonding is obtained for small loading levels before crack initiation. A DIC residual threshold equal to 1% of the grey level was chosen here. Figures 10a and10b show qualitative comparisons between the debonding angle measurements using the visual and the DIC residual methods.

The third method utilized for the determination of debonding angle consists of combining the simulation and experimental results. In this method (denoted as 'FEA' in Figure 10c), the debonding angle was determined by matching the angle at which the experimentally measured and numerically evaluated values of Δ t equated for a given applied stress. The variation of debonding angle with the applied global stress, determined using the three approaches, is shown in Figure 10c. A significantly smaller debonding angle after initiation is determined based on the FEA approach than the other two measurement methods. This difference could be explained by the consequential role of the aforementioned process zone, which is not considered in the calculation of Δ t but could induce a sufficiently large displacement jump to be detected visually or by DIC. On the other hand, at crack initiation, the difference between the debonding angle measurements by DIC residual and visual approaches remains smaller than 16 degrees (20% relative error), which enables the determination of the debonding initiation angle by accounting for measurement uncertainty.

Although some of the main limitations associated with the experimental nature of the study are presented in the following section, there are a few fundamental issues that must be addressed when using approaches like those presented herein. Perhaps the most significant of these challenges is the lack of a universal definition for debonding length, as it relates to the nature of the modeling analyses conducted. In particular, the modeling results can be prone to significant variations when considering linear elastic fracture mechanics (LEFM) with no process zone or a cohesive zone model. In addition, the accurate identification of the debonding 'tip' is not straightforward. The latter, in turn, will add to the challenging task of defining a debonding length. Finally, from a practical perspective, what can be measured by DIC are the displacement fields around the fiber-matrix interface, which is a consequence of a debonding state, possibly in the presence of a process zone (e.g., see the strain distribution along the debonded interface in Supplementary Information). The authors anticipate that a coupled experimental-simulations approach could be more realistic if direct comparisons are made on the displacement/strain fields rather than on a subjective metric, i.e., debonding length or angle. Nonetheless, while only a first step towards the quantitative experimental characterization of fiber-matrix debonding, the results presented in this work highlight the significance of coupled image-based and FE analyses in the understanding of the process.

Current Limitations and Future Directions

The most significant limitation of the current experimental approach is its inability to provide information on the out-of-plane deformation behavior of the matrix in the fiber vicinity. Previous reports have shown that highly localized out-of-plane deformations can occur in areas near the fiber-matrix interface due to Poisson's effect, especially in the case of ductile matrix materials [START_REF] Tabiai | In-situ full field out of plane displacement and strain measurements at the micro-scale in single reinforcement composites under transverse load[END_REF]. While the epoxy matrix used herein was brittle enough to mitigate potential errors due to such out-of-plane deformations, we acknowledge the absence of out-of-plane deformation characteristics in the present analyses. Efforts are currently underway to develop techniques that allow the measurement of surface displacements and strains in three dimensions (using stereo DIC). To take it to the next step, a single camera viewing the specimen thickness could also be used to characterize tunneling of the debonded interface along the fiber length (see schematic in

Supplementary Information

). The latter can only be possible by taking advantage of the semi-translucent nature of the matrix and fibers that lead to light refraction upon crack tunneling [START_REF] Koyanagi | Mixed-mode interfacial debonding simulation in single-fiber composite under a transverse load[END_REF][START_REF] Tandon | Evaluation of interfacial normal strength in a SCS-0/epoxy composite with cruciform specimens[END_REF]. An alternative approach will be the application of in-situ microcomputed tomography (micro-CT) [START_REF] Martyniuk | 3D in situ observations of glass fiber/matrix interfacial debonding[END_REF].

Another limitation of the current study originates from the application of unsized glass macro fibers.

The unsized fibers used in this study led to weak fiber-matrix interface strengths [START_REF] Thomason | Glass fibre sizing: A review[END_REF]. The consequence of such weak fiber-matrix interfaces is manifested in the lack of concrete evidence for crack kink-out and matrix cracking mechanisms.

Finally, the proposed approach could also be applied to estimate a traction-separation profile that accounts for mode mixity by evaluating the normal and tangential displacement jumps at several locations along the fiber-matrix interface, eventually allowing for the concurrent estimation of the interface shear strength and critical energy release rate.

Conclusions

Single fiber model composites were fabricated using glass macro fibers embedded in an epoxy matrix. Induced by the application of a far-field tensile load, strain fields developed in the vicinity of the fiber-matrix interface were characterized using high magnification optical digital image correlation. Local strains and debonded interface propagation characteristics were quantified. It was shown that the strain evolution in certain locations around the fiber-matrix interface can be used as effective indicators for the interface debonding initiation. Coupled with a theoretical solution for the interfacial stress model, the local inplane stress associated with the debonding initiation was identified. The local strains, interface opening displacement, and stresses were used as input to develop a simple bi-linear traction-separation model for the composite system examined. The visual characterization of debonding evolution was performed and supplemented by results obtained from simulations. It was shown that the experimental quantifications (consisting of purely visual and DIC residual methods) lead to relatively consistent measures for the debonding angle and stress-dependent growth of the debonded interface. On the other hand, the finite element simulations predict a significantly smaller debonding initiation angle. The approach presented here provides a systematic experimental methodology for the quantitative characterization of fiber-matrix interface debonding mechanisms in single fiber composite specimens.

S2. FEA Informed Specimen Size Determination Based on Energy Release Rate (ERR)

The ERR is defined as the opposite of the derivative of the elastic strain energy with respect to the crack surface. The elastic strain energy is calculated as 𝑊𝑊 = 1 𝑉𝑉 ∫ 𝜎𝜎 𝑉𝑉 : 𝜀𝜀𝑑𝑑𝑑𝑑 for several debonding angles by successively releasing nodes from the fiber pole at the fiber-matrix interface. This calculation is performed for several w/r ratios, which enables computing the relative differences between the calculated ERRs and that obtained under the assumption of a fiber embedded in an infinite matrix [3]. Figure S2 shows the variation of the relative difference (in percentage) as a function of w/r ratio, highlighting a relative error smaller than 2%

for the w/r ratios used in this work.

Figure 5 .

 5 Figure 5. Local strain variations shown in Figure 5 indicate an initial rise, leading to a peak at the global stress of 5.35 MPa. Upon the progression of debonding beyond this critical stress, the strain curves show a decline up to global stresses of ca. 8.82 MPa, followed by a constant increase until the final failure. These trends can also be inferred from the strain derivative (𝑑𝑑𝜀𝜀 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 /𝑑𝑑𝜎𝜎 ∞ ) data shown in the top panel of the same figure. The second critical stress (i.e., global stress of 8.82 MPa) corresponds to the instant at which the unloading of the local area of interest has been completed, thus, representing the end of the interface separation process. The opening displacement corresponding to this stress value can be readily extracted from the data shown in Figure 7.

FiguresFigure 1 -

 1 Figures

Figure 2 -

 2 Figure 2-(a) Schematic of the single-fiber test specimen. Tensile specimen thickness: 6 mm. (b) speckle patterned area of interest with its grayscale intensity shown in (c). (d) Tensile stress-strain behavior of the epoxy matrix. Scatter bars indicate the range of variation within three independent measurements.

Figure 3 -

 3 Figure 3-Evolution of in-plane strain components, (a) εyy, (b) εxx, and (c) εxy with global stress. Tensile load is applied in y-direction. All scale bars represent 2 mm.

Figure 4 - 2 :

 42 Figure 4-Schematic sequence of the fiber-matrix interfacial debonding process. 1: Initial (stress-free) state, 2: After the application of a tensile load and before the onset of debonding, and 3: Upon further increase of the tensile load and the formation of a debonded interface. The debonding leads to the partial unloading of the matrix located near the interface, as shown in an exaggerated fashion by the darker rectangular regions in 2 and 3.

Figure 5 -

 5 Figure 5-Variations of global strain (bottom), local longitudinal strains (middle) and

Figure 6 -

 6 Figure 6-Variation of (a) 𝑑𝑑𝜀𝜀 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 / 𝑑𝑑𝜎𝜎 ∞ as a function of the debonding angle, and (b) Local longitudinal strain as a function of global stress for various debonding angles, comparing DIC (experimental) measurements and FE simulation results.

Figure 7 -

 7 Figure 7-(a) Schematic view and (b) experimentally measured variation of vertical gauge displacement, Δ t , across the fiber-matrix interface measured by an optical extensometer with an initial length of 𝐿𝐿 0 t . 𝐿𝐿 0 f and 𝐿𝐿 0 m

Figure 8 -

 8 Figure 8-Vertical gauge displacement as a function of debonding angle for three extensometer gauge lengths, 𝐿𝐿 0 t , of 320, 440, and 660 µm. Interface opening displacement due to debonding, δ n , is also included for reference.

Figure 9 -

 9 Figure 9-Traction-separation model identified for the in-plane normal debonding between the glass macro fiber and epoxy matrix in this work. Horizontal scatter bars represent the variations in the top and bottom poles of the interface.

Figure 10 -

 10 Figure 10-Determination of the debonding angle by (a) visual and (b) DIC residual approaches. The global stress corresponding to both images is approximately 7 MPa. The debonding angle in (b) is identified based on a 1% threshold. (c) Variation of debonding angle with global stress, determined using three different approaches: visual, DIC residual, and debonding opening (FEA) calculations. The average of the first two approaches is shown by a continuous dashed line.

Figure S1 -

 S1 Figure S1-(a) Schematic of the 2D finite element model used in this study. (b) Mesh sensitivity analysis showing the variation of normalized local strain with respect to global strain with the number of degrees of freedom (DOF).

Table 1 -

 1 Mechanical properties of the matrix and fiber materials examined in this work. Determined using 0.2% offset method; b Based on at least 3 independent in-house measurements; c Data provided by the manufacturer

	Tables
	Material	Young's	Poisson's	Yield strength	Tensile strength	Failure
		Modulus	ratio	(MPa) a	(MPa)	strain, εf
		(GPa)
	Epoxy b	2.24±0.10	0.40±0.04	31.30±1.11	45.25±0.74	0.035±0.001
	Glass Fiber c	63	0.3	-	-	-
	a

Available on: H. Girard et al., Composites Pat A 2023; 171: 107573 https://doi.org/10.1016/j.compositesa.2023.107573

Acknowledgments B.K. gratefully acknowledges the financial support from the Advanced Materials and Manufacturing Institute (AMMI) at Rowan University. This work was partially supported by Rowan University faculty start-up funds.

Supplementary Information

Additional details regarding the finite element model, the sensitivity of ERR to specimen dimensions, strain distribution along the debonded interface, stress-strain response of the single-fiber specimens, and proposed experimental setup for the concurrent measurement of debonding and tunneling are provided as Supplementary Information.

Supplementary Information

Experimental Characterization of In-Plane Debonding at Fiber-Matrix Interface Using Single Glass Macro Fiber Samples

Hugo Girard a, † , Behrad Koohbor

S1. Finite Element Analyses

Figure S1a illustrates the 2D plane stress model created for finite element analyses. The geometry consists of a fiber of radius r embedded in a matrix (fiber to free edge distance w, length 30r to ensure a prescribed loading at infinity) and is meshed with four-node linear elements. Mesh size sensitivity analysis was performed based on a 1% convergence of strain fields criterion (Figure S1b). Dirichlet boundary conditions were prescribed on the top and bottom side of the domain, corresponding to a remotely applied stress (𝜎𝜎 ∞ ).

The influence of the fiber to free edge distance was studied by computing the stress fields for several w/r ratios. The stress at the fiber-matrix interface was then compared to the analytical solution corresponding to a fiber embedded in an infinite elastic matrix [1] or to similar FE configurations to evaluate the perturbation induced by the free edges.

The assumption of plane stress state is valid until debonding initiation. The debonding initiation leads to the generation of stress singularity and a 3D state of stress. Therefore, upon debonding initiation, the assumption of plane stress can no longer be valid. Nevertheless, a 2D analysis can provide reliable results, as verified by the comparison of critical stress intensity factors obtained from 3D and 2D (surface analyses) discussed in Doitrand et al. [2].

S3. Stress-Strain Response of Single Fiber Specimens

The stress-strain response of the single-fiber specimens, obtained from two independent measurements, is shown in Figure S3. The curves are slightly different than those of the neat epoxy (see Figure 2 in the main article). Due to the debonding-induced stress singularity, the ultimate strength of the single-fiber specimens is slightly lower than that of the epoxy material.

S4. Strain Distribution along the Debonded Interface

Figure S4 shows the Variation of in-plane strain components, εyy, εxx, and εxy, along the upper half of the interface and at different global stress values. Note that the local strains in this graph are extracted from a curved line segment located 250 µm away from the interface inside the matrix. The local strains show significant variations at stresses approaching the failure strength of the specimens. However, at lower global stresses, no discernible variation is observed in the local strain curves. This observation suggests that tracking of local strain alone may not be sufficient for an accurate assessment of the debonding propagation process.

Another point worth mentioning here is the apparent asymmetry in the strain curves. Symmetry in strain patterns is expected where (1) the fiber is embedded precisely at the center of the dog bone tensile specimens,

(2) the tensile load is perfectly uniaxial and parallel to the specimen mid-plane, and (3) before debonding initiation. The experimental nature of this study makes it impractical to place the fiber exactly at the centroid of the tensile test pieces, as also evidenced by the slight asymmetry in the strain contour plots depicted in Figure 3 in the main article. The slight asymmetry in the strain distributions herein most likely originates from the inevitable imperfections associated with the experimental nature of the work. In addition, the symmetry in the strain field cannot be guaranteed after debonding.

S5. Experimental Setup for Concurrent Measurement of Debonding and Tunneling

The experimental setup proposed for the concurrent measurement of fiber-matrix interface debonding and tunneling is shown schematically in Figure S5. The application of a dual imaging system allows for the identification of full-field deformation distribution (including out-of-plane motions) on the specimen front.

The side-view imaging will facilitate the coupling of in-plane debonding with the tunneling effect along the length of the single macro fiber.