#  On Some Estimate for the Norm
of an Interpolation Projector

Mikhail Nevskii111 Department of Mathematics, P.G. Demidov Yaroslavl State
University,
Sovetskaya str., 14, Yaroslavl, 150003, Russia,
<EMAIL_ADDRESS>orcid.org/0000-0002-6392-7618

(May 7, 2022)

###### Abstract

Let $Q_{n}=[0,1]^{n}$ be the unit cube in ${\mathbb{R}}^{n}$ and let
$C(Q_{n})$ be a space of continuous functions $f:Q_{n}\to{\mathbb{R}}$ with
the norm $\|f\|_{C(Q_{n})}:=\max_{x\in Q_{n}}|f(x)|.$ By
$\Pi_{1}\left({\mathbb{R}}^{n}\right)$ denote a set of polynomials of degree
$\leq 1$, i. e., a set of linear functions on ${\mathbb{R}}^{n}$. The
interpolation projector $P:C(Q_{n})\to\Pi_{1}({\mathbb{R}}^{n})$ with the
nodes $x^{(j)}\in Q_{n}$ is defined by the equalities
$Pf\left(x^{(j)}\right)=f\left(x^{(j)}\right)$, $j=1,$ $\ldots,$ $n+1$. Let
$\|P\|_{Q_{n}}$ be the norm of $P$ as an operator from $C(Q_{n})$ to
$C(Q_{n})$. If $n+1$ is an Hadamard number, then there exists a nondegenerate
regular simplex having the vertices at vertices of $Q_{n}$. We discuss some
approaches to get inequalities of the form $||P||_{Q_{n}}\leq c\sqrt{n}$ for
the norm of the corresponding projector $P$.

Keywords: Hadamard matrix, regular simplex, linear interpolation,projector,
norm

MSC: 41A05, 52B55, 52C07

##  1 Introduction

Let $K$ be a convex body in ${\mathbb{R}}^{n}$. i. e., a compact convex subset
of ${\mathbb{R}}^{n}$ with nonempty interior. Denote by $C(K)$ a space of
continuous functions $f:K\to{\mathbb{R}}$ with the uniform norm

$\|f\|_{C(K)}:=\max\limits_{x\in K}|f(x)|.$

By $\Pi_{1}\left({\mathbb{R}}^{n}\right)$ we mean a set of polynomials in $n$
variables of degree $\leq 1$, i. e., of linear functions on
${\mathbb{R}}^{n}$. For $x^{(0)}\in{\mathbb{R}}^{n},R>0$, by $B(x^{(0)};R)$ we
denote the $n$-dimensional Euclidean ball given by the inequality
$\|x-x^{(0)}\|\leq R$. Here

$\|x\|:=\sqrt{(x,x)}=\left(\sum\limits_{i=1}^{n}x_{i}^{2}\right)^{1/2}.$

By definition, $B_{n}:=B(0;1),Q_{n}:=[0,1]^{n},Q_{n}^{\prime}:=[-1,1]^{n}$.

Let $S$ be a nondegenerate simplex in ${\mathbb{R}}^{n}$ with the vertices
$x^{(j)}=\left(x_{1}^{(j)},\ldots,x_{n}^{(j)}\right),$ $1\leq j\leq n+1.$
Consider the following vertex matrix of this simplex:

${\bf S}:=\left(\begin{array}[]{cccc}x_{1}^{(1)}&\ldots&x_{n}^{(1)}&1\\\
x_{1}^{(2)}&\ldots&x_{n}^{(2)}&1\\\ \vdots&\vdots&\vdots&\vdots\\\
x_{1}^{(n+1)}&\ldots&x_{n}^{(n+1)}&1\\\ \end{array}\right).$

Let ${\bf S}^{-1}$ $=(l_{ij})$. Linear polynomials
$\lambda_{j}(x):=l_{1j}x_{1}+\ldots+l_{nj}x_{n}+l_{n+1,j}$ have the property
$\lambda_{j}\left(x^{(k)}\right)$ $=$ $\delta_{j}^{k}$. We call $\lambda_{j}$
the basic Lagrange polynomials corresponding to $S$. For an arbitrary
$x\in{\mathbb{R}}^{n}$,

$x=\sum_{j=1}^{n+1}\lambda_{j}(x)x^{(j)},\quad\sum_{j=1}^{n+1}\lambda_{j}(x)=1.$

These equalities mean that $\lambda_{j}(x)$ are the barycentric coordinates of
$x$. Also we have $\lambda_{j}(x)=\frac{\Delta_{j}(x)}{\Delta}$, where
$\Delta=\det(\bf S)$ and $\Delta_{j}(x)$ is obtained from $\Delta$ by
replacing the $j$th row by the row $(x_{1}\ \ldots\ x_{n}\ 1)$. For details,
see [12, §1.1].

We say that an interpolation projector $P:C(K)\to\Pi_{1}({\mathbb{R}}^{n})$
corresponds to a simplex $S\subset K$ if the nodes of $P$ coincide with the
vertices of $S$. This projector is defined by the equalities
$Pf\left(x^{(j)}\right)=f\left(x^{(j)}\right).$ The following analogue of the
Lagrange interpolation formula holds:

$Pf(x)=\sum\limits_{j=1}^{n+1}f\left(x^{(j)}\right)\lambda_{j}(x).$ (1)

Denote by $\|P\|_{K}$ the norm of $P$ as an operator from $C(K)$ in $C(K)$.
From (1), it follows that

$\|P\|_{K}=\max_{x\in K}\sum_{j=1}^{n+1}|\lambda_{j}(x)|=\max_{x\in
K}\sum_{j=1}^{n+1}\frac{|\Delta_{j}(x)|}{|\Delta|}.$ (2)

If $K$ is a convex polytope, then the maxima in (2) also can be taken only
over $x\in{\rm ver}(K)$, where ${\rm ver}(K)$ denotes the vertex set of $K$.

In this paper, we consider the case when $n+1$ is an Hadamard number,i. e.,
there exists an Hadamard matrix of order $n+1$. An Hadamard matrix of order
$m$ is a square binary matrix $\bf H$ with entries either $1$ or $-1$ which
satisfies

${\bf H}^{-1}=\frac{1}{m}\,{\bf H}^{T}.$

This means that the rows of ${\bf H}$ are pairwise orthogonal, with respect to
the standard scalar product on ${\mathbb{R}}^{m}.$

The order of an Hadamard matrix is 1 or 2 or some multiple of 4 (see [4]).But
it is still unknown whether an Hadamard matrix exists for every order of the
form $m=4k$. This is one of the longest lasting open problems in Mathematics
called the Hadamard matrix conjecture. The orders below 1500 for which
Hadamard matrices are not known are 668, 716, 892, 956, 1132, 1244, 1388, and
1436 (for links, see, e.g., [5], [8]).

Two Hadamard matrices are called equivalent iff one can be obtained from
another by performing a finite sequence of the following operations:
multiplication of some rows or columns by $-1$, or permutation of rows or
columns. Up to equivalence, there is a unique Hadamard matrix of orders 1, 2,
4, 8, and 12. There are 5 equivalence classes of Hadamard matrices of order
16, 3 of order 20, 60 of order 24, and 487 of order 28. For orders 32, 36, and
40, the number of equivalence classes is much greater. For $n=32$ there are at
least 3,578,006 equivalence classes; for $n=36$, at least 4,745,357 (see [5]).

Let $Q$ be an arbitrary $n$-dimensional cube. If $n+1$ is an Hadamard number,
then there exists a nondegenerate regular simplex having the vertices
coinciding with vertices of $Q$. In the paper, we discuss some approaches to
obtain estimates of the form $||P||_{Q}$ $\leq$ $c\sqrt{n}$ for the norm of
the corresponding interpolation projector. For convenience, we will take
$Q=Q_{n}$ or $Q=Q_{n}^{\prime}$.

##  2 Making Use of the Numbers $h_{n}$

Denote by $h_{n}$ the maximum value of a determinant of order $n$ with entries
$0$ or $1$. Let $\nu_{n}$ be the maximum volume of an $n$-dimensional simplex
contained in $Q_{n}$. These numbers satisfy the equality $h_{n}=n!\nu_{n}$
(see [6]). For any $n$, there exists in $Q_{n}$ a maximum volume simplex with
some vertex coinciding with a vertex of the cube. For such simplices, we have
the following theorem.

Theorem 1. Suppose $S\subset Q_{n}$ is a maximum volume simplex with some
vertex coinciding with a vertex of the cube, $P$ is the corresponding
interpolation projector. Then

$\|P\|_{Q_{n}}\leq\frac{2h_{n+1}}{h_{n}}+1.$ (3)

Proof. Let $x^{(1)},\ldots,x^{(n+1)}$ be the vertices of $S$. We may assume
that $x^{(n+1)}=0,$ i. e., the vertex matrix of the simplex has the form

${\bf S}=\left(\begin{array}[]{cccc}x_{1}^{(1)}&\ldots&x_{n}^{(1)}&1\\\
\vdots&\vdots&\vdots&\vdots\\\ x_{1}^{(n)}&\ldots&x_{n}^{(n)}&1\\\
0&\ldots&0&1\\\ \end{array}\right).$

Since $S$ is a maximum volume simplex in $Q_{n}$, we have ${\rm
vol}(S)=\nu_{n}$ and $|\Delta|=|\det({\bf S})|=n!\nu_{n}=h_{n}.$ Hence, if
$x\in{\rm ver}(Q_{n})$, then $|\Delta_{j}(x)|\leq|\Delta|$, $1\leq j\leq n+1.$

Fix $x\in{\rm ver}(Q_{n})$ and consider determinants $\Delta_{j}(x)$, $1\leq
j\leq n$. From definition of $\Delta_{j}(x)$ and properties of determinant,

$|\Delta_{j}(x)|=\left|\begin{array}[]{cccc}x_{1}^{(1)}&\ldots&x_{n}^{(1)}&1\\\
\ldots&\ldots&\ldots&\ldots\\\ x_{1}&\ldots&x_{n}&1\\\
\ldots&\ldots&\ldots&\ldots\\\ x_{1}^{(n)}&\ldots&x_{n}^{(n)}&1\\\
0&\ldots&0&\pm 1\\\
\end{array}\right|=\left|\begin{array}[]{cccc}x_{1}^{(1)}&\ldots&x_{n}^{(1)}&0\\\
\ldots&\ldots&\ldots&\ldots\\\ x_{1}&\ldots&x_{n}&0\\\
\ldots&\ldots&\ldots&\ldots\\\ x_{1}^{(n)}&\ldots&x_{n}^{(n)}&0\\\
0&\ldots&0&\pm 1\\\ \end{array}\right|=$

$=\left|\begin{array}[]{cccc}x_{1}^{(1)}&\ldots&x_{n}^{(1)}&0\\\
\ldots&\ldots&\ldots&\ldots\\\ x_{1}&\ldots&x_{n}&0\\\
\ldots&\ldots&\ldots&\ldots\\\ x_{1}^{(n)}&\ldots&x_{n}^{(n)}&0\\\
x_{1}^{(j)}&\ldots&x_{n}^{(j)}&\pm 1\\\
\end{array}\right|=\left|\begin{array}[]{cccc}x_{1}&\ldots&x_{n}&0\\\
x_{1}^{(1)}&\ldots&x_{n}^{(1)}&0\\\ \ldots&\ldots&\ldots&\ldots\\\
x_{1}^{(j)}&\ldots&x_{n}^{(j)}&\pm 1\\\ \ldots&\ldots&\ldots&\ldots\\\
x_{1}^{(n)}&\ldots&x_{n}^{(n)}&0\\\ \end{array}\right|.$

Here we mark the $j$th row. Consequently, there exist numbers $u_{j}=\pm 1$
such that

$\sum_{j=1}^{n}|\Delta_{j}(x)|=\left|\begin{array}[]{cccc}x_{1}&\ldots&x_{n}&0\\\
x_{1}^{(1)}&\ldots&x_{n}^{(1)}&u_{1}\\\ \vdots&\vdots&\vdots&\vdots\\\
x_{1}^{(n)}&\ldots&x_{n}^{(n)}&u_{n}\end{array}\right|.$

Define $n$-dimensional binary vectors $v$ and $w$ as follows. If $u_{j}=1$,
then $v_{j}=1$ and $w_{j}=0$; if $u_{j}=-1,$ then $v_{j}=0$ and $w_{j}=1$. We
have $u=v-w$, hence

$\sum_{j=1}^{n}|\Delta_{j}(x)|=\left|\begin{array}[]{cccc}x_{1}&\ldots&x_{n}&0\\\
x_{1}^{(1)}&\ldots&x_{n}^{(1)}&v_{1}\\\ \vdots&\vdots&\vdots&\vdots\\\
x_{1}^{(n)}&\ldots&x_{n}^{(n)}&v_{n}\end{array}\right|\quad-\quad\left|\begin{array}[]{cccc}x_{1}&\ldots&x_{n}&0\\\
x_{1}^{(1)}&\ldots&x_{n}^{(1)}&w_{1}\\\ \vdots&\vdots&\vdots&\vdots\\\
x_{1}^{(n)}&\ldots&x_{n}^{(n)}&w_{n}\end{array}\right|.$

Denote these determinants of order $n+1$ by $d_{1}$ and $d_{2}$. Since their
entries are in $[0,1]$, then $|d_{1}|,|d_{2}|\leq h_{n+1}$,

$\sum_{j=1}^{n}|\Delta_{j}(x)|=d_{1}-d_{2}\leq|d_{1}|+|d_{2}|\leq 2h_{n+1}.$

Using $|\Delta|=h_{n},$ we get

$\sum_{j=1}^{n}\frac{|\Delta_{j}(x)|}{|\Delta|}\leq\frac{2h_{n+1}}{h_{n}}.$

It remains to include the case $j=n+1$. Recall that

$\frac{|\Delta_{n+1}(x)|}{|\Delta|}\leq 1,\quad x\in{\rm ver}(Q_{n}),$

therefore, for any vertex of the cube,

$\sum_{j=1}^{n+1}\frac{|\Delta_{j}(x)|}{|\Delta|}\leq\frac{2h_{n+1}}{h_{n}}+1.$

We obtain

$\|P\|_{Q_{n}}=\max\limits_{x\in{\rm
ver}(Q_{n})}\sum\limits_{j=1}^{n+1}|\lambda_{j}(x)|=\max\limits_{x\in{\rm
ver}(Q_{n})}\sum\limits_{j=1}^{n+1}\frac{|\Delta_{j}(x)|}{|\Delta|}\leq\frac{2h_{n+1}}{h_{n}}+1.$

$\Box$

Corollary 1.

$\theta_{n}\leq\frac{2h_{n+1}}{h_{n}}+1.$ (4)

This follows immediately from (3).

Corollary 2. There exist $c,c_{1}>0$ does not depend on $n$ such that

$\frac{h_{n+1}}{h_{n}}\geq c\sqrt{n},\quad\frac{\nu_{n}}{\nu_{n+1}}\leq
c_{1}\sqrt{n}.$ (5)

Proof. The author proved that $\theta_{n}\geq\gamma\sqrt{n}$, $\gamma>0$ not
depends on $n$ (see [12]). Hence, the left-hand inequality in (5) follows from
(4). Since $h_{n}=n!\nu_{n}$, the right-hand inequality appears from the left-
hand one. $\Box$

Let us note more explicit forms of (4). For any $n$,

$\theta_{n}>\frac{\sqrt{n-1}}{e}.$

Consequently,

$\frac{h_{n+1}}{h_{n}}\geq\frac{\theta_{n}-1}{2}>\frac{\sqrt{n-1}}{2e}-\frac{1}{2},$
$\frac{\nu_{n+1}}{\nu_{n}}=\frac{h_{n+1}}{(n+1)h_{n}}>\frac{\sqrt{n-1}}{2e(n+1)}-\frac{1}{2(n+1)}.$

The above proof of Theorem 1 follows the approach of [10]. Unfortunately,
later it was noticed that the arguments in [10] perfomed for an arbitrary $n$
have a lacuna. Yet in the case when $n+1$ is an Hadamard number Theorem 1
allows to obtain for $\theta_{n}$ upper bounds which cannot be improved in
order of $n$. To show this, let us utilize the estimates

$h_{n}\leq\frac{\left(n+1\right)^{(n+1)/2}}{2^{n}},\quad n\in{\mathbb{N}};$
(6) $h_{n}\leq\frac{n^{n/2}\sqrt{2n+1}}{2^{n}},\quad n\ \mbox{ is even}.$ (7)

The equality in (6) takes place if and only if $n+1$ is an Hadamard number.
Relation (6) was proved by Hadamard [3]; inequality (7) was obtained by Barba
[1] (see [6]).

Corollary 3. Suppose $n+1$ is an Hadamard number. If $S$ is a regular simplex
with the vertices coinciding with vertices of $Q_{n},$ then for the
corresponding interpolation projector $P$

$\|P\|_{Q_{n}}\leq\sqrt{2n+3}+1.$

Proof. If $n+1$ is an Hadamard number, then the considered simplex $S$ has
maximum volume among all simplices contained in $Q_{n}$. Hence, we can apply
to $S$ Theorem 1:

$\|P\|_{Q_{n}}\leq\frac{2h_{n+1}}{h_{n}}+1.$

As it was noted above, in this case

$h_{n}=\frac{\left(n+1\right)^{(n+1)/2}}{2^{n}}.$

Since $n+1$ is even, we can apply to $h_{n+1}$ inequality (7):

$h_{n+1}\leq\frac{(n+1)^{(n+1)/2}\sqrt{2n+3}}{2^{n+1}}.$

From this, it follows

$\|P\|_{Q_{n}}\leq\frac{2h_{n+1}}{h_{n}}+1\leq\frac{2(n+1)^{(n+1)/2}\sqrt{2n+3}}{2^{n+1}}\cdot\frac{2^{n}}{\left(n+1\right)^{(n+1)/2}}+1=$
$=\sqrt{2n+3}+1.$

$\Box$

Corollary 4. If $n+1$ is an Hadamard number, then
$\theta_{n}\leq\sqrt{2n+3}+1.$

Upper bounds given in Corollaries 3 and 4 are unimprovable in order of $n$,
but can be made more precise. We will show this in next two sections.

##  3 The Proof Based on Hadamard Matrices

Recall (see [6]) that $n+1$ is an Hadamard number iff it is possible to
inscribe an $n$-dimensional regular simplex into an $n$-dimensional cube in
such a way that the vertices of the simplex coincide with vertices of the
cube. This can be easily shown for the cube $Q_{n}^{\prime}=[-1,1]^{n}$. There
is a simple correspondence between Hadamard matrices of order $n+1$ with the
last column consisting of $1$’s and $n$-dimensional regular simplices whose
vertices are at the vertices of $Q_{n}^{\prime}$. Since we need this
correspondence, we will dwell on it in more detail.

If $S$ is a regular simplex of the specified type, then its vertex matrix
${\bf S}$ is an Hadamard matrix of order $n+1$. Indeed, let
$x^{(1)},\ldots,x^{(n+1)}$ be the vertices of $S$. Since $x^{(j)}$ coincide
with vertices of the cube, the entries of $\bf S$ are $\pm 1$ and its last
column consists of $1$’s. Denote by $h^{(j)}$ the rows of ${\bf S}$ and make
sure that these vectors are pairwise orthogonal in ${\mathbb{R}}^{n+1}$.
Equalities $||x^{(j)}||=\sqrt{n}$ mean that simplex $S$ is inscribed into an
$n$-dimensional ball of radius $\sqrt{n}$. The edge-length $d$ of a regular
simplex and radius $R$ of the circumscribed ball satisfy the relation

$d=R\sqrt{2}\sqrt{\frac{n+1}{n}}$ (8)

(see, e. g., [12]). If $R=\sqrt{n},$ then $d=\sqrt{2(n+1)}.$ Hence,

$2(n+1)=||x^{(j)}-x^{(k)}||^{2}=2||x^{(j)}||^{2}-2(x^{(j)},x^{(k)})=2n-2(x^{(j)},x^{(k)}).$

We see that $(x^{(j)},x^{(k)})=-1.$ Since the rows of $\bf S$ in the first $n$
columns contain coordinates of vertices of the simplex and the last element of
any row is $1$, then

$(h^{(j)},h^{(k)})=(x^{(j)},x^{(k)})+1=0.$

Thus, $h^{(j)}$ are pairwise orthogonal in ${\mathbb{R}}^{n+1}$. This means
that $\bf S$ is an Hadamard matrix of order $n+1$. We have

${\bf S}^{-1}=\frac{1}{n+1}\,{\bf S}^{T}.$ (9)

Conversely, let ${\bf H}$ be an Hadamard matrix of order $n+1$ with the last
column consisting of $1$’s. Consider the simplex $S$ whose vertex matrix is
$\bf H$. In other words, the vertices of $S$ are given by the rows of ${\bf
H}$ (excepting the last component). The simplex $S$ is inscribed into the cube
$Q_{n}^{\prime},$ moreover, ${\rm ver}(S)\subset{\rm ver}(Q_{n}^{\prime})$. If
$h^{(j)}$ are the rows of ${\bf H}$ and $x^{(j)}$ are the vertices of $S$,
then $(h^{(j)},h^{(k)})=0,$ $(x^{(j)},x^{(k)})=-1,$ $\|x^{(j)}\|^{2}=n.$ It
follows that $||x^{(j)}-x^{(k)}||^{2}=2(n+1)$. So, $S$ is a regular simplex
with the edge-length $\sqrt{2(n+1)}$.

Notice that any Hadamard matrix of order $n+1$ is equivalent to the vertex
matrix of some $n$-dimensional regular simplex inscribed in $Q_{n}^{\prime}$.
This vertex matrix can be obtained from the given Hadamard matrix after
multiplication of some rows by $-1$.

Theorem 2. Suppose $n+1$ is an Hadamard number, $S$ is an $n$-dimensional
regular simplex having the vertices at vertices of $Q_{n}^{\prime}$. Then for
the corresponding interpolation projector
$P:C(Q_{n}^{\prime})\to\Pi_{1}({\mathbb{R}}^{n})$

$\|P\|_{Q_{n}^{\prime}}\leq\sqrt{n+1}.$ (10)

Proof. Let $x^{(1)},\ldots,x^{(n+1)}$ be the vertices and let
$\lambda_{1},\ldots,\lambda_{n+1}$ be the basic Lagrange polynomials of the
simplex. Under our propositions, the vertex matrix $\bf S$ is an Hadamard
matrix of order $n+1$ with the last column consisting of $1$’s.

Let’s show that for $x\in{\mathbb{R}}^{n}$

$\sum_{j=1}^{n+1}\lambda_{j}(x)^{2}=\frac{||x||^{2}+1}{n+1}.$ (11)

Denote $y=(x_{1},\ldots,x_{n},1)\in{\mathbb{R}}^{n+1}.$ The coefficients of
$\lambda_{j}$ form the $j$th column of ${\bf S}^{-1}$. Since $\bf S$ is an
Hadamard matrix, it satisfies (9), whence

$\lambda_{j}(x)=\frac{1}{n+1}(h^{(j)},y).$

Further, $\displaystyle\frac{h^{(j)}}{\sqrt{n+1}}$ form an orthonormalized
basis of ${\mathbb{R}}^{n+1},$ therefore,

$y=\sum_{j=1}^{m}\frac{(h^{(j)},y)}{\sqrt{n+1}}\frac{h^{(j)}}{\sqrt{n+1}},\quad(y,y)=\sum_{j=1}^{n+1}\frac{(h^{(j)},y)^{2}}{n+1}.$

From this, we have

$||x||^{2}+1=(y,y)=\sum_{j=1}^{n+1}\frac{(h^{(j)},y)^{2}}{n+1}=(n+1)\sum_{j=1}^{n+1}\lambda_{j}(x)^{2}.$

Last equalities yield (11).

If $x$ is a vertex of $Q_{n}^{\prime}$, then $||x||^{2}=n$ and (11) implies
$\sum\lambda_{j}(x)^{2}=1.$ Applying Cauchy inequality, we get for $x\in{\rm
ver}(Q_{n}^{\prime})$

$\sum_{j=1}^{n+1}|\lambda_{j}(x)|\leq\left(\sum_{j=1}^{n+1}\lambda_{j}(x)^{2}\right)^{\frac{1}{2}}\cdot\sqrt{n+1}=\sqrt{n+1}.$

Suppose $P:C(Q_{n}^{\prime})\to\Pi_{1}({\mathbb{R}}^{n})$ is the interpolation
projector corresponding to $S$. Then we have

$||P||_{Q_{n}^{\prime}}=\max_{x\in{\rm
ver}(Q_{n}^{\prime})}\sum_{j=1}^{n+1}|\lambda_{j}(x)|\leq\sqrt{n+1}.$

Theorem is proved. $\Box$

By similarity reasons, Theorem 2 overcomes to any $n$-dimensional cube, for
instance, to the cube $Q_{n}=[0,1]^{n}$.

Corollary 5. Suppose $n+1$ is an Hadamard number. Then
$\theta_{n}\leq\sqrt{n+1}.$

The result of Corollary 5 is known (see [9], [12]), but the given proof is
more clearly related to Hadamard matrices.

Note that $n$-dimensional regular simplices with vertices at vertices of the
cube can be located differently with respect to vertices and faces of the
cube. This is quite evident when the norms of the corresponding projectors are
different. But this is also possible if regular simplices generate the same
norms. In more detail, we will describe an approach based on comparison of
$\mu$-vertices of the cube with respect to various simplices.

A notion $\mu$-vertex of the cube $Q_{n}$ with respect to a simplex $S\subset
Q_{n}$ was introduced by the author in [11]. If in place of $Q_{n}$ we
consider an arbitraty $n$-dimensional cube $Q$, we obtain equivalent results.

Assume $1\leq\mu\leq n.$ We say that $x\in{\rm ver}(Q)$ is a $\mu$-vertex of a
cube $Q$ with respect to a simplex $S\subset Q,$ iff for the interpolation
projector $P:C(Q)\to\Pi_{1}\left({\mathbb{R}}^{n}\right)$ with the nodes at
vertices of $S$ holds

$\sum_{j=1}^{n+1}|\lambda_{j}(x)|=\|P\|_{Q}$

and among numbers $\lambda_{j}(x)$ there are exactly $\mu$ negatives. This
notion is closely connected with some relation between $\|P\|_{Q}$ and the
value

$\xi(Q;S):=\min\\{\sigma\geq 1:Q\subset\sigma S\\}.$

Here $\sigma S$ denotes the homothetic copy of $S$ with the center of
homothety at the center of gravity of $S$ and the ratio of homothety $\sigma$.
We call $\xi(Q;S)$ the absorption index of the cube $Q$ by the simplex $S$.
For a projector $P:C(Q)\to\Pi_{1}\left({\mathbb{R}}^{n}\right)$ and the
corresponding simplex $S$, it is proved in [11] that

$\frac{n+1}{2n}\left(\|P\|_{Q}-1\right)+1\leq\xi(Q;S)\leq\frac{n+1}{2}\left(\|P\|_{Q}-1\right)+1.$
(12)

The right-hand equality in (12) takes place if and only if there exists an
$1$-vertex of $Q$ with respect to $S.$ If for some $\mu$ there is a
$\mu$-vertex of $Q$ with respect to $S,$ then

$\frac{n+1}{2\mu}\left(\|P\|_{Q}-1\right)+1\leq\xi(Q;S).$

Of course, simplices which have different sets of $\mu$-vertices with respect
to the containing cube, are differently located in the cube, even if have the
same projectors’ norms. These simplices are non-equivalent in the following
sense: one of them cannot be mapped into another by an orthogonal transform
which maps the cube into itself. Let us give some examples concerning
$Q=Q_{n}^{\prime}.$

While obtaining estimates for minimal norms of projectors, in [13] various
$n$-dimensional regular simplices arising from different Hadamard matrices of
order $n+1$ were calculated. For $n=15$, the order of matrices is equal to
$16$. Up to equivalence, there are exactly five these Hadamard matrices. They
correspond to five simplices described in Table 1. The results of calculations
were kindly delivered to the author by A. Yu. Ukhalov.

Table 1:  Regular simplices for $n=15$

$\begin{array}[]{|c|c|c|c|}\hline\cr S&{\|P\|}_{Q_{15}^{\prime}}&\mbox{\bf
Values }\ \mu&\mbox{\bf Number of}\ \mu\mbox{\bf-vertices}\\\ \hline\cr
S_{1}&4&6&m_{6}=448\\\ S_{2}&4&6&m_{6}=192\\\ S_{3}&4&6&m_{6}=64\\\
S_{4}&\frac{7}{2}&4,5,6,8&m_{4}=896,m_{5}=1344,m_{6}=5376,m_{8}=1344\\\
S_{5}&\frac{7}{2}&4,5,6,8&m_{4}=896,m_{5}=1344,m_{6}=5376,m_{8}=1344\\\
\hline\cr\end{array}$

By $m_{\mu}$ we denote a number of $\mu$-vertices of the cube $Q_{n}^{\prime}$
with respect to every simplex. For the rest $1\leq\mu\leq 15$, excepting the
given in Table 1, values $m_{\mu}$ are zero. Every simplex $S_{1}$, $S_{2}$,
and $S_{3}$ generates the same projector norm and has only $6$-vertices. But
the numbers of $6$-vertices for them are different, hence, these simplices are
pairwise non-equivalent. Each of them also is non-equivalent as to $S_{4}$, so
as to $S_{5}.$ The latter simplices generates equal norms and have the same
sets of $\mu$-vertices. Also we have $\theta_{15}\leq\frac{7}{2}=3.5$. This is
more exact than the estimate $\theta_{15}\leq 4$ of Corollary 5 for $n=15$.

Another example is related to $n=23.$ In [7], the results for 60 regular
simplices inscribed into $Q_{23}^{\prime}$ were mentioned. The simplices were
built from the available 60 pairwise non-equivalent Hadamard matrices of order
$24$. For all simplices, excepting ones with the numbers 16, 53, 59, and 60,
the projector norm equals $\frac{14}{3}=4.6666\ldots$, while for each of these
four simplices the projector norm is $\frac{9}{2}=4.5$. In particular, we have
$\theta_{23}\leq 4.5$ (what is noted also in [13]). This inequality is more
exact than the estimate $\theta_{23}\leq\sqrt{24}=4.8989\ldots$ of Corollary 5
for $n=23$. Each of the four exceptional simplices is not equivalent to any of
the 56 others.

Despite the possible differs, for all regular simplices with verices at
vertices of the cube, inequality (10) holds. The inscribed regular simplices
satisfying $\|P\|_{Q_{n}^{\prime}}=\sqrt{n+1}$ exist at least for $n=1$,
$n=3$, and $n=15.$ The problem of full description of dimensions $n$ with such
property is still open.

##  4 Connection with Interpolation on a Ball

A regular simplex inscribed into an $n$-dimensional ball has maximal volume
among all simplices contained in this ball. There are no another simplices
having this property (see [2], [15], and [16]). In the case when $n+1$ is an
Hadamard number, the similar property with respect to an $n$-dimensional cube
holds for a regular simplex inscribed into the cube. Furthermore, the upper
bound for projectors’ norm corresponding to these simplices both on a cube and
on a ball is equal to $\sqrt{n+1}$.

In this section, we give another proof of Theorem 2. It seems to be shorter
than previous one from Section 3, but is based on some relations which were
obtained in [14] in a rather complicated way. This approach is related with
interpolation on a Euclidean ball.

At first, suppose $S$ is a regular simplex inscribed into an $n$-dimensional
ball $B=B(x^{(0)};R)$ and $P:C(B)\to\Pi_{1}({\mathbb{R}}^{n})$ is the
corresponding interpolation projector. Obviously, $\|P\|_{B}$ does not depend
neither on the center $x^{(0)}$ and the radius $R$ of the ball, nor on the
chosen regular simplex inscribed into this ball. In other words, $\|P\|_{B}$
depends only on dimension $n$. For $0\leq t\leq n+1$, consider the function

$\psi(t):=\frac{2\sqrt{n}}{n+1}\Bigl{(}t(n+1-t)\Bigr{)}^{1/2}+\left|1-\frac{2t}{n+1}\right|.$
(13)

Let us denote $\displaystyle
a:=\left\lfloor\frac{n+1}{2}-\frac{\sqrt{n+1}}{2}\right\rfloor.$ It was proved
in [14] that

$\|P\|_{B}=\max\\{\psi(a),\psi(a+1)\\}.$

This equality yields

$\sqrt{n}\leq\|P\|_{B}\leq\sqrt{n+1}.$ (14)

Moreover, $\|P\|_{B}=\sqrt{n}$ only when $n=1$, and $\|P\|_{B}=\sqrt{n+1}$ if
and only if $\sqrt{n+1}$ is an integer.

Theorem 2 follows immediately from the right-hand inequality (14). Assume
$n+1$ is an Hadamard number. Consider an interpolation projector $P$ with the
nodes at those vertices of the cube $Q_{n}^{\prime}$ that form a regular
simplex $S$. Since $Q_{n}^{\prime}$ is inscribed into the unit ball $B_{n}$,
simplex $S$ is also inscribed into $B_{n}$. It remains to apply formula (2)
for the projector’s norm both on the cube and on the ball and also inequality
(14) in the case $B=B_{n}$:

$\|P\|_{Q_{n}^{\prime}}=\max_{x\in
Q_{n}^{\prime}}\sum_{j=1}^{n+1}|\lambda_{j}(x)|\leq\max_{x\in
B_{n}}\sum_{j=1}^{n+1}|\lambda_{j}(x)|=\|P\|_{B_{n}}\leq\sqrt{n+1}.$

Inequality (10) is proved.

The upper bound $\sqrt{n+1}$ of the projector’s norm is delivered in
interpolation on the ball or on the cube in different ways. For dimensions
$n=m^{2}-1$ , and only in these cases, the equality $\|P\|_{B_{n}}=\sqrt{n+1}$
holds for any regular simplex inscribed into $B_{n}$. If $n+1$ is an Hadamard
number, then the equality $\|P\|_{Q_{n}^{\prime}}=\sqrt{n+1}$ may hold as for
all regular simplices having vertices at vertices of the cube $(n=1,n=3)$, as
for some of them $(n=15),$ or may not be executed at all.

The connection with constructions on a ball can be also seen in the proof
given in Section 3 (see (8)).

## References

  * [1] G. Barba. Intorno al teorema di Hadamard sui determinanti a valore massimo, Glornale Mat. Battaglini (3), 71 (1933), 70–86.
  * [2] L. Fejes Tót. Regular Figures, New York: Macmillan/Pergamon, 1964.
  * [3] J. Hadamard. Résolution d’une question relative aux déterminants, Bull. Sciences Math. (2), 17 (1893), 240–246.
  * [4] M. Hall, Jr. Combinatorial Theory, Blaisdall Publishing Company, Waltham (Massachusetts)–Toronto–London, 1967.
  * [5] K. J. Horadam. Hadamard Matrices and Their Applications, Princeton University Press, 2007.
  * [6] M. Hudelson, V. Klee, and D. Larman. Largest $j$-simplices in $d$-cubes: some relatives of the Hadamard maximum determinant problem, Linear Algebra Appl., 241–243 (1996), 519–598.
  * [7] I. S. Kudryavcev, E. A. Ozerova, and A. Yu. Ukhalov. New estimates for the norms of minimal projectors, In: Sovremennye Problemy Matematiki i Informatiki. Vypusk 17 (Modern Problems in Mathematics and Informatics. Issue 17), P. G. Demidov Yaroslavl State University, Yaroslavl, 2017, 74–81 (in Russian).
  * [8] P. K. Manjhi and M. K. Rama. Some new examples of circulant partial Hadamard matrices of type $4-H(k\times n)$, Advances and Applications in Mathematical Sciences, 21:5 (2022), 2559–2564.
  * [9] M. V. Nevskii. Estimates for minimal norm of a projector in linear interpolation over vertices of an $n$-dimensional cube, Modeling and Analysis of Information Systems, 10:1 (2003), 9–19.
  * [10] M. V. Nevskii. Minimal projectors and largest simplices, Modeling and Analysis of Information Systems, 14:1 (2007), 3–10.
  * [11] M. V. Nevskii. On a certain relation for the minimal norm of an interpolation projector, Modeling and Analysis of Information Systems, 16:1 (2009), 24–43.
  * [12] M. V. Nevskii. Geometricheskie Ocenki v Polinomialnoi Interpolyacii (Geometric Estimates in Polynomial Interpolation), Yaroslavl’: Yarosl. Gos. Univ., 2012 (in Russian).
  * [13] M. V. Nevskii and A. Yu. Ukhalov. On optimal interpolation by linear functions on an $n$-dimensional cube, Modeling and Analysis of Information Systems, 25:3 (2018), 291–311
  * [14] M. V. Nevskii and A. Yu. Ukhalov. Linear interpolation on a Euclidean ball in ${\mathbb{R}}^{n}$, Modeling and Analysis of Information Systems, 26:2 (2019), 279–296.
  * [15] D. Slepian. The content of some extreme simplices, Pacific J. Math., 31 (1969), 795–808.
  * [16] D. Vandev. A minimal volume ellipsoid around a simplex, C. R. Acad. Bulg. Sci., 45:6 (1992), 37–40.


# Morphology of graphene flakes in Ni-graphene nanocomposites and its
influence on hardness: an atomistic study

Vardan Hoviki Vardanyan  Physics Department, University Kaiserslautern, Erwin-
Schrödinger-Straße, D-67663 Kaiserslautern, Germany  Research Center OPTIMAS,
University Kaiserslautern, Erwin-Schrödinger-Straße, D-67663 Kaiserslautern,
Germany  Herbert M. Urbassek<EMAIL_ADDRESS>http://www.physik.uni-
kl.de/urbassek/ Physics Department, University Kaiserslautern, Erwin-
Schrödinger-Straße, D-67663 Kaiserslautern, Germany  Research Center OPTIMAS,
University Kaiserslautern, Erwin-Schrödinger-Straße, D-67663 Kaiserslautern,
Germany

###### Abstract

The effect of graphene flakes on the strength of Ni-graphene composites is
investigated using molecular dynamics simulation. Rather than introducing
flakes as flat structures into the Ni matrix, as it is common in available
studies, we introduce them into the heated liquid Ni and let the structures
equilibrate in a 14-ns molecular-dynamics run; these structures are then
quenched to obtain the composites. By varying the interaction of flake edge
atoms with the Ni matrix, two different flake morphologies – wrinkled vs flat
– are obtained. The mechanical properties, and in particular the composite
hardness, are investigated by a simulated nanoindentation test. The flake
morphology affects the plastic activity and the hardness of the composites.
Wrinkled flakes show a higher potential in absorbing dislocations than flat
flakes, resulting in a considerably reduced hardness of the composite. No
effect of the changed interaction between graphene edge atoms and the Ni
matrix on the dislocation activity and hardness could be observed.
Furthermore, we demonstrate that a high graphene content in the plastic zone
leads to an increased dislocation absorption and weakens the composite.

nickel-graphene composites; flake morphology; nanoindentation; hardness;
molecular dynamics

##  I Introduction

Graphene-metal composites have excellent mechanical properties Xiong  _et al._
(2015); Ramanathan  _et al._ (2008); Zhang  _et al._ (2014); Yang  _et al._
(2016); Guo  _et al._ (2019); Feng  _et al._ (2017a); Liu  _et al._ (2016) as
they combine the superior mechanical properties of graphene Lee  _et al._
(2008); Cao  _et al._ (2020) with the ductility of the matrix metal. The
increasing interest in graphene-based metal nanocomposites Feng  _et al._
(2017b); Hwang  _et al._ (2013); Khan  _et al._ (2015) requires a better
understanding of the hardening mechanisms in such materials.

The dislocation-interface interaction is key to understanding the mechanical
properties of graphene-metal composites. Experimental studies indicate the
importance of dislocation pile-up and blocking at the metal-graphene interface
as the source of the graphene-induced hardening in the composite Hwang  _et
al._ (2017); Li  _et al._ (2015); Kim  _et al._ (2013).

Molecular dynamics simulations have been proven as an excellent method to
study the dislocation-interface interaction, such as for metal-amorphous
interfaces Zhou and Chen (2016); Sterwerf  _et al._ (2017); Song  _et al._
(2017); Avila  _et al._ (2021), metallic multilayers Wang and Misra (2011) or
poly-crystalline systems Shabib and Miller (2009); Schiøtz  _et al._ (1998);
Van Swygenhoven and Weertman (2006); Zhou  _et al._ (2020). Several molecular
dynamics studies were devoted to investigating graphene-metal nanocomposites
Chang  _et al._ (2013); Muller  _et al._ (2018); Liu  _et al._ (2016); _et
al._ (2017); Vardanyan and Urbassek (2019, 2020); Shuang and Aifantis (2021);
Shuang  _et al._ (2021); Zhang  _et al._ (2020a). The considerable improvement
of the composite hardness Peng and Sun (2020); Rezaei (2018) was explained by
the hindering of dislocation propagation in the metal matrix by graphene Chang
_et al._ (2013); _et al._ (2017); Ma  _et al._ (2020). In addition, the
importance of the absorption of dislocations by the interface was reported
Shuang and Aifantis (2021, 2020) and it was found that dislocation absorption
can weaken the system Vardanyan and Urbassek (2019, 2020).

In further detail, Shuang and Aifantis Shuang and Aifantis (2021) pointed out
that the edges of graphene flakes are important for the composite strength as
they act as dislocation sources in compression experiments. They also argued
that the previously reported improvement of hardening Zhang  _et al._ (2019a);
Rezaei  _et al._ (2016); Zhang  _et al._ (2019b) may be an artifact of the
simulation design, if graphene flakes crossing the simulation cell with
periodic boundary are used. The authors recommend to use free graphene flakes
in future simulations.

Xia _et al._ Xia  _et al._ (2021) demonstrated that the strength of graphene-
metal nanocomposites depends on the metal grain size and the graphene volume
concentration. Wei _et al._ Wei  _et al._ (2021) showed that there is an
optimum concentration of graphene for better mechanical performance such that
any further increase of the graphene content lowers the mechanical resistivity
of the composite.

The previously published research in metal-graphene composites Han  _et al._
(2021); Rong  _et al._ (2018); Zhu  _et al._ (2021); Wang  _et al._ (2021);
Zhang  _et al._ (2020b); Kumar (2018); Shuang and Aifantis (2021); Shuang  _et
al._ (2021); Zhang  _et al._ (2020a, 2021); Vardanyan and Urbassek (2019,
2020) lacks a discussion of the influence of the graphene concentration on the
composite hardness, and in particular on a transition from hardening to
weakening of the composite structures. As a rule, if the simulated composite
shows hardening by graphene reinforcement, a further increase of the graphene
content further increases the hardness Weng  _et al._ (2018); Zhang  _et al._
(2019a); Ma  _et al._ (2020); Shuang  _et al._ (2021); and analogous for a
weakening effect Chang  _et al._ (2013); Zhang  _et al._ (2020b). This
contradicts numerous experimental findings Yousefi  _et al._ (2013); Chu and
Jia (2014); Hu  _et al._ (2016); Zhang and Zhan (2016); Chen  _et al._
(2016a); Gao  _et al._ (2016); Li and Xiong (2017); Chen  _et al._ (2016b)
which indicate that there exists an optimum graphene fraction in the component
to maximize the hardness.

Most previous studies use a flat morphology of graphene flakes. However, it
well known that graphene has a wrinkled structure in composites Gao  _et al._
(2016); Song  _et al._ (2020) and molecular dynamics simulations also show
easy wrinkling of graphene Varillas and Frank (2021). The bending stiffness of
a single graphene layer was reported Wei  _et al._ (2013); Deng and Berry
(2016); Chen  _et al._ (2019); Varillas and Frank (2021) and it was shown how
the difference between in-plane stiffness and out-of-plane bending stiffness
allows crumpling of graphene flakes and the formation of stretchable
structures Blees  _et al._ (2015); Cao  _et al._ (2020). Moreover,
experimental studies show that graphene crumpling has a significant effect on
softening a free-standing graphene flake Nicholl  _et al._ (2015). Recently,
Kumar _et al._ Kumar  _et al._ (2019) studied aggregation of graphene flakes
in a liquid aluminum matrix and showed a way for simulations to obtain more
complex flake structures than the usually used flat shape.

In the present study, we focus on the effect of the morphology of graphene
flakes on the composite hardness. By modifying the graphene-Ni interaction,
graphene flakes change their morphology in the liquid Ni melt and keep this
morphology during quenching into the solid state. Simulated nanoindentation of
these nanocomposite structures shows that the morphology strongly influences
the composite hardness. In addition, we evaluate the effect of the local
graphene content in the plastic zone under the indenter, and of pre-existing
dislocations on the hardening and weakening of the composite with respect to a
polycrystalline Ni sample.

##  II Method

###  II.1 Preparation of the systems

The system preparation starts with a liquid block of nickel at 2300 K with
dimensions of $25.6\times 25.6\times 21.3$ nm, containing approximately
$1.1\times 10^{6}$ atoms. 64 Graphene flakes are evenly distributed in the
liquid, see Fig. 1. Each flake has a square shape with edge length 26 Å
containing 288 carbon atoms. Thus the atomic fraction of graphene in the
composite amounts to $1.6$ %.

Figure 1: Preparation of the Ni-graphene nanocomposite system. The cuboid
simulation volume is filled with liquid Ni atoms (yellow) and graphene flakes
(blue); graphene edge atoms with a modified C-Ni interaction are colored red.
Half of the nickel atoms are removed for clarity. The inset shows
schematically the difference in our models of graphene: at the top graphene
with strongly bound edge atoms; at the bottom a graphene flake with unmodified
edge atoms see text.

In our simulation the interaction between Ni atoms is modeled by the
interaction potential developed by Mishin _et al._ Mishin  _et al._ (1999),
and the interaction between carbon atoms by the bond-order potential of Stuart
_et al._ Stuart  _et al._ (2000). The Ni-C interaction is modeled by a
Lennard-Jones (LJ) potential. Here, we consider two choices:

  1. (i)

All C atoms interact via the same LJ potential. We choose the parameters as
given by Huang _et al._ Huang  _et al._ (2003): $\varepsilon=23.049$ meV,
$\sigma=2.852$ Å. In the following, we will denote this interaction and the
embedded graphene flakes as _graphene with unmodified edges_.

  2. (ii)

It is known that edge atoms of the graphene flakes can interact covalently
with Ni. Tavazza _et al._ Tavazza  _et al._ (2015) calculated the Ni-C
interaction by DFT and fitted their results to a LJ potential with the
parameters $\varepsilon=200$ meV, $\sigma=1.514$ Å. We will denote graphene
flakes as _graphene with strongly bound edge atoms_ , in which edge atoms
underlie this increased interaction, while the inner graphene atoms feel the
weak LJ potential.

The difference of our models of weak and strong graphene is sketched in the
inset of Fig. 1.

We note that in previous studies of the Ni-graphene system Kumar  _et al._
(2019); Vardanyan and Urbassek (2019, 2020), the effects of a changed
interaction potential between C and Ni atoms were studied, but it was changed
for all C atoms, including those in the interior of the flakes. Here for the
first time, we differentiate between interior flake atoms and edge atoms, and
investigate to what extent their different interaction with Ni influences the
composite preparation and the mechanical properties of the composite.

Besides the Ni-graphene composites, we also prepared a pure Ni system with the
same size and following the same preparation procedure for comparison
purposes.

###  II.2 Equilibration

We let the system equilibrate in an NPT ensemble at 2300 K for 14 ns. This
time is sufficient to equilibrate the systems, see Supplementary Material (SM)
NiG . During the equilibration process, the graphene flakes changed their
positions and even their shapes in the liquid metal. Because of the strong
attraction of C atoms modeled by the AIREBO potential Stuart  _et al._ (2000),
graphene flakes tend to aggregate and form covalent bonds between each other.
However, distinct differences show up between graphene with strongly bound
edge atoms and unmodified edges, see Fig. 2.

Figure 2: Graphene structures obtained after equilibration at 2300 K for (a)
graphene flakes with strongly bound edge atoms and (b) with unmodified edges.
The red dashed lines show a cut through the substrate as displayed in the
insets. Color code as in Fig. 1.  Figure 3: Crystallographic structure of Ni-
graphene composites after quenching for (left) graphene flakes with strongly
bound edge atoms and (right) with unmodified edges. Atoms are colored
according to common-neighbor analysis. Green: fcc nickel; dark blue: SFs in
nickel; yellow: grain boundaries; light blue: graphene; red: edge atoms of
graphene flakes. The yellow arrow shows the propagation of partial
dislocations leaving behind a twin structure.

In the case of graphene with strongly bound edge atoms, Fig. 2(a), the flakes
are flat and contact each other at their edges, since the AIREBO potential
describing the C-C interaction shows an increased attraction between C atoms
there. The final system shows a flat morphology, since the increased C-Ni
attraction at the edges prevents further flakes from approaching.

As we show in the SM, NiG , graphene flakes with unmodified edges tend to roll
up in the liquid Ni metal; the mutual attraction of edge atoms stabilizes the
tubular structures. When other flakes approach such a tube, it may attach to
the edges resulting in the wrinkled structures we observe. If the Ni-C
attraction is increased in strongly bound edge atoms, the flake edges are
surrounded by a quite stable Ni shell; this Ni shell reduces the amount of Ni
flakes rolling up while not completely erasing it. In addition, it may be
supposed that the Ni shell enhances the flexural stiffness of the flakes and
thus impedes the roll-up process. As a consequence, we observe a tendency
towards flat structures in flakes with strongly bound edge atoms.

Also in graphene with unmodified edges, the C-C interaction induces the flakes
to approach each other at the edges. However, further structural changes occur
that are driven by reactive interactions of the edges of graphene flakes. The
flakes continue further aggregation resulting in a wrinkled morphology with
rough surfaces and encapsulated voids.

In experiments, it is known that the structure of graphene flakes inside
metal-graphene composites is wrinkled rather than flat Yousefi  _et al._
(2013); Chu and Jia (2014); Hu  _et al._ (2016); Zhang and Zhan (2016); Chen
_et al._ (2016a); Gao  _et al._ (2016); Li and Xiong (2017); Chen  _et al._
(2016b); Song  _et al._ (2020).

We conclude that by modifying the interaction of graphene edge atoms with Ni,
it is possible to modify the morphology of the graphene flakes within the
liquid Ni metal. In the following we will discern the _flat morphology_ from
the _wrinkled morphology_ of composites.

Figure 4: Creation of 4 Ni-graphene composites: Two systems with different Ni-
graphene interaction were equilibrated at 2300 K, resulting in a flat and
wrinkled flake morphology, respectively. After quenching the composite systems
were subjected to graphene interactions with strongly bound edge atoms and
with unmodified edges and then relaxed, such that 4 different composite
systems were created. Color code as in Fig. 1.

###  II.3 Quenching

After equilibration of the liquid, the systems are quenched by decreasing the
temperature from 2300 K to 0.1 K with a rate of 1 K/ps. Due to the high
quenching rate, Ni does not crystallize as a single crystal, but in a poly-
crystalline form, as was discussed in detail for the case of Al in Ref. Hou
_et al._ (2016). We observe an average grain size of 4.2 (4.8) nm for the
composite with wrinkled (flat) morphology. For pure Ni, which underwent the
same quenching procedure, the grain size is comparable, 4.1 nm. In both
composites, the flake conformation – wrinkled and flat, respectively – was
conserved during the quenching process, see Fig. 3, and the graphene flakes
are positioned inside the grain boundaries. In the case of wrinkled graphene,
encapsulated voids can be found within the graphene, and all quenched systems
show a high amount of defects – in particular stacking fault planes – in the
grains that originated during the quenching process.

In order to clearly differentiate between the effects of graphene morphology
and C-Ni interaction, we create two further systems from these quenched
structures by changing the Ni-C interaction of edge atoms and then relaxing
the systems.

We thus have the following composite systems:

  1. (i)

A composite with flat morphology with strongly bound edge atoms.

  2. (ii)

A composite with flat morphology with unmodified edges. It is created from
system (i) by changing the interaction of all edge atoms.

  3. (iii)

A composite with wrinkled morphology with unmodified edges.

  4. (iv)

A composite with wrinkled morphology with strongly bound edge atoms. It is
created from system (iii) by modifying the Ni-C interaction of edge atoms.

These four composites allow us to differentiate the effects of graphene
morphology and Ni-C interaction on mechanical properties. The generation of
these composite systems and their morphologies are summarized in Fig. 4.

Figure 5: (a) Evolution of the total dislocation length with indentation
depth. Data are for a representative simulation for each system. (b)
Dependence of the length of indentation-induced dislocations on the initial
dislocation length. Each data point is the result of an individual indentation
simulation.

###  II.4 Creation of surfaces

During equilibration and quenching we apply periodic boundary conditions in
all cartesian directions. In order to perform indentation simulations, a free
surface is created. For increasing statistics, we create 5 replicas of each
quenched system and introduce a free surface in each system at a different
position.

After introducing free boundary conditions at these surfaces, the
nanocomposite systems are relaxed in the isothermal-isobaric NPT ensemble.

###  II.5 Indentation simulations

Indentation simulations are performed by modeling the indenter by a repulsive
potential Kelchner  _et al._ (1998),

$V(r)=\left\\{\begin{array}[]{ll}k(R-r)^{3},&r<R,\\\ 0,&r\geq
R.\end{array}\right.$ (1)

The indenter radius is $R=40$ Å and the indenter stiffness $k=10$ eVÅ-3
Kelchner  _et al._ (1998); Ziegenhain  _et al._ (2009). Using a displacement-
controlled algorithm Ruestes  _et al._ (2017), the indenter is moved
perpendicular into the surface to a final depth of 40 Å with a velocity of 20
m/s.

For each system replica, we perform 5 independent indentation simulations
which differ from each other by the exact lateral positioning of the indenter;
it was moved randomly to another position by around $\pm$ 5 Å. Thus we have 25
indentation simulations for each system.

The simulations are performed with the open-source code LAMMPS Plimpton (1995)
using a constant time step of 1 fs. Common-neighbor analysis (CNA) Stukowski
(2012) is used to identify the crystalline structure and the dislocation
detection algorithm Stukowski (2012); Stukowski  _et al._ (2012); Stukowski
and Arsenlis (2012) to determine the length of dislocations. OVITO Stukowski
(2010) is used to visualize the simulation results.

##  III Results

###  III.1 Dislocation generation

Figure 6:  Side view of the wrinkled composite with unmodified edges after
indentation. (a) VMSS map. White dashed line shows the boundary of the plastic
zone. Graphene atoms are colored light blue. Ni atoms are colored according to
their VMSS. (b) Crystallographic structure of this system, where atoms are
colored according to common-neighbor analysis. Color code as in Fig. 3.

Figure 7:  a) Top view onto the $x$–$y$ plane and b) side view onto the
$x$–$z$ of the plastic zone in the composite with flat morphology and strongly
bound edges. c) Top view onto the $x$–$y$ plane and d) side view onto the
$x$–$z$ of the plastic zone in the composite with wrinkled morphology and
unmodified edges. Dislocations are colored according to their Burgers vector.
Green: $\frac{1}{2}\langle 112\rangle$; dark blue: $\frac{1}{2}\langle
110\rangle$; pink: $\frac{1}{6}\langle 110\rangle$; yellow:
$\frac{1}{3}\langle 001\rangle$; red: other. Yellow shading shows indented
surface. Gray contours highlight Ni grains.

Fig. 5(a) shows the evolution of the total dislocation length during
indentation for representative cases. In the following, we consider only
dislocations generated inside the grains and leave out those inside the grain
boundaries or very close to them. As a result of the solidification of
composites, all systems already have a dislocation network before indentation,
which we denote as ‘initial dislocations’. During indentation, new
dislocations are nucleated mainly from the grain boundaries as well as from
high-stress regions in the grains close to the indenter; this can be seen in
Fig. 5(a) as a general steady increase of the total dislocation length with
increasing indentation depth. This increase is dominant in the case of the
flat graphene morphology, while the wrinkled morphology also shows strong
reductions in the dislocation length. Unlike in the case of single-crystalline
metals, where the total dislocation length exhibits a characteristic steady
increase Vardanyan and Urbassek (2019), here the absorption of dislocations by
grain boundaries and graphene flakes acts to reduce the dislocation length.
Fig. 5(a) thus shows that the wrinkled morphology is more effective in
absorbing dislocations.

Fig. 5(b) shows the dependence of the length of generated dislocations on the
initial dislocation length for all 100 simulations. Due to the various
surfaces cut out of the quenched composite system and the inaccuracies of the
dislocation detection algorithm Stukowski (2012); Stukowski  _et al._ (2012);
Stukowski and Arsenlis (2012), the length of the initial dislocations shows
some spread. The spread in the generated dislocations is even larger, since
for each generated surface, 5 random indentation points were chosen and the
generated plasticity depends on the exact indentation point. Here and in the
following, we determine the length of generated dislocations by the difference
between the final and initial length of dislocations for each system. Note
that the length of generated dislocations does not correlate with the initial
dislocation length. On average, systems with a flat flake morphology show an
increase of the dislocation network by on average 80 nm, while in the case of
the wrinkled morphology the increase is smaller, 58 nm. In selected cases, a
wrinkled geometry may even lead to a reduction of the dislocation network due
to the strong absorption activity of the wrinkled flakes. In no cases, an
effect of the C-Ni interaction of the edge atoms on the length of the
generated dislocations could be observed.

The plastic deformation is mapped by the von-Mises shear strain (VMSS) in Fig.
6 (a). The plastic deformation of a polycrystalline sample is a consequence of
dislocation nucleation inside the grains and shear activation in grain
boundaries. We identify the boundary of the plastic zone in this figure such
that all shear strains are situated in the plastic zone. We find that in all
composites studied, the radius of the plastic zone amounts to roughly 10–15
nm. Fig. 6 (b) displays a CNA analysis of the same system as in Fig. 6 (a). It
shows that the two wrinkled graphene flakes existing in the plastic zone
create large defective regions inside the plastic zone.

The dislocation structure in the grains after indentation is shown in Fig. 7
for two different composites.

  * •

Fig. 7 (a,b) shows the dislocation structure in the case of the flat
morphology with strongly bound edge atoms. In this particular simulation,
there are 2690 C atoms inside the plastic zone. The total length of
dislocations amounts to 1290 Å. This large value indicates pile-up of
dislocations inside the plastic zone. Dislocations can be seen around the
graphene flakes. Thus, the graphene-nickel interface in case of flat
morphology is an additional barrier for dislocations.

  * •

Fig. 7 (c,d) shows the dislocation structure in the case of the wrinkled
morphology with unmodified edges. There are 6301 C atoms inside the plastic
zone and the total length of dislocations amounts to 540 Å. This morphology
results in a more inhomogeneous dispersion of graphene inside the plastic
zone. Because graphene is not flat, the interface between graphene and nickel
has a higher tendency of absorbing dislocations. Dislocation nucleation occurs
only in grains far from graphene.

Fig. 8 provides quantitative details about the dependence of the dislocation
generation and the distribution of graphene in the plastic zone and around it.
To this end we determine the length of dislocations before and after
indentation in spherical shells of width 1 nm around the indentation point.
Two scenarios have been selected. Fig. 8a shows a case, where the amount of
graphene within the plastic zone is small, but more massive outside, at a
distance of 15 nm from the indentation point. As a consequence, a considerable
increase of the dislocation length in the plastic zone and a pile-up of
dislocations in front of the graphene peak at 15 nm can be observed. In the
second case, Fig. 8b, the amount of graphene within the plastic zone is high,
and not so massive outside of the zone. Here the length of dislocations within
the plastic zone is reduced during indentation; this plot directly shows the
effectivity of graphene to absorb dislocations.

We conclude that the dislocations that are invariably created during the
indentation are absorbed by the graphene flakes; a large amount of graphene in
the plastic zone may even reduce the initial length of dislocations. The
effect of the flake morphology shows up by a higher absorption activity of
wrinkled flakes compared to flat flakes. The effect of C-Ni interaction on
dislocation activity is negligible.

###  III.2 Hardness

Figure 8:  Radial dependence of the initial and final dislocation length as
well as the number of carbon atoms for (a) the flat morphology with strongly
bound edges and (b) the wrinkled morphology with unmodified edges. Data show
quantities in radial shells of 1 nm width around the indentation point. The
vertical dashed line indicates the extension of the plastic zone.

Hardness was determined as the contact pressure averaged over indentation
depths between 20 Å and 40 Å, where the hardness already showed its stable,
saturated value, see SM NiG . The results of all 100 indentations are shown in
Fig. 9. Fig. 9 (a) shows the dependence of the hardness on the length of
generated dislocations. Even though the data show a considerable scatter,
several tends can be observed. Thus, a common trend is seen in that all
systems show a steady increase of hardness upon increasing of dislocation
length. This is caused by dislocation pileup which increases hardness Misra
_et al._ (2005); Wang and Misra (2011); Beyerlein  _et al._ (2015).

Polycrystalline nickel – with a grain size of 4.1 nm – shows a hardness of
9.95 GPa. This is in reasonable agreement with literature data Liu  _et al._
(2013) that obtained around 12 GPa for 5-nm grained polycrystalline Ni; the
slight difference may be explained by the different interatomic interaction
potentials used.

Composites with flat morphology show an enhanced hardness as compared to the
wrinkled morphology, see Fig. 9a; the difference amounts to around 1 GPa on
average. The hardening of polycrystalline Ni by flat flakes is particularly
well seen by comparing the hardness results with those of polycrystalline Ni:
the hardness of teh flat morphology is in all cases above that of
polycrystalline nickel. This means a definite strengthening of the composite.
The reason lies in the way how the 2D graphene flake can sustain the stress
field exerted by the indenter on the composite material. In the case of the
flat morphology, in-plane deformation of the flake will result in hardening,
while in the case of the wrinkled morphology, in addition out-of-plane
deformations will assist to sustain the load. As the stiffness of 2D materials
is considerably stronger for in-plane stresses than for out-of-plane stresses,
the flat material will be more effective in sustaining these stresses and thus
enable a hardening of the composite as seen in our Fig. 9a.

In general, wrinkled graphene shows less resistance against deformation than
flat graphene. Nicholl _et al._ Nicholl  _et al._ (2015) find that the in-
plane stiffness of graphene is strongly reduced upon crumpling. In our
simulations, we find that wrinkled graphene also more easily deforms during
loading of the composite resulting in out-of-plane deformations and thus
offering an additional mode of reducing the hardness of the composite (see SM,
NiG , for an example). In addition, wrinkled graphene allows for a similar or
even better way to dislocation absorption than flat graphene.

The effect of flat graphene on the mechanical properties of Ni has been
repeatedly studied in the past Chang  _et al._ (2013); Muller  _et al._
(2018); Liu  _et al._ (2016); _et al._ (2017); Vardanyan and Urbassek (2019,
2020); Shuang and Aifantis (2021); Shuang  _et al._ (2021); Zhang  _et al._
(2020a). There it has been shown that graphene mainly absorbs dislocations and
hence weakens the nickel matrix as compared to a single-crystalline Ni block
whose hardness amounts to 15 GPa Vardanyan and Urbassek (2019). Yazdandoost
_et al._  _et al._ (2017) investigated the effect of graphene flakes inserted
into the grain boundaries of a polycrystalline Ni matrix and found almost no
effect on hardness, arguing that dislocation nucleation and blocking cancel in
the case of a random orientation of the flat graphene flakes. In the present
work, we find that flat graphene flakes embedded in the Ni grain boundaries
harden the polycrystalline Ni matrix. We presume that this occurs because our
Ni grain boundaries and the enclosed flakes have been carefully relaxed in our
sample preparation, Sect. II, in contrast to the procedure based on Voronoi
tesselation in Ref.  _et al._ (2017). In the case of wrinkled flakes, besides
the obvious effect of dislocation absorption which weakens the matrix, the
out-of-plane deformation of the flakes adds a second effect such that the
overall hardness of the wrinkled composite lies below that of polycrystalline
Ni. The present work shows that the morphology of graphene affects the
hardness of the composite in a more refined manner, such that both weakening
and softening of the matrix material becomes possible.

Fig. 9 (b) shows the dependence of the composite hardness on the number of
carbon atoms in the plastic zone. In all cases, the hardness decreases with
increasing graphene content. Thus, graphene flakes are not per se
strengthening the composite. Rather, graphene flakes act as defect structures
that absorb dislocations.

Fig. 9 (c) shows the effect of graphene atoms on dislocation formation inside
the plastic zone. An increase of the graphene concentration leads on average
to less dislocation generation in the plastic zone; this trend is visible for
all systems studied. It is caused by the absorption of dislocations by
graphene flakes in the plastic zone.

In all our correlations, Figs. 9a–c, no difference in hardness based on the
Ni-C edge atom interaction shows up. In all cases, flakes with strongly bound
edges and unmodified edge atoms feature on average the same hardness. We thus
conclude that while the effect of the graphene morphology on the mechanical
properties is strong, the effect of the C-Ni interaction of the edge atoms
itself is minor.

Figure 9: (a) Dependence of composite hardness on (a) the length of
dislocations generated in grains and (b) the number of carbon atoms in the
plastic zone. (c) Dependence of the length of dislocations generated in grains
on the number of carbon atoms in the plastic zone. Lines are to guide the eye.
The gray line gives the average hardness of polycrystalline Ni.

##  IV Summary

Using molecular dynamics simulation, we studied the effect of the morphology
of graphene flakes in a Ni-graphene nanocomposite on the mechanical
properties, in particular the hardness. The flake morphology was established
by the self-organization of graphene flakes in a liquid Ni matrix. By varying
the interaction potential of graphene edge atoms with Ni, two different flake
morphologies – flat vs wrinkled – could be generated. The modification of the
interaction potential is justified by the fact that edge atoms have
unsaturated bonds, in contrast to interior flake atoms, and are therefore
prone to generating strong Ni-C bonds. In the quenched solid structures of the
Ni-graphene nanocomposites, graphene is located in the grain boundaries of Ni.

The mechanical properties of the composite were evaluated in a simulated
nanoindentation test. As known from previous studies Misra  _et al._ (2005);
Wang and Misra (2011); Beyerlein  _et al._ (2015), the dislocations generated
by the indentation process will pile up in front of the Ni grain boundaries
and the graphene flakes, leading to a hardening of the composite compared to a
pure poly-crystalline Ni specimen. We could verify by selected case studies
that a high concentration of graphene in the plastic zone of the nanoindenter
leads to a strong absorption of dislocations, such that even some of the
dislocations present by the solidification process before the indentation were
annihilated. However, this effect only occurs for the wrinkled flake
morphology; there, the absorption of dislocations at the flake dominates the
pile-up and leads to a strong reduction in hardness, amounting to roughly10 %
in the case investigated.

While the effect of the flake morphology on the hardness is strong, the
changed interaction between graphene edge atoms and the Ni matrix shows only
negligible effects on the dislocation activity and hardness. Our results thus
suggest that future studies of mechanical properties of metal-graphene
composites need to take the flake morphology in due account.

###### Acknowledgements.

We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) – project number 252408385 – IRTG 2057. Simulations were
performed at the High Performance Cluster Elwetritsch (RHRK, TU
Kaiserslautern, Germany).

## References

  * Xiong  _et al._ (2015) D.-B. Xiong, M. Cao, Q. Guo, Z. Tan, G. Fan, Z. Li, and D. Zhang, ACS Nano 9, 6934 (2015).
  * Ramanathan  _et al._ (2008) T. Ramanathan, A. A. Abdala, S. Stankovich, D. A. Dikin, M. Herrera-Alonso, R. D. Piner, D. H. Adamson, H. C. Schniepp, X. Chen, R. S. Ruoff, S. T. Nguyen, I. A. Aksay, R. K. Prud’Homme, and L. C. Brinson, Nature nanotechnology 3, 327 (2008).
  * Zhang  _et al._ (2014) P. Zhang, L. Ma, F. Fan, Z. Zeng, C. Peng, P. E. Loya, Z. Liu, Y. Gong, J. Zhang, X. Zhang, P. M. Ajayan, T. Zhu, and J. Lou, Nature Communications 5, 3782 EP (2014).
  * Yang  _et al._ (2016) Z. Yang, D. Wang, Z. Lu, and W. Hu, Applied Physics Letters 109, 191909 (2016).
  * Guo  _et al._ (2019) Q. Guo, K. Kondoh, and S. M. Han, MRS Bulletin 44, 40 (2019).
  * Feng  _et al._ (2017a) Q. Feng, X. Song, H. Xie, H. Wang, X. Liu, and F. Yin, Materials & Design 120, 193 (2017a).
  * Liu  _et al._ (2016) X. Liu, F. Wang, W. Wang, and H. Wu, Carbon 107, 680 (2016).
  * Lee  _et al._ (2008) C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008).
  * Cao  _et al._ (2020) K. Cao, S. Feng, Y. Han, L. Gao, T. Hue Ly, Z. Xu, and Y. Lu, Nature Communications 11, 284 (2020).
  * Feng  _et al._ (2017b) S. Feng, Q. Guo, Z. Li, G. Fan, Z. Li, D.-B. Xiong, Y. Su, Z. Tan, J. Zhang, and D. Zhang, Acta Mater. 125, 98 (2017b).
  * Hwang  _et al._ (2013) J. Hwang, T. Yoon, S. H. Jin, J. Lee, T.-S. Kim, S. H. Hong, and S. Jeon, Advanced Materials 25, 6724 (2013).
  * Khan  _et al._ (2015) M. Khan, M. N. Tahir, S. F. Adil, H. U. Khan, M. R. H. Siddiqui, A. A. Al-warthan, and W. Tremel, J. Mater. Chem. A 3, 18753 (2015).
  * Hwang  _et al._ (2017) B. Hwang, W. Kim, J. Kim, S. Lee, S. Lim, S. Kim, S. H. Oh, S. Ryu, and S. M. Han, Nano Letters 17, 4740 (2017).
  * Li  _et al._ (2015) Z. Li, Q. Guo, Z. Li, G. Fan, D.-B. Xiong, Y. Su, J. Zhang, and D. Zhang, Nano Lett. 15, 8077 (2015).
  * Kim  _et al._ (2013) Y. Kim, J. Lee, M. S. Yeom, J. W. Shin, H. Kim, Y. Cui, J. W. Kysar, J. Hone, Y. Jung, S. Jeon, and S. M. Han, Nat. Commun. 4, 2114 (2013).
  * Zhou and Chen (2016) X. Zhou and C. Chen, International Journal of Plasticity 80, 75 (2016).
  * Sterwerf  _et al._ (2017) C. Sterwerf, T. Kaub, C. Deng, G. B. Thompson, and L. Li, Thin Solid Films 626, 184 (2017).
  * Song  _et al._ (2017) H. Y. Song, J. J. Xu, Y. G. Zhang, S. Li, D. H. Wang, and Y. L. Li, Materials & Design 127, 173 (2017).
  * Avila  _et al._ (2021) K. E. Avila, V. H. Vardanyan, S. Küchemann, and H. M. Urbassek, Applied Surface Science 551, 149285 (2021).
  * Wang and Misra (2011) J. Wang and A. Misra, Current Opinion in Solid State and Materials Science 15, 20 (2011).
  * Shabib and Miller (2009) I. Shabib and R. E. Miller, Acta Materialia 57, 4364 (2009).
  * Schiøtz  _et al._ (1998) J. Schiøtz, F. D. Di Tolla, and K. W. Jacobsen, Nature 391, 561 (1998).
  * Van Swygenhoven and Weertman (2006) H. Van Swygenhoven and J. R. Weertman, Materials Today 9, 24 (2006).
  * Zhou  _et al._ (2020) X. Zhou, Z. Feng, L. Zhu, J. Xu, L. Miyagi, H. Dong, H. Sheng, Y. Wang, Q. Li, Y. Ma, H. Zhang, J. Yan, N. Tamura, M. Kunz, K. Lutker, T. Huang, D. A. Hughes, X. Huang, and B. Chen, Nature 579, 67 (2020).
  * Chang  _et al._ (2013) S.-W. Chang, A. K. Nair, and M. J. Buehler, Philosophical Magazine Letters 93, 196 (2013).
  * Muller  _et al._ (2018) S. E. Muller, R. R. Santhapuram, and A. K. Nair, Computational Materials Science 152, 341 (2018).
  *  _et al._ (2017) F. , A. Yari Boroujeni, and R. Mirzaeifar, Phys. Rev. Materials 1, 076001 (2017).
  * Vardanyan and Urbassek (2019) V. H. Vardanyan and H. M. Urbassek, Computational Materials Science 170, 109158 (2019).
  * Vardanyan and Urbassek (2020) V. H. Vardanyan and H. M. Urbassek, Materials 13, 1683 (2020).
  * Shuang and Aifantis (2021) F. Shuang and K. E. Aifantis, Carbon 172, 50 (2021).
  * Shuang  _et al._ (2021) F. Shuang, Z. Dai, and K. E. Aifantis, ACS Applied Materials & Interfaces 13, 26610 (2021).
  * Zhang  _et al._ (2020a) S. Zhang, P. Huang, and F. Wang, Materials & Design 190, 108555 (2020a).
  * Peng and Sun (2020) W. Peng and K. Sun, Mechanics of Materials 141, 103270 (2020).
  * Rezaei (2018) R. Rezaei, Computational Materials Science 151, 181 (2018).
  * Ma  _et al._ (2020) Y. Ma, S. Zhang, Y. Xu, X. Liu, and S.-N. Luo, Phys. Chem. Chem. Phys. 22, 4741 (2020).
  * Shuang and Aifantis (2020) F. Shuang and K. E. Aifantis, Scripta Materialia 181, 70 (2020).
  * Zhang  _et al._ (2019a) C. Zhang, C. Lu, L. Pei, J. Li, and R. Wang, Composites Part B: Engineering , 107610 (2019a).
  * Rezaei  _et al._ (2016) R. Rezaei, C. Deng, H. Tavakoli-Anbaran, and M. Shariati, Philosophical Magazine Letters 96, 322 (2016).
  * Zhang  _et al._ (2019b) C. Zhang, C. Lu, L. Pei, J. Li, R. Wang, and K. Tieu, Carbon 143, 125 (2019b).
  * Xia  _et al._ (2021) X. Xia, Z. Du, J. Zhang, J. Li, and G. J. Weng, International Journal of Engineering Science 162, 103476 (2021).
  * Wei  _et al._ (2021) C. Wei, N. Ye, L. Hong, J. Yao, W. Xia, J. Mao, Y. Wang, Y. Zhao, and J. Tang, ACS Applied Materials & Interfaces 13, 21714 (2021).
  * Han  _et al._ (2021) R. Q. Han, H. Y. Song, and M. R. An, Computational Materials Science 197, 110604 (2021).
  * Rong  _et al._ (2018) Y. Rong, H. P. He, L. Zhang, N. Li, and Y. C. Zhu, Computational Materials Science 153, 48 (2018).
  * Zhu  _et al._ (2021) J. Zhu, X. Liu, Z. Wang, and Q. Yang, Modelling and Simulation in Materials Science and Engineering 29, 055017 (2021).
  * Wang  _et al._ (2021) L. Wang, J. Jin, P. Yang, S. Li, S. Tang, Y. Zong, and Q. Peng, Computational Materials Science 191, 110309 (2021).
  * Zhang  _et al._ (2020b) Y. Zhang, Q. An, J. Li, B. Lu, W. Wu, and R. Xia, Journal of Molecular Modeling 26, 335 (2020b).
  * Kumar (2018) S. Kumar, Materials Chemistry and Physics 208, 41 (2018).
  * Zhang  _et al._ (2021) S. Zhang, F. Wang, and P. Huang, Journal of Materials Science & Technology 87, 176 (2021).
  * Weng  _et al._ (2018) S. Weng, H. Ning, T. Fu, N. Hu, Y. Zhao, C. Huang, and X. Peng, Scientific Reports 8, 3089 (2018).
  * Yousefi  _et al._ (2013) N. Yousefi, X. Lin, Q. Zheng, X. Shen, J. R. Pothnis, J. Jia, E. Zussman, and J.-K. Kim, Carbon 59, 406 (2013).
  * Chu and Jia (2014) K. Chu and C. Jia, physica status solidi (a) 211, 184 (2014).
  * Hu  _et al._ (2016) Z. Hu, G. Tong, D. Lin, Q. Nian, J. Shao, Y. Hu, M. Saeib, S. Jin, and G. J. Cheng, Journal of Materials Processing Technology 231, 143 (2016).
  * Zhang and Zhan (2016) D. Zhang and Z. Zhan, Journal of Alloys and Compounds 654, 226 (2016).
  * Chen  _et al._ (2016a) F. Chen, J. Ying, Y. Wang, S. Du, Z. Liu, and Q. Huang, Carbon 96, 836 (2016a).
  * Gao  _et al._ (2016) X. Gao, H. Yue, E. Guo, H. Zhang, X. Lin, L. Yao, and B. Wang, Powder Technology 301, 601 (2016).
  * Li and Xiong (2017) G. Li and B. Xiong, Journal of Alloys and Compounds 697, 31 (2017).
  * Chen  _et al._ (2016b) Y. Chen, X. Zhang, E. Liu, C. He, C. Shi, J. Li, P. Nash, and N. Zhao, Scientific Reports 6, 19363 (2016b).
  * Song  _et al._ (2020) G. Song, Q. Wang, L. Sun, S. Li, Y. Sun, Q. Fu, and C. Pan, Materials & Design 191, 108629 (2020).
  * Varillas and Frank (2021) J. Varillas and O. Frank, Carbon 173, 301 (2021).
  * Wei  _et al._ (2013) Y. Wei, B. Wang, J. Wu, R. Yang, and M. L. Dunn, Nano Letters 13, 26 (2013).
  * Deng and Berry (2016) S. Deng and V. Berry, Materials Today 19, 197 (2016).
  * Chen  _et al._ (2019) W. Chen, X. Gui, L. Yang, H. Zhu, and Z. Tang, Nanoscale Horiz. 4, 291 (2019).
  * Blees  _et al._ (2015) M. K. Blees, A. W. Barnard, P. A. Rose, S. P. Roberts, K. L. McGill, P. Y. Huang, A. R. Ruyack, J. W. Kevek, B. Kobrin, D. A. Muller, and P. L. McEuen, Nature 524, 204 (2015).
  * Nicholl  _et al._ (2015) R. J. T. Nicholl, H. J. Conley, N. V. Lavrik, I. Vlassiouk, Y. S. Puzyrev, V. P. Sreenivas, S. T. Pantelides, and K. I. Bolotin, Nature Communications 6, 8789 (2015).
  * Kumar  _et al._ (2019) S. Kumar, S. K. Pattanayek, and S. K. Das, The Journal of Physical Chemistry C 123, 18017 (2019).
  * Mishin  _et al._ (1999) Y. Mishin, D. Farkas, M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. B 59, 3393 (1999).
  * Stuart  _et al._ (2000) S. J. Stuart, A. B. Tutein, and J. A. Harrison, J. Chem. Phys. 112, 6472 (2000).
  * Huang  _et al._ (2003) S.-P. Huang, D. S. Mainardi, and P. B. Balbuena, Surface Science 545, 163 (2003).
  * Tavazza  _et al._ (2015) F. Tavazza, T. P. Senftle, C. Zou, C. A. Becker, and A. C. T. van Duin, J. Phys. Chem. C 119, 13580 (2015).
  * (70) See Supplementary Material for additional data and analyses.
  * Hou  _et al._ (2016) Z. Y. Hou, K. J. Dong, Z. A. Tian, R. S. Liu, Z. Wang, and J. G. Wang, Phys. Chem. Chem. Phys. 18, 17461 (2016).
  * Kelchner  _et al._ (1998) C. L. Kelchner, S. J. Plimpton, and J. C. Hamilton, Phys. Rev. B 58, 11085 (1998).
  * Ziegenhain  _et al._ (2009) G. Ziegenhain, A. Hartmaier, and H. M. Urbassek, J. Mech. Phys. Sol. 57, 1514 (2009).
  * Ruestes  _et al._ (2017) C. J. Ruestes, E. M. Bringa, Y. Gao, and H. M. Urbassek, in  _Applied Nanoindentation in Advanced Materials_, edited by A. Tiwari and S. Natarajan (Wiley, Chichester, UK, 2017) Chap. 14, pp. 313–345.
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995),  _http://lammps.sandia.gov/_.
  * Stukowski (2012) A. Stukowski, Model. Simul. Mater. Sci. Eng. 20, 045021 (2012).
  * Stukowski  _et al._ (2012) A. Stukowski, V. V. Bulatov, and A. Arsenlis, Model. Simul. Mater. Sci. Eng. 20, 085007 (2012).
  * Stukowski and Arsenlis (2012) A. Stukowski and A. Arsenlis, Model. Simul. Mater. Sci. Eng. 20, 035012 (2012).
  * Stukowski (2010) A. Stukowski, Model. Simul. Mater. Sci. Eng. 18, 015012 (2010),  _http://www.ovito.org/_.
  * Misra  _et al._ (2005) A. Misra, J. P. Hirth, and R. G. Hoagland, Acta Mater. 53, 4817 (2005).
  * Beyerlein  _et al._ (2015) I. J. Beyerlein, M. J. Demkowicz, A. Misra, and B. P. Uberuaga, Progress in Materials Science 74, 125 (2015).
  * Liu  _et al._ (2013) X. Liu, F. Yuan, and Y. Wei, Applied Surface Science 279, 159 (2013).


# Travelling waves in a minimal go-or-grow model of cell invasion

Carles Falcó<EMAIL_ADDRESS>Rebecca M. Crossley
<EMAIL_ADDRESS>Ruth E. Baker<EMAIL_ADDRESS>
###### Abstract

We consider a minimal go-or-grow model of cell invasion, whereby cells can
either proliferate, following logistic growth, or move, via linear diffusion,
and phenotypic switching between these two states is density-dependent. Formal
analysis in the fast switching regime shows that the total cell density in the
two-population go-or-grow model can be described in terms of a single
reaction-diffusion equation with density-dependent diffusion and
proliferation. Using the connection to single-population models, we study
travelling wave solutions, showing that the wave speed in the go-or-grow model
is always bounded by the wave speed corresponding to the well-known Fisher-KPP
equation.

###### keywords:

travelling waves, go-or-grow, reaction-diffusion, fast reaction limit

[1]organization=Mathematical Institute, addressline=University of Oxford,
postcode=OX2 6GG, city=Oxford, country=United Kingdom

##  1 Introduction

Reaction-diffusion equations are often used to describe the dynamics of
invading species, such as cells, and often exhibit travelling wave solutions.
A classic example of this is the Fisher-Kolmogorov–Petrovsky–Piskunov (FKPP)
equation, which assumes that cells move via linear diffusion and proliferate
following logistic growth [1, 2]. In the last decade however, so-called go-or-
grow models [3, 4, 5, 6, 7] have been suggested to describe the dichotomous
behaviour observed in experiments whereby an individual cell either
proliferates slowly and migrates quickly, or proliferates quickly and migrates
slowly [8, 9]. Under the go-or-grow hypothesis, cells are also able to
transition between the migratory and proliferative states via a phenotypic
switch. The mechanism underlying this switching behaviour is not yet fully
understood [10, 11], but could be due to a variety of factors such as
surrounding cell density, pressure or environmental features, amongst others.
One way to model heterogeneous populations of cells under the go-or-grow
hypothesis is to use systems of reaction-diffusion equations. For specific
forms of the phenotypic switching functions, previous works study
approximations of the travelling wave solutions by combining the method
introduced by Canosa [12] with qualitative theory and numerical simulations
[5, 6]. However, these studies also revealed notable differences between the
minimum travelling wave speed predicted analytically and the speed observed
numerically. In addition, minimal progress has been made to understand how
invasion phenomena depend on the switching functions more generally [7].

In this work, we study invasion phenomena in a minimal go-or-grow model,
whereby cells either proliferate or migrate, and switching between phenotypic
states occurs via general switching functions. After presenting this model in
Section 2, we first demonstrate, via a formal fast reaction limit, the
connection between a model involving heterogeneous populations and a model
consisting of a single population of cells, with density-dependent diffusion
and proliferation terms. Interestingly, this reveals a connection between
multiple phenotype models with density-dependent switching functions, and
degenerate diffusion models consisting of a single phenotype, for which the
standard linear stability analysis does not yield accurate travelling wave
speed predictions [13, 14]. We then proceed in Section 3 to perform a
travelling wave analysis of the minimal go-or-grow model presented in this
work and show that the minimum travelling wave speed emitted from the model
does not exceed that of the FKPP equation. Furthermore, we demonstrate that
for specific forms of the phenotypic switching function, solutions of the
multiple phenotype model do not attain this minimum travelling wave speed.
This provides an analytical explanation for the observed discrepancies in
prior studies between the minimum wave speed predicted analytically and the
numerically observed travelling wave speed [8].

##  2 A minimal go-or-grow model

We consider a heterogeneous cell population consisting of cells which are
either proliferative or invasive. The density of the proliferative population
is represented by $\rho_{1}(\mathbf{x},t)$ and the density of the invasive
population is given by $\rho_{2}(\mathbf{x},t)$, for
$(\mathbf{x},t)\in\mathbb{R}^{d}\times[0,+\infty)$. For the sake of
simplicity, we assume that invasive cells move via linear diffusion, and can
switch to the proliferative phenotype at a rate given by the switching
function $\Gamma_{1}(\rho)$, which depends on the total cell density
$\rho=\rho_{1}+\rho_{2}$. On the other hand, we assume that the proliferative
cells follow logistic growth, and can switch to the invasive phenotype at a
rate given by $\Gamma_{2}(\rho)$. With this, the model reads

$\begin{split}\partial_{t}\rho_{1}&=\Delta\rho_{1}-\rho_{1}\Gamma_{1}(\rho)+\rho_{2}\Gamma_{2}(\rho)\,,\\\
\partial_{t}\rho_{2}&=\rho_{2}\left(1-{\rho}\right)+\rho_{1}\Gamma_{1}(\rho)-\rho_{2}\Gamma_{2}(\rho)\,.\end{split}$
(1)

Following from standard rescaling of space and time, the diffusion coefficient
and the intrinsic proliferation rate in Eqs. (1) have been set to one without
loss of generality. We also note that this model does not account for cell
death, however apoptosis can be readily incorporated into the model, as
demonstrated in [5, 6], without altering the analysis. The results in this
work generalise directly to this scenario.

The model given in Eqs. (1) is studied in [5, 6] for specific choices of
switching functions, which facilitate analytic tractability. Here, and unless
stated otherwise, we keep $\Gamma_{1}(\rho)$ and $\Gamma_{2}(\rho)$ general.
Following the biological intuition that regions of high cell density lead to
contact inhibition or nutrient limitation, and hence to inhibition of
proliferation [15, 16], typical switching functions are based on the
assumption that $\Gamma_{1}(\rho)$ is non-increasing, and $\Gamma_{2}(\rho)$
is non-decreasing. While the model presented here is simpler, it remains
closely related to those examined in [3, 4], which consider density-dependent
diffusion, and in [7], which further incorporates the influence of the
extracellular matrix on cell invasion.

###  2.1 Connecting heterogeneous to single-population models

We investigate the fast phenotypic switching limit of Eqs. (1). In other
words, we assume $\Gamma_{1}(\rho)=\tilde{\Gamma}_{1}(\rho)/\epsilon$, and
$\Gamma_{2}(\rho)=\tilde{\Gamma}_{2}(\rho)/\epsilon$ where
$\tilde{\Gamma}_{1}(\rho),\,\tilde{\Gamma}_{2}(\rho)\sim O(1)$ for
$\epsilon\ll 1$. Upon formally taking the limit $\epsilon\rightarrow 0$ in
Eqs. (1), we obtain
$\rho_{1}\tilde{\Gamma}_{1}(\rho)=\rho_{2}\tilde{\Gamma}_{2}(\rho)$. Consider
now the equation governing the evolution of the total cell density over time,

$\partial_{t}\rho=\Delta\rho_{1}+\rho_{2}(1-\rho)\,.$ (2)

In the limit $\epsilon\rightarrow 0$, the two cell densities $\rho_{1}$ and
$\rho_{2}$ can be related to the total cell density via

$\rho_{1}=\frac{\Gamma_{2}(\rho)}{\Gamma_{1}(\rho)+\Gamma_{2}(\rho)}\,\rho\,,\quad\rho_{2}=\frac{\Gamma_{1}(\rho)}{\Gamma_{1}(\rho)+\Gamma_{2}(\rho)}\,\rho\,.$
(3)

Combining Eqs. (2) and (3) we obtain the following equation for the evolution
of the total cell density

$\partial_{t}\rho=\nabla\cdot\left(D(\rho)\nabla\rho\right)+r(\rho)\rho\,,$
(4)

with the density-dependent diffusion coefficient $D(\rho)$, and proliferation
rate $r(\rho)$ given by

$D(\rho)=\frac{\Gamma_{1}(\rho)\Gamma_{2}^{\prime}(\rho)-\Gamma_{1}^{\prime}(\rho)\Gamma_{2}(\rho)}{\left(\Gamma_{1}(\rho)+\Gamma_{2}(\rho)\right)^{2}}\,\rho+\frac{\Gamma_{2}(\rho)}{\Gamma_{1}(\rho)+\Gamma_{2}(\rho)},\qquad
r(\rho)=\frac{\Gamma_{1}(\rho)(1-\rho)}{\Gamma_{1}(\rho)+\Gamma_{2}(\rho)}\,.$

In this case, assuming that $\Gamma_{1}$ is a non-increasing function of the
total cell density is sufficient to show that the diffusion coefficient does
not become negative, i.e. $D(\rho)\geq 0$ for any $\rho\geq 0$; and also that
the per capita proliferation rate is a non-increasing function of density,
i.e. $r^{\prime}(\rho)\leq 0$. As mentioned earlier, this assumption on
$\Gamma_{1}$ aligns with basic biological intuition. Moreover, under the
assumption that cells remain proliferative at low densities, so that
$\Gamma_{2}(0)=0$, we obtain that the diffusion coefficient follows
$D(\rho)\sim 2\Gamma_{2}^{\prime}(0)/\Gamma_{1}(0)\rho+O(\rho^{2})$, which
resembles typical non-linear equations with degenerate diffusion.

Finally, we highlight that in the case of constant switching rates
($\Gamma_{1}^{\prime}=\Gamma_{2}^{\prime}=0$) we recover the well-known FKPP
equation [13]

$\partial_{t}\rho=\theta\Delta\rho+(1-\theta)\rho(1-\rho)\,,$ (5)

where $\theta\in(0,1)$ is determined by the constant ratio
$\Gamma_{1}/\Gamma_{2}$. Equations of this type often exhibit travelling wave
solutions $\rho(\mathbf{x},t)=U(\mathbf{r}\cdot\mathbf{x}-ct)$, where
$\mathbf{r}$ is the direction of propagation, and $c$ corresponds to the wave
speed. In particular, the FKPP equation (Eq. (5)) admits travelling wave
solutions with constant speed $c\geq c_{\text{min}}=2\sqrt{\theta(1-\theta)}$.
We note that $c_{\text{min}}\leq 1$ for $\theta\in(0,1)$, which we use as a
reference to compare with the speed of propagation of travelling wave
solutions to Eqs. (1).

##  3 Travelling wave analysis for heterogeneous populations

We now revisit Eqs. (1) to analyse travelling wave solutions of the form:
$\rho_{1}(\mathbf{x},t)=U_{1}(\mathbf{r}\cdot\mathbf{x}-ct)$,
$\rho_{2}(\mathbf{x},t)=U_{2}(\mathbf{r}\cdot\mathbf{x}-ct)$. Fig. 1 shows
travelling wave solutions of Eqs. (1) subject to different choices of the
switching functions. By defining
$z=\mathbf{r}\cdot\mathbf{x}-ct\in(-\infty,+\infty)$, and writing Eqs. (1) in
travelling wave coordinates we obtain

$\begin{split}c\frac{\mathrm{d}U_{1}}{\mathrm{d}z}+\frac{\mathrm{d}^{2}U_{1}}{\mathrm{d}z^{2}}-U_{1}\Gamma_{1}(U)+U_{2}\Gamma_{2}(U)&=0\,,\\\
c\frac{\mathrm{d}U_{2}}{\mathrm{d}z}+U_{2}(1-U)+U_{1}\Gamma_{1}(U)-U_{2}\Gamma_{2}(U)&=0\,,\end{split}$
(6)

where $U(z)=U_{1}(z)+U_{2}(z)$. Next, we look at steady states of the system
(6).

Figure 1: Numerical solutions to Eqs. (1) subject to the initial conditions:
$\rho_{1}(x,0)=\rho_{2}(x,0)=0.2$ for $x<100$, and
$\rho_{1}(x,0)=\rho_{2}(x,0)=0$ for $x\geq 100$, and for different choices of
the switching functions $\Gamma_{1},\,\Gamma_{2}$. Solutions are plotted at
times $t=200,300,400,500$. Chosen switching functions are as follows: (a)
$\Gamma_{1}(\rho)=0.5,\,\Gamma_{2}(\rho)=1$; (b)
$\Gamma_{1}(\rho)=0.5,\,\Gamma_{2}(\rho)=1.5\rho$; (c)
$\Gamma_{1}(\rho)=0.5(1-\rho),\,\Gamma_{2}(\rho)=1.5\rho$;
(d)$\Gamma_{1}(\rho)=0.5(1-\rho^{2}/(0.5^{2}+\rho^{2})),\,\Gamma_{2}(\rho)=1.5\rho^{2}/(0.5^{2}+\rho^{2})$.

First note that the origin $(U_{1},U_{2})=(0,0)$ is always a steady state, and
we assume $U_{1},U_{2}\rightarrow 0$ as $z\rightarrow+\infty$. The other
steady states, corresponding to $z\rightarrow-\infty$, depend on the choice of
phenotypic switching functions (see Fig. 1) and belong to the set
$\left\\{(U_{1},U_{2})\in[0,1]^{2}:U_{1}\Gamma_{1}(U)=U_{2}\Gamma_{2}(U),\;U_{2}(1-U)=0\right\\}$.
In particular, a steady state with $U_{2}=0$ can only occur if there exists a
$U^{*}>0$ such that $\Gamma_{1}(U^{*})=0$. In this case, we obtain

$(U_{1},U_{2})=(U^{*},0)\,.$

If $\Gamma_{1}(U)>0$ for all $U>0$, then we obtain

$(U_{1},U_{2})=\left(\frac{\Gamma_{2}(1)}{\Gamma_{1}(1)+\Gamma_{2}(1)},\frac{\Gamma_{1}(1)}{\Gamma_{1}(1)+\Gamma_{2}(1)}\right)\,.$

###  3.1 Dispersion relation and minimum travelling wave speed

When seeking an expression for the travelling wave speed, $c$, it is often
commonplace to perform a linear stability analysis centred around the origin,
which often yields a lower bound on the speed: $c\geq c_{\text{min}}$, for
some minimum wave speed $c_{\text{min}}\geq 0$. Unfortunately, for the model
presented in this work, analytical progress following this approach is limited
and instead we shift the investigation to identify the dispersion relation,
which relates the travelling wave speed with the dynamics at the leading
front. We observe that at the leading front, cell densities are small, and
hence we can linearise Eqs. (6) to give

$\begin{split}c\frac{\mathrm{d}U_{1}}{\mathrm{d}z}+\frac{\mathrm{d}^{2}U_{1}}{\mathrm{d}z^{2}}-\gamma_{1}U_{1}+\gamma_{2}U_{2}&=0\,,\\\
c\frac{\mathrm{d}U_{2}}{\mathrm{d}z}+\gamma_{1}U_{1}+(1-\gamma_{2})U_{2}&=0\,,\end{split}$
(7)

where we denote $\gamma_{i}=\Gamma_{i}(0)$ for $i=1,2$. By substituting the
ansatz $U_{1}(z)\sim Ae^{-\sigma z}$, $U_{2}(z)\sim Be^{-\sigma z}$, for $z\gg
1$ and $A,B\geq 0$ constants, we obtain

$\left(\;\smash[b]{\underbrace{\begin{matrix}-c\sigma+\sigma^{2}-\gamma_{1}&\gamma_{2}\\\
\gamma_{1}&-c\sigma+1-\gamma_{2}\end{matrix}}_{M}}\;\right)\vphantom{\underbrace{\begin{matrix}-c\sigma+\sigma^{2}-\gamma_{1}&\gamma_{2}\\\
\gamma_{1}&-c\sigma+1-\gamma_{2}\end{matrix}}_{M}}\begin{pmatrix}A\\\
B\end{pmatrix}=0\,.$

In order to ensure $A,\,B\neq 0$, we require $\det M=0$, which gives

$c^{2}-c\left(\sigma+\frac{1-\gamma_{1}-\gamma_{2}}{\sigma}\right)+1-\gamma_{2}-\frac{\gamma_{1}}{\sigma^{2}}=0\,.$

The travelling wave speed, $c$, can then be solved for as a function of
$\sigma$, giving the dispersion relation:

$c(\sigma)=\frac{1}{2}\left(f(\sigma)+\sqrt{f(\sigma)^{2}-4(1-\gamma_{2})+\frac{4\gamma_{1}}{\sigma^{2}}}\,\right),\quad
f(\sigma)=\sigma+\frac{1-\gamma_{1}-\gamma_{2}}{\sigma}.$ (8)

The expression above predicts the existence of a minimum travelling wave
speed, $c_{\text{min}}$ for $\sigma=\sigma^{*}$, such that
$c_{\text{min}}=c(\sigma^{*})$. Finding the minimum speed analytically seems
challenging; instead we observe that $c(1)=1$, and by differentiating
$c(\sigma)$ with respect to $\sigma$ and evaluating at $\sigma=1$, we obtain

$\left(\frac{\mathrm{d}c}{\mathrm{d}\sigma}\right)_{\sigma=1}=\frac{\gamma_{2}-\gamma_{1}}{\gamma_{1}+\gamma_{2}}\,.$

Thus, for $\gamma_{1}=\gamma_{2}$, $\sigma^{*}=1$ yields the minimum wave
speed, and in this case it is given by $c_{\text{min}}=1$. On the other hand,
when $\gamma_{1}>\gamma_{2}$ we have
$({\mathrm{d}c}/{\mathrm{d}\sigma})_{\sigma=1}<0$, and $\sigma^{*}>1$ with
$c_{\text{min}}<1$. An analogous argument yields $\sigma^{*}<1$ and
$c_{\text{min}}<1$ for $\gamma_{1}<\gamma_{2}$. Hence, the minimum travelling
wave speed displayed by Eqs. (1) is, at most, the minimum travelling wave
speed corresponding to the FKPP model, Eq. (5). Fig. 2a shows an excellent
agreement between the numerically observed speed and the prediction from Eq.
(8) for a fully linear model with constant switching functions. Next, we
explore whether this agreement holds for more general nonlinear models.

Figure 2: Relationship between the estimated numerical speed (dots) of
travelling wave solutions to Eqs. (1), and the minimum wave speed predicted by
Eq. (8) (dashed lines). Numerical wave speed estimated using a trapezoidal
rule applied to the expression
$c=\int_{-\infty}^{+\infty}U_{2}(z)(1-U(z))\,\mathrm{d}z$, which follows from
Eq. (6). Eqs. (1) are solved on a domain of length 7000, for a maximum
simulation time of 6500, and initial conditions given by
$\rho_{1}(x,0)=\rho_{2}(x,0)=0.2$ if $x<100$ and
$\rho_{1}(x,0)=\rho_{2}(x,0)=0$ if $x\geq 100$. Chosen switching functions are
as follows: (a) $\Gamma_{1}(\rho)=\gamma_{1},\,\Gamma_{2}(\rho)=\gamma_{2}$;
(b) $\Gamma_{1}(\rho)=\gamma_{1},\,\Gamma_{2}(\rho)=\beta\rho$; (c)
$\Gamma_{1}(\rho)=\gamma_{1}(1-\rho^{2}/(0.5^{2}+\rho^{2})),\,\Gamma_{2}(\rho)=\beta\rho^{2}/(0.5^{2}+\rho^{2})$.

###  3.2 Explicit expressions for the minimum travelling wave speed

Fully analytical expressions for $c_{\text{min}}$ can be obtained under the
assumption that $\gamma_{1}=0$ or $\gamma_{2}=0$. We present first the
$\gamma_{2}=0$ case as it represents a more biologically realistic scenario
where cells remain proliferative at low densities. In this case, Eq. (8)
simplifies to

$c(\sigma)=\frac{1}{2}\left(\sigma+\frac{1-\gamma_{1}}{\sigma}+\frac{\left|\sigma^{2}-1-\gamma_{1}\right|}{\sigma}\right)\,.$

Note then that $c(\sigma)$ increases for $\sigma>\sqrt{1+\gamma_{1}}$, and
decreases for $\sigma<\sqrt{1+\gamma_{1}}$. Hence,
$\sigma^{*}=\sqrt{1+\gamma_{1}}$, and

$c_{\text{min}}=\frac{1}{\sqrt{1+\gamma_{1}}}\,.$ (9)

When $\gamma_{1}=0$, a similar argument yields the minimum wave speed

$\displaystyle c_{\text{min}}=\begin{cases}\sqrt{1-\gamma_{2}}&\mbox{for
}\gamma_{2}<1\,,\\\ 0&\mbox{for }\gamma_{2}>1\,.\end{cases}$

###  3.3 The minimum wave speed is not attainable when $\Gamma_{2}(0)=0$

Based on the behaviour observed in Fig. 2b-c, we focus on the case
$\gamma_{2}=\Gamma_{2}(0)=0$, with $\Gamma_{2}\neq 0$ to ensure that
travelling wave solutions exist. By linearising Eqs. (6) subject to
$\Gamma_{2}(0)=0$, representing the region near the invading front where cell
densities are low, we obtain

$\begin{split}c\frac{\mathrm{d}U_{1}}{\mathrm{d}z}+\frac{\mathrm{d}^{2}U_{1}}{\mathrm{d}z^{2}}-\gamma_{1}U_{1}&=0\,,\\\
c\frac{\mathrm{d}U_{2}}{\mathrm{d}z}+\gamma_{1}U_{1}+U_{2}&=0\,.\end{split}$

Solving this system of equations, subject to $U_{1},\,U_{2}\rightarrow 0$ as
$z\rightarrow\infty$, yields

$U_{1}(z)=Ae^{-\lambda(c)z},\quad
U_{2}(z)=\frac{\gamma_{1}}{c\lambda(c)-1}\,Ae^{-\lambda(c)z},\quad\lambda(c)=\frac{1}{2}\left(c+\sqrt{c^{2}+4\gamma_{1}}\right)\,,$

for an arbitrary constant $A>0$. Positivity of travelling wave solutions (i.e.
$U_{1},U_{2}\geq 0$ for all $z\in\mathbb{R}$) requires $c\lambda(c)-1>0$,
which is only satisfied when $c>c_{\text{min}}$ (as in Eq. (9)). In
particular, we obtain that travelling wave solutions for $\Gamma_{2}(0)=0$
exhibit a wave speed exceeding $c_{\text{min}}$, i.e. $c>c_{\text{min}}$ when
$\Gamma_{2}(0)=0$.

##  4 Discussion and open problems

In this work, we examined a minimal go-or-grow model, where cells can either
proliferate following logistic growth or move through linear diffusion. Via a
formal fast reaction limit in the fast phenotypic switching regime, we
revealed a connection between heterogeneous models and single-population
models. This relationship implies an approximate equivalence between discrete-
state representations of cell populations and models where diffusion and
proliferation rates are determined by density. Moreover, our analysis suggests
that the invasion speed in travelling wave solutions of heterogeneous models
is bounded by the wave speed of solutions to single-populations models. We
highlight that this phenomenon is likely to change when energy or resource
constraints associated to diffusion and proliferation are taken into account
[7]. More research is needed, however, to better understand how to model such
energy constraints and their impact on cell invasion.

Furthermore, our travelling wave analysis also reveals a lower bound on the
invasion speed of travelling wave solutions to go-or-grow models. This
estimate, based on a linear stability analysis, breaks down for nonlinear
models (see Fig. 2b-c), as observed in similar models [5, 14, 17, 18].
Numerically, we observe that $c_{\text{min}}\leq c\leq 1$ for the model given
by Eqs. (6), and that $c\rightarrow 1$ as the nonlinear switching rates
increase. Although this behaviour intuitively makes sense when considering the
connection to single-population models, we do not currently have a rigorous
proof. We leave these conjectures for further investigation.

## Acknowledgments

RMC is supported by the Engineering and Physical Sciences Research Council
(EP/T517811/1) and the Oxford-Wolfson-Marriott scholarship at Wolfson College,
University of Oxford. This work was supported by a grant from the Simons
Foundation (MP-SIP-00001828, REB).

## References

  * [1] R. A. Fisher.  The wave of advance of advantageous genes.  Annals of Eugenics, 7(4):355–369, 1937.
  * [2] A. N. Kolmogorov, I. Petrovskii, and N. S. Piskunov.  A study of the equation of diffusion with increase in the quantity of matter, and its application to a biological problem.  Moskovskogo University Bulletin of Mathematics, 1(6):1–25, 1937\.
  * [3] P. Gerlee and S. Nelander.  The impact of phenotypic switching on glioblastoma growth and invasion.  PLoS Computational Biology, 8(6):e1002556, 2012.
  * [4] P. Gerlee and S. Nelander.  Travelling wave analysis of a mathematical model of glioblastoma growth.  Mathematical Biosciences, 276:75–81, 2016.
  * [5] T. L. Stepien, E. M. Rutter, and Y. Kuang.  Traveling waves of a go-or-grow model of glioma growth.  SIAM Journal on Applied Mathematics, 78(3):1778–1801, 2018.
  * [6] A. Tursynkozha, A. Kashkynbayev, B. Shupeyeva, E. M. Rutter, and Y. Kuang.  Traveling wave speed and profile of a “go or grow” glioblastoma multiforme model.  Communications in Nonlinear Science and Numerical Simulation, 118:107008, 2023.
  * [7] R. M. Crossley, K. J. Painter, T. Lorenzi, P. K. Maini, and R. E. Baker.  Phenotypic switching mechanisms determine the structure of cell migration into extracellular matrix under the ’go-or-grow’ hypothesis.  arXiv preprint arXiv:2401.07279, 2024.
  * [8] A. M. Stein, T. Demuth, D. Mobley, M. Berens, and L. M. Sander.  A mathematical model of glioblastoma tumor spheroid invasion in a three-dimensional in vitro experiment.  Biophysical Journal, 92(1):356–365, 2007.
  * [9] K. S. Hoek, O. M. Eichhoff, N. C. Schlegel, U. Döbbeling, N. Kobert, L. Schaerer, S. Hemmi, and R. Dummer.  In vivo switching of human melanoma cells between proliferative and invasive states.  Cancer Research, 68(3):650–656, 2008.
  * [10] D. Tosh and J. M. W. Slack.  How cells change their phenotype.  Nature Reviews Molecular Cell Biology, 3(3):187–194, 2002.
  * [11] P. J. Quesenberry and J. M. Aliotta.  Cellular phenotype switching and microvesicles.  Advanced Drug Delivery Reviews, 62(12):1141–1148, 2010.
  * [12] J. Canosa.  On a nonlinear diffusion equation describing population growth.  IBM Journal of Research and Development, 17(4):307–313, 1973.
  * [13] J. D. Murray.  Mathematical Biology I: An Introduction.  Springer New York, 2001.
  * [14] M. El-Hachem, S. W. McCue, and M. J. Simpson.  Travelling wave analysis of cellular invasion into surrounding tissues.  Physica D: Nonlinear Phenomena, 428:133026, 2021.
  * [15] A. Puliafito, L. Hufnagel, P. Neveu, S. Streichan, A. Sigal, D. K. Fygenson, and B. I. Shraiman.  Collective and single cell behavior in epithelial contact inhibition.  Proceedings of the National Academy of Sciences, 109(3):739–744, 2012.
  * [16] C. Falcó, D. J. Cohen, J. A. Carrillo, and R. E. Baker.  Mechanical constraints and cell cycle regulation in models of collective cell migration.  arXiv preprint arXiv:2401.08805, 2024.
  * [17] C. Colson, F. Sánchez-Garduño, H. M. Byrne, P. K. Maini, and T. Lorenzi.  Travelling-wave analysis of a model of tumour invasion with degenerate, cross-dependent diffusion.  Proceedings of the Royal Society A, 477(2256):20210593, 2021.
  * [18] R. M. Crossley, P. K. Maini, T. Lorenzi, and R. E. Baker.  Traveling waves in a coarse-grained model of volume-filling cell invasion: Simulations and comparisons.  Studies in Applied Mathematics, 151(4):1471–1497, 2023.