# Higher order assortativity for directed weighted networks and Markov chains

Alberto Arcagni1 Roy Cerqueti2,3 Rosanna Grassi4,∗

###### Abstract

This paper proposes a new class of assortativity measures for weighted and
directed networks. We extend the classical Newman’s degree-degree
assortativity by considering nodes’ attributes different from the degree.
Moreover, we propose connections among the nodes through directed paths of
length greater than one – thus, obtaining higher order assortativity. We
provide an empirical application of these measures for the paradigmatic case
of the trade network. Importantly, we show how this global network indicator
is strongly related to the autocorrelations of the states of a Markov chain.

Keywords: Networks, Higher order Assortativity, Random walks, Markov Chains.

1 MEMOTEF Department, Sapienza University of Rome, Via del Castro Laurenziano,
9 - 00161, Rome, Italy
<EMAIL_ADDRESS>
2 Department of Economic and Social Sciences, Sapienza University of Rome,
Piazzale Aldo Moro, 5 - 00185, Rome, Italy
<EMAIL_ADDRESS>
3 GRANEM, Université d’Angers, SFR CONFLUENCES, F-49000 Angers, France
∗ Corresponding author
4 Department of Statistics and Quantitative Methods, University of Milano -
Bicocca, Via Bicocca degli Arcimboldi, 8 - 20126 Milan, Italy
<EMAIL_ADDRESS>
##  1 Introduction

Asymmetric interactions characterize many complex systems in nature. Directed
networks are suitable tools to represent these complex situations. Despite the
remarkable persistence of this kind of systems in real-world phenomena, a
formal representation of the directed networks is still not used, probably due
to the difficulty in the mathematical treatment of such models. Indeed, some
network indicators are not easy to manage in the directed framework;
sometimes, and in some cases, a proper definition of the network measures for
the directed networks is still missing. Under this perspective, this paper
focuses on the issue of providing a generalized concept of the assortativity
measure for directed networks. The theme is relevant because assortativity is
a global indicator that provides useful insights about the network structure.
In the classical definition of Newman (2002), the assortativity is represented
by a global measure based on the Pearson correlation between the degrees of
nodes. This measure ranges in the interval $[-1,1]$, having positive values in
assortative networks and negative values in disassortative ones. The
interpretation of the degree-degree assortativity measure is simple:
undirected networks, high assortativity means that high-degree nodes tend to
connect to other high-degree nodes. On the contrary, in disassortative
networks, high-degree nodes tend to be connected to low-degree nodes. The
assortativity measure can be extended into two different directions. On the
one side, one can consider other quantitative attributes of the nodes
different from the degree; on the other side, one can move from the adjacency
of the nodes – which is the basis of Newman’s degree-degree assortativity –
and propose more general ways to connect them. We move in the direction of
meeting both such important requirements. Indeed, in the former case, degree-
based assortativity provides a restrictive view of the accordance among nodes.
A generalization to different nodes’ attributes overcomes this issue,
extending the meaning of the coefficient to describe
similarities/dissimilarities between nodes other than the degree (see
Meghanathan (2016), Arcagni et al. (2019), Yuan et al. (2021)). In the latter
case, it is relevant to mention the generalization of Newman’s assortativity
index proposed in Arcagni et al. (2017). The authors provide a closed formula
that leads to a unifying approach of the assortativity for undirected and
unweighted networks; in particular, they consider pairs of vertices not
necessarily adjacent but also connected through paths, shortest paths and
random walks. Importantly, connecting nodes through paths of length greater
than one allows us to gain insights on the similarities between connected, but
not adjacent nodes. Indeed, high-degree nodes may not be directly connected to
each other, but remain equally connected through low-degree ones.

In Arcagni et al. (2017), the new concept of assortativity coefficient based
on paths is labelled as higher order assortativity. The higher order
assortativity allows us to recover, in particular cases, the Newman index as
well as other measures based on assortativity (Estrada et al., 2008; Mayo et
al., 2015; Van Mieghem et al., 2010), using suitable definitions of the matrix
governing the connections. In the same line, Arcagni et al. (2019) extended
the higher order assortativity to a new class of higher order measures
specific for weighted undirected networks, recovering as a particular case the
weighted assortativity index introduced in Leung and Chau (2007). The problem
of assortativity in directed networks has also been investigated (Yuan et al.,
2021). In particular, Foster et al. (2010); Piraveenan et al. (2010) propose a
natural modification of the Newman formula, encoding the assortativity measure
into four directed measures on the basis of the combination of nodes’ in- and
out-degrees. The coefficients quantify the tendency of nodes with high in-
degree to be connected to nodes with high in-degree – and the same applies to
the other cases in-degree with out-degree, out-degree with in-degree and out-
degree with out-degree. However, a higher order measure of assortativity
tailored for weighted and direct networks is still missing.

This paper contributes to this debate by filling this gap. We introduce a new
concept of assortativity measure of order $h$ – with $h\geq 1$ – for the
general case of weighted and directed networks. According to the previous
literature and when needed, we briefly refer to such a concept as the higher
order assortativity measure. We here complete the definitions introduced in
Arcagni et al. (2017, 2019) by including also the weighted directed networks.
Moreover, the proposed concept is formulated on a node attribute that is not
necessarily the degree or the strength. In doing so, we extend the Newman
assortativity coefficient in both the directions described above.

It is important to notice that the proposed assortativity measure offers a
relevant interpretation in the field of the Markov chains (see, e.g., Norris
et al., 1998, for a survey on this class of stochastic processes). Indeed, the
directed and weighted networks are particularly suitable for the description
of random moving flows over a discrete set of states identified by nodes. In
the graph representation of a Markov chain, the existence of an arc connecting
state/node $i$ with state/node $j$ is associated with a non-null probability
of jumping from $i$ to $j$. Therefore, the (suitably normalized) weights of
the arcs in a complex network are naturally associated with the probability to
observe a flow from a state to another one in the corresponding Markov chain.
In stating the bridge between complex networks and Markov chains, we follow
the route traced by some relevant contributions in the literature. It is worth
mentioning Gómez et al. (2010), where the authors deal with a disease
propagation model by merging complex networks and Markov chains. In the same
application environment, Iannelli et al. (2017) elaborate on the informative
content of the general network-based measures derived from Markov chain
theory.

Our novelty here is to connect the higher order assortativity – a topological
property of the network involving nodes beyond the adjacent ones – and
homogeneous discrete-time Markov chains.

Specifically, let us consider the nodes’ assortativity referred to a given
vertex centrality measure. The assortativity index of order $h$ ranges in
$[-1,1]$; it approaches $1$ ($-1$, resp.) when the nodes connected through
walks of length $h$ have a high level of concordance (discordance,
respectively) in the considered centrality measure. Such a variation range
suggests that the assortativity index of order $h$ can be interpreted as an
autocorrelation term at $h$ lags in the context of Markov chains. This is
exactly our case. Indeed, the order $h$ is the length of a walk from one node
to another. If we imagine that each link takes one unit of time, then the
considered walk takes $h$ units of time. In this respect, the order $h$ can be
seen as a temporal variable for a discrete-time Markov chain with states given
by the nodes of the network along with their centrality measures and
transition probabilities opportunely related to the weights of the network’s
links.

This interpretation of the assortativity measure of order $h$ allows stating a
natural bridge between complex networks and stochastic processes. In so doing,
we are able to move from the information content of the higher order
assortativity of a network to the dynamical properties of the underlying
Markov chain. Specifically, the temporal dimension of the network and the
regularities captured by the autocorrelation – which are hidden in the network
structure – become clear in moving to the Markov chain theory.

The paper is structured as follows. In Section 2 we describe the mathematical
notations. In Section 3 we introduce a novel definition of higher order
assortativity extended to weighted digraphs. In Section 4 we show an
application of the proposed index to the trade network. Section 5 illustrates
how the higher-order assortativity measure can be interpreted in the context
of the discrete-time homogeneous Markov chains with finite states. Conclusions
follow in Section 6.

##  2 Preliminaries

In this section we outline of some basic definitions and notation related to
directed complex networks that are used to present our methodological
proposal.

For a more detailed treatment we refer, for instance, to Bang-Jensen and Gutin
(2008); Newman (2010); Gross and Yellen (2003).

A directed graph (or digraph) $N=(V,A)$ is a pair of sets $V$ and $A$, where
$V$ is the set of $n$ vertices (or nodes) and $A$ is the set of $m$ ordered
pairs (arcs) of vertices of $V$; if $(i,j)$ and/or $(j,i)$ is an element of
$A$, then vertices $i$ and $j$ are adjacent. A $i-j$ directed walk is a
sequence of vertices and arcs from $i$ to $j$ such that all arcs have the same
direction.

A weight $w_{ij}>0$ can be associated with each arc $(i,j)$ so that a weighted
digraph is obtained. Moreover, it is assumed that $w_{ij}=0$ if and only if
$(i,j)\notin A$. The weights $w$’s are collected in a real $n$-square matrix W
(the weighted adjacency matrix). By definition, the elements of such a matrix
describe both the adjacency relationships between the vertices in $V$ and the
weights on the arcs. In this context, a weighted directed network is a
directed graph $N=(V,A)$ equipped by a weighted not symmetric matrix W. In the
unweighted case, non-null weights can take only unitary value, so that matrix
W provides only information about the adjacency relationships. In the case of
unweighted networks, we denote W by A (the adjacency matrix).

The $(i,j)$ element of the $k-$power of the matrix A is the number of directed
walks of length $k$ from $i$ to $j$.

Since in-degree, out-degree, in-strength and out-strength are the most popular
centrality measures in the context of the assortativity – as the original
definition of Newman (2002) suggests – we introduce a particular notation for
them, along with their definition.

We define the in-degree $d^{I}_{i}$ (respectively out-degree $d^{O}_{i}$) of a
node $i\in V$ as the number of arcs pointing towards (respectively starting
from) $i$. Then, the in-degree and out-degree vectors are, respectively
$\mathbf{d}^{I}$, and $\mathbf{d}^{O}$. The vector of the vertices’ degrees is
given by $\mathbf{d}=\mathbf{d}^{I}+\mathbf{d}^{O}$. Analogously, definitions
of the in-strength and out-strength of node $i\in V$ ($s^{I}_{i}$ and
$s^{O}_{i}$, respectively) and in-strength, out-strength vectors
($\mathbf{s}^{I}$ and $\mathbf{s}^{O}$, respectively) are provided. A vertex
$i\in V$ with $d^{I}_{i}=0$ and $d^{O}_{i}>0$ is called _source_ , as it is
the origin of its outcoming arcs. Similarly, a vertex $i\in V$ with
$d^{O}_{i}=0$ and $d^{I}_{i}>0$ is called _sink_ since it is the end of each
incoming arc.

##  3 A novel definition of higher order assortativity

We consider two centrality measures for the vertices of a weighted directed
network and collect their value the definition of higher order assortativity
for weighted digraphs, thus extending the related concepts introduced by
Arcagni et al. (2017, 2019) as a particular application of the Pearson
correlation index in the $n$-dimensional vectors x and y. The higher order
assortativity index of length $h\geq 1$ for a directed graph and associated
with x and y can be defined as follows

$r_{h}=r(\textbf{x},\textbf{y},\mathbf{E}_{h})=\frac{\mathbf{x}^{\top}\left(\mathbf{E}_{h}-\mathbf{p}_{h}\mathbf{q}_{h}^{\top}\right)\mathbf{y}}{\sqrt{\left[\mathbf{x}^{\top}\left(\mathbf{D_{p}}_{h}-\mathbf{p}_{h}\mathbf{p}_{h}^{\top}\right)\mathbf{x}\right]\left[\mathbf{y}^{\top}\left(\mathbf{D_{q}}_{h}-\mathbf{q}_{h}\mathbf{q}_{h}^{\top}\right)\mathbf{y}\right]},}$
(1)

where $\textbf{E}_{h}$ is a normalized $n$-squared matrix whose element on the
$i$-th row and $j$-th column is denoted by $e^{(h)}_{ij}$, with
$\|\textbf{E}_{h}\|_{1}=\sum_{i,j=1}^{n}\left|e_{ij}\right|=1$ while
$\textbf{p}_{h}=\textbf{E}_{h}\textbf{1}$,
$\textbf{q}_{h}=\textbf{E}_{h}^{\top}\textbf{1}$, where 1 is the conformable
vector of elements equal to 1, and $\mathbf{D_{p}}_{h}$ and
$\mathbf{D_{q}}_{h}$ are the diagonal matrices with the vectors
$\textbf{p}_{h}$ and $\textbf{q}_{h}$, respectively, on the main diagonal.

In particular we assume that
$\textbf{E}_{1}=\frac{\textbf{W}}{\|\textbf{W}\|_{1}}$ and, as we will see in
Section 5.2, matrix $\textbf{E}_{h}$ is also related to Markov chains
$\forall\,h\geq 1$. Therefore, the assortativity measure in equation (1)
depends on suitable centrality measures x and y, that are related to the
initial and terminal state, respectively, of a stochastic process of $h$
steps.

Accordingly to Newman (2002), we point out that the in-strength and out-
strength of the nodes plays a relevant role when dealing with assortativity
also in the context of weighted digraphs. Indeed, such centrality measures
provide distinct information about the network topology. Therefore, we
describe such cases in detail and – at the same time – we also give the
interpretation of the higher-order assortativity measure in this particular
context.

Specifically, for a complete treatment of the higher-order assortativity
measure for weighted and directed networks, we consider all the possible
combinations of in-strength and the out-strength and, starting from equation
(1), we define the four higher-order assortativity indexes for weighted
digraph.

  1. 1.

$r^{OI}_{h}=r(\textbf{s}^{O},\textbf{s}^{I},\mathbf{E}_{h})$ (out-in)
This coefficient evaluates the correlation between the out-strength and in-
strength vectors using formula (1), for nodes that are connected through a
directed walk of length $h$. Specifically, a high value of $r^{OI}_{h}$ is
associated to high concordance between the out-strength and the in-strength of
such vertices, while a value of $r^{OI}_{h}$ close to -1 relies to high
discordance.

  2. 2.

$r^{IO}_{h}=r(\textbf{s}^{I},\textbf{s}^{O},\mathbf{E}_{h})$ (in-out)
This index measures the concordance and discordance between the in-strength of
the starting vertices with the out-strength of the arrival vertices after a
directed walk of length $h$.

  3. 3.

$r^{II}_{h}=r(\textbf{s}^{I},\textbf{s}^{I},\mathbf{E}_{h})$ (in-in)
This coefficient provides an evaluation of the correlation between the in-
strength of the starting vertices with the in-strength of the ending vertices
in presence of a directed walk of length $h$.

  4. 4.

$r^{OO}_{h}=r(\textbf{s}^{O},\textbf{s}^{O},\mathbf{E}_{h})$ (out-out)
This higher order assortativity index captures the correlation between the
out-strength of the starting vertices with the out-strength of the terminal
ones when they are connected through a directed walk of length $h$.

Notice that, given $h$, not necessarily $r^{OI}_{h}=r^{IO}_{h}$ because the
indexes depend the distribution of directed walks described by the matrix
$\textbf{E}_{h}$ that is not necessarily symmetric.

##  4 Higher order assortativity in the world trade network

To illustrate the advantages of our definition we show an application of the
higher order assortativity to the trade network.
Data refer to traded goods between countries, the volume of trades are
expressed in US dollars111Data are free and available for the download by
ComTrade - see UN COMTRADE (2022). We constructed both export and import
weighted networks where vertices are countries and weights on arcs represent
the volume of imported or exported goods between 222 countries in the year
2021. Notice that the number of vertices is greater than the recognized
countries, because ComTrade includes even almost uninhabited territories.
Therefore, in the following, the word “country” covers a broad definition that
also includes these regions.

Concerning links between countries, the export network is expected to be
complementary to the import network, as an import from $i$ to $j$ should be
equal to the export from $j$ to $i$. This is not the case because data
describe trades officially declared by _reporters_ countries with their
_partners_ , and some declarations are missed or mismatch. This problem still
persists also considering years before 2021.

Here we followed an approach similar to the one adopted by Antonietti et al.
(2022), that constructed a symmetric trade network aggregating import and
export data. However, unlike Antonietti et al. (2022), in our network we
preserve the arc directions, working in this way with a directed network.

More specifically, let $\mathbf{W}_{I}$ and $\mathbf{W}_{E}$ be weighted
adjacency matrices representing the volume of import and export trades among
countries. Thus, the trade network is defined as $N=(V,\mathbf{W})$ where the
elements $w_{ij}$ of the matrix $\mathbf{W}$ are:

$w_{ij}=\max(w_{ij}^{E},w_{ji}^{I}),$

being $w_{ij}^{E}$ and $w_{ji}^{I}$ the $ij-$ entry of matrices
$\mathbf{W}_{E}$ and $\mathbf{W}_{I}^{\top}$, respectively. In this way, the
asymmetric matrix $\mathbf{W}$ expresses the flow of goods from country $i$ to
$j$. It is worth noting that in constructing $\mathbf{W}$ we choose the
maximum between the import and export. Indeed, in line with Antonietti et al.
(2022), we assume that the higher value declared, the higher the quality of
the information that can be contained.

###  4.1 Analysis of the main characteristics of the network

At first, we analyse some classical network characteristics.

The number of vertices is $n=222$. Let $\mathbf{A}$ be the adjacency matrix
obtained by $\mathbf{W}$ by neglecting arcs’ weights. The network _density_ is
computed as the ratio between the total number of arcs and the maximum number
of potential arcs222$N$ is a directed network without self-loops $(n-1)n$.

Therefore, the density is

$\frac{\|\textbf{A}\|_{1}}{(n-1)n}=\frac{\sum_{i,j=1}^{n}a_{ij}}{(n-1)n}=0.4996.$

Although one might expect the trade network to be complete, this is not true
in this case. This is due to the presence of small countries and regions
representing peripheral nodes, that probably communicate only a part of their
exports or imports.

Distributions of centrality measures, in particular degree and strength,
provide other important information. The degree and strength of a vertex
represent the number and volume of a country’s exchanges, then we get insight
of the diversification and intensity of trades. Figure 1 shows the marginal
distributions of in-out degree and in-out strength. In- and out-strength are
graphically represented using the logarithmic scale with base 10.

Figure 1: Marginal distributions of the considered centrality measures.
Histogram classes represent the deciles.

The distributions of “in” and “out” measures are quite similar, but it is
worth stressing some characteristics of the degree and log-strength
distributions:

  1. 1.

on the $x-$ axis, the magnitude of the degrees corresponding to the number of
countries is few hundreds, whereas that of the strength ranges from $10^{6}$
to $10^{12}$;

  2. 2.

the log-strength distribution is bell-shaped;

  3. 3.

the degree distribution is polarized, as more than half of the countries have
connections with less than 100 other countries, and one-tenth of the countries
are connected with almost all the other ones.

As we will see later, this last characteristic is meaningful in searching for
clusters of central and peripheral countries.

We now compare countries through the four centrality measures: in-degree, out-
degree, in-log-strength and out-log-strength.

Table 1 shows the strong correlation (computed using the Pearson coefficient)
between the considered variables.

| in-deg. | out-deg. | in-log-str. | out-log-str.
---|---|---|---|---
in-deg. | 1 |  |  |
out-deg. | 0.883 | 1 |  |
in-log-str. | 0.792 | 0.785 | 1 |
out-log-str. | 0.742 | 0.771 | 0.947 | 1
Table 1: Correlation matrix between centrality measures.

As expected, “in” and “out” directions are strongly correlated ($0.883$ and
$0.947$ for degree and log-strength, respectively). However, there is also a
strong _linear_ correlation between degree and log-strength centralities with
values ranging from $0.771$ to $0.792$. Although the correlation is not
causation, these high values suggest that an increasing number of connections
increases volumes exponentially.

We now analyse the principal components from the standardized variables. The
first component represents 87% of the whole variability (considering the first
two together, the variability rises to 96%), and its loading is all positive.
Therefore, its scores can be used to get a mixed and rough order of the
countries in terms of centrality. The top 10 of the central countries are, in
increasing order, Canada, Belgium, Japan, Italy, Great Britain, Netherlands,
France, Germany, China and the United States.

We performed a hierarchical cluster analysis using the four principal
components instead of the strongly correlated original variables. We used
Euclidean distances between countries, and we defined distances between groups
with the complete method.

The number of groups suggested by the hierarchical method is two. One cluster
of 116 is composed of small territories and countries almost peripheral in
terms of volumes, and the other cluster of the remaining 106 countries
controls most of the volumes. We refer to this cluster as “main network”.

Figure 2: World map of the trade network. Countries are coloured in green and
pink, according to the cluster they belong to.

Figure 2 provides a graphic representation of the trade network, showing the
groups that emerged from the cluster analysis, namely the main network and
peripheral countries.
However, how these groups relate to each other is a question that still needs
to be investigated. As we will see in the next subsection, assortativity
provides an answer in terms of centrality correlation.

###  4.2 Higher order assortativity

Evidently, there are several combinations of centrality measures and adjacency
matrices to be used for computing the higher order assortativity coefficient.
In this section, we focus only on the cases offering a more intuitive
interpretation of the results in the context of world trade333We also computed
the higher order assortativity measures in all the remaining cases and results
are available upon request..

It is important to observe that, working on directed networks, in the
assortativity coefficient defined in formula (1) there is not symmetry between
vectors x and y. Indeed, in this application x and y are centrality measures
and, as pointed out in Section 3, the resulting coefficient depends on the out
and in-centralities we stand for x and y. We give the following economical
interpretation of the used centrality measures:

  * •

the out-going centrality measures the supply, i.e., out-degree and out-
strength of vertex $i$ is the supply of good of country $i$.

  * •

the in-going centrality measures the demand, i.e., in-degree and in-strength
of vertex $i$ is the demand of good of country $i$.

  * •

the direction of the edges represents the export, i.e., an arrow from $i$ to
$j$ represents the export from country $i$ to country $j$.

In view of the economic meaning of centrality measures given above, the most
suitable directed assortativity in this application is the out-in
assortativity $r^{OI}_{h}$ in order to effectively track the paths of goods’
exchanges between countries.

Figure 3 depicts the higher order direct assortativity for increasing $h$
computed on both the unweighted and weighted network.

Figure 3: Export-import higher order assortativity for the unweighted (left)
and weighted (right) network.

In both cases, the network is disassortative at order $h=1$; this means that
the exporters with a large number of connections and a large volume tend to
trade with small importers. Note that, in absolute value, the assortativity is
immediately lower than $0.5$, and it is due to big traders that also need to
trade one each other (for instance, The United States and China) and maybe
also small traders.
At order $h=2$ the assortativity becomes positive. A similar behaviour has
been also observed for undirected networks (see Arcagni et al., 2017, 2019),
and it may be the consequence of the disassortative behaviour at order $h=1$.
Indeed, if, at the first step, exporters and importers tend to belong to
different groups, then in step two, they tend to belong to the same one. This
bounce between disassortative and assortative behaviour continues for odd and
even orders, respectively.

Considering the weights the behaviour of higher order assortativity is similar
because of the correlations between centrality measures reported in Table 1.
However, considering weights, we observe a strong reduction in the amplitude
of variation. In fact, weights on arcs separate the main network, identified
with the cluster analysis, from the remaining, and the main network of trades
is almost complete as expected.

Finally, we observe that assortativity tends to $0$ when the order increases
because of the increment of randomness in the expected flow of goods in the
network.

##  5 A Markov chains-based interpretation of the higher order assortativity

In this section, we formalize the informative content of the higher-order
assortativity measure in terms of the transition probabilities of discrete-
time homogeneous Markov chains with finite states. As a premise, we introduce
some basic definitions and notations for the Markov chains. For a complete
treatment, the reader can refer, for instance, to Norris et al. (1998).

###  5.1 Markov chains: basic definitions and notation

A homogeneous discrete-time Markov chain with finite states is a stochastic
process $\mathcal{X}=(X(t):t\in\mathbb{N})$, such that the following
conditions are true.

  * $(P1)$

There exists a set $V$ with $n$ elements such that $X(t)\in V$, for each
$t\in\mathbb{N}$. The elements of $V$ are the states of $\mathcal{X}$, and the
set $V$ is the states space of the Markov chain.

  * $(P2)$

The Markov property holds, i.e.:

$P(X(t+1)=j|X(t)=i,X(t-1)=i_{t-1},\dots,X(0)=i_{0})=P(X(t+1)=j|X(t)=i),$

for each $t\in\mathbb{N}$ and $i,j,i_{0},\dots,i_{t-1}\in V$.

  * $(P3)$

The one-step transition probabilities are invariant with respect to time,
i.e.:

$P(X(t+1)=j|X(t)=i)=P(X(1)=j|X(0)=i),$

for each $t\in\mathbb{N}$ and $i,j\in V$. We collect all the one-step
transition probabilities in a $n$-square matrix $\mathbf{P}$, whose generic
element is $p_{ij}=P(X(1)=j|X(0)=i)$. $\mathbf{P}$ is the transition
probability matrix of the Markov chain $\mathcal{X}$.

As we will see below, the nodes set $V$ of the network coincides with the
states space of the Markov chain.

Properties $(P1)-(P3)$ allow us to identify a Markov chain through three
elements: the states space $V$, the transition probability matrix $\mathbf{P}$
and the initial probability distribution
$\pi_{0}=(\pi_{0}(1),\dots,\pi_{0}(n))$, where $\pi_{0}(i)=P(X(0)=i)$, for
each $i\in V$. In particular, condition $(P3)$ states that the Markov chain
$\mathcal{X}$ is homogeneous.

Furthermore, property $(P3)$ can be extended to the case of $h$-steps
transition probabilities that are invariant with respect to time for each
$h\in\mathbb{N}$. A straightforward computation shows that

$P(X(t+h)=j|X(t)=i)=P(X(h)=j|X(0)=i),$ (2)

for each $t,h\in\mathbb{N}$ and $i,j\in V$. We collect the $h$-steps
transition probabilities in the $n$-square matrix $\mathbf{P}_{h}$. The
element of such a matrix in the $i$-th row and $j$-th column is
$p^{(h)}_{ij}=P(X(h)=j|X(0)=i)$.

Since the states space $V$ is finite, the matrix $\mathbf{P}_{h}$ is
$\mathbf{P}^{h}$, i.e., the $h$-steps transition probability matrix is the
$h$-th power of the one-step matrix $\mathbf{P}$.

###  5.2 Higher-order assortativity in the context of Markov chains

Given $h\geq 1$, we recall that $\|\textbf{E}_{h}\|_{1}=1$; therefore, by
referring to the Markov chains theory, the generic element $e^{(h)}_{ij}$ may
represent the joint probability that $i$ and $j$ are the position states of
$\mathcal{X}$ at time $t$ and $t+h$, respectively, for each $t$, i.e.:

$e_{ij}^{(h)}=P(X(h)=j,X(0)=i).$ (3)

By considering the Markov chain $\mathcal{X}$ for which (3) is true and by
using the conditioned probability definition, we can write

$\mathbf{E}_{h}=\mathbf{D}_{\pi_{0}}\mathbf{P}^{h}$ (4)

where $\mathbf{D}_{\pi_{0}}$ is a diagonal $n$-squared matrix whose diagonal
entries are the components of the initial distribution $\pi_{0}$. Being matrix
$\mathbf{P}$, and then $\mathbf{P}^{h}$ with $h>1$, stochastic by row,
$\mathbf{P}^{h}\textbf{1}=\textbf{1}$, hence:

$\textbf{E}_{h}\textbf{1}=(\mathbf{D}_{\pi_{0}}\mathbf{P}^{h})\textbf{1}=\mathbf{D}_{\pi_{0}}(\mathbf{P}^{h}\textbf{1})=\mathbf{D}_{\pi_{0}}\textbf{1}=\pi_{0}$
(5)

then $\pi_{0}=\textbf{p}_{h}$, $h\geq 1$, where $\textbf{p}_{h}$ is the vector
appearing in formula (1). Indeed, being $e^{(h)}_{ij}$ joint probabilities, by
equation (3):

$\sum_{j=1}^{n}e_{ij}^{(h)}=\sum_{j=1}^{n}P(X(h)=j,X(0)=i)=P(X(0)=i)=\pi_{0}(i)\,\forall
i=1,...,n$

that is, $\textbf{p}_{h}$ is the initial probability distribution.

It is clear that, in general, there is not a unique Markov chain $\mathcal{X}$
with state space $V$ for which (4) is valid. However, having set
$\textbf{E}_{1}=\frac{\textbf{W}}{\|\textbf{W}\|_{1}}$, in virtue of equation
(5), there is a unique vector $\pi_{0}$. So, once we choose the network
$N=(V,A)$, there is a unique Markov chain $\mathcal{X}=(X(t):t\in\mathbb{N})$,
with $V,\mathbf{P}$ and $\pi_{0}$ fixed.

Note that for each node $i\in V$ that in the network is a sink, matrix
$\mathbf{P}$ has to be stochastic; therefore, a $1$ is added in the main
diagonal to associate to the sink an absorbing state in the Markov chain, and
$\pi_{0}(i)$ has to be $0$ to do not modify the topology of the network and
comply with equation (4) for $h=1$. On the contrary, $\pi_{0}(i)>0\,\forall
i\in V$ if $i$ is not a sink because (4) returns a zero line in
$\textbf{E}_{1}$ no matter the values in P.

Analogously to equation (5),

$\textbf{E}_{h}^{\top}\textbf{1}=(\mathbf{D}_{\pi_{0}}\mathbf{P}^{h})^{\top}\textbf{1}=(\mathbf{P}^{h})^{\top}\mathbf{D}_{\pi_{0}}\textbf{1}=(\mathbf{P}^{h})^{\top}\pi_{0}=\pi_{h}$
(6)

then $\pi_{h}=\textbf{q}_{h}$, $h\geq 1$, where $\textbf{q}_{h}$ is the vector
appearing in formula (1). Again, by eq. (3):

$\sum_{i=1}^{n}e_{ij}^{(h)}=\sum_{i=1}^{n}P(X(h)=j,X(0)=i)=P(X(h)=j)=\pi_{h}(j)\,\forall
j=1,...,n$

that is, $\textbf{q}_{h}$ is the probability distribution after $h$ steps.

Keeping this in mind, we can rewrite the higher order assortativity index
(defined in equation 1) in Markov chains context:

$r_{h}=r(\textbf{x},\textbf{y},\mathbf{D}_{\pi_{0}}\mathbf{P}^{h})=\frac{\mathbf{x}^{\top}\left(\mathbf{D}_{\pi_{0}}\mathbf{P}^{h}-\pi_{0}\pi_{h}^{\top}\right)\mathbf{y}}{\sqrt{\left[\mathbf{x}^{\top}\left(\mathbf{D}_{\pi_{0}}-\pi_{0}\pi_{0}^{\top}\right)\mathbf{x}\right]\left[\mathbf{y}^{\top}\left(\mathbf{D}_{\pi_{h}}-\pi_{h}\pi_{h}^{\top}\right)\mathbf{y}\right]},}$
(7)

Then, the higher order assortativity index of length $h$ represents the
autocorrelation at time-distance $h$ of a Markov chain, not evaluated on the
states $v_{i}\in V$ but on their transformations $\mathbf{x}$ and
$\mathbf{y}$. Indeed, in a network context, centrality measures can be the
image of a real function on the node set $V$.

###  5.3 Some illustrative examples

We present some toy examples to illustrate how the newly introduced
assortativity measure works in the context of Markov chains. To this aim, we
consider some instances of weighted networks whose related Markov chains have
a meaningful state classification. We enter the details.

We consider four networks $N_{k}=(V,\textbf{W}_{k})$, with $k=1,\dots,4$,
sharing the same set of nodes $V=\\{a,b,c,d,e\\}$, and with weighted adjacency
matrices $\textbf{W}_{k}$ defined as follows:

$\textbf{W}_{1}=\left(\begin{array}[]{ccccc}0&2&\mathbf{2}&2&0\\\
396&0&0&4&0\\\ 0&0&0&0&0\\\ 0&4&0&0&396\\\ 0&2&\mathbf{2}&2&0\\\
\end{array}\right),\qquad\textbf{W}_{2}=\left(\begin{array}[]{ccccc}0&2&\mathbf{196}&2&0\\\
396&0&0&4&0\\\ 0&0&0&0&0\\\ 0&4&0&0&396\\\ 0&2&\mathbf{196}&2&0\\\
\end{array}\right)$
$\textbf{W}_{3}=\left(\begin{array}[]{ccccc}0&\mathbf{2}&198&0&0\\\
0&0&0&100&0\\\ 150&0&0&0&150\\\ 0&200&0&0&0\\\ 200&0&200&0&0\\\
\end{array}\right),\qquad\textbf{W}_{4}=\left(\begin{array}[]{ccccc}0&\mathbf{198}&198&0&0\\\
0&0&0&100&0\\\ 150&0&0&0&150\\\ 0&200&0&0&0\\\ 200&0&200&0&0\\\
\end{array}\right)$

The networks are shown in Figure 4.

Figure 4: Graphical representations of the four illustrative examples.

By inspection of the figure, it is worth noting that $N_{1}$ and $N_{2}$ have
the same topology and differ only by two weights. In particular, there are no
outflows from $c$; that is, $c$ is a sink. However, while in the former case
the inflows in $c$ have remarkably low weights, in the latter, the inflows in
$c$ are particularly high. Such evidence is true in terms of weights’ absolute
values (i.e. considering matrices $\textbf{W}_{1}$ and $\textbf{W}_{2}$) but
also under a relative perspective (i.e. considering the transition probability
matrices $\textbf{P}_{1}$ and $\textbf{P}_{2}$ shown below).

Also, $N_{3}$ and $N_{4}$ share the same topology, and there are no outflows
from the class of nodes $\\{b,d\\}\subset V$. In both cases, it is possible to
have an inflow in such a class from the set $\\{a,c,e\\}$, through the arc
from $a$ to $b$. Such an arc has a low weight in network $N_{3}$ and
remarkably high for $N_{4}$.

We now consider the centrality measures. In Table 2, we report the in-strength
$\textbf{s}^{I}_{k}$ and out-strength $\textbf{s}^{O}_{k}$ distributions for
each network $N_{k},\,k=1,...,4$ as well the average values. Notice that the
average values for the in-strengths and the out-strengths are the same, due to
the property of associativity of the sum.

$V$ | $\textbf{s}^{I}_{1}$ | $\textbf{s}^{O}_{1}$ | $\textbf{s}^{I}_{2}$ | $\textbf{s}^{O}_{2}$ | $\textbf{s}^{I}_{3}$ | $\textbf{s}^{O}_{3}$ | $\textbf{s}^{I}_{4}$ | $\textbf{s}^{O}_{4}$
---|---|---|---|---|---|---|---|---
a | $396$ | $6$ | $396$ | $200$ | $350$ | $200$ | $350$ | $396$
b | $8$ | $400$ | $8$ | $400$ | $202$ | $100$ | $398$ | $100$
c | $4$ | $0$ | $392$ | $0$ | $398$ | $300$ | $398$ | $300$
d | $8$ | $400$ | $8$ | $400$ | $100$ | $200$ | $100$ | $200$
e | $396$ | $6$ | $396$ | $200$ | $150$ | $400$ | $150$ | $400$
mean |  | $162.4$ |  | $240$ |  | $240$ |  | $279.2$
Table 2: In-strength $\textbf{s}^{I}_{k}$ and out-strength
$\textbf{s}^{O}_{k}$ for each example $k=1,\dots,4$.

We classify the nodes as big or small with respect to the average values of
the corresponding centrality measure. From Table 2, it is possible to observe
that in network $N_{1}$, nodes $a$ and $e$ are big receivers but small
spreaders. On the contrary, $b$ and $d$ are small receivers and big spreaders
and node $c$ is peripheral. In network $N_{2}$, the node centrality is the
same as the one in network $N_{1}$, with the exception of node $c$ that
becomes a big receiver. In network $N_{3}$, nodes of the absorbing class
$\\{b,d\\}$ are peripheral, and the node $a$ is a big receiver but a small
spreader. Node $e$ is a small receiver but a big spreader, and $c$ is central.
Centralities in the network $N_{4}$ differ only for node $a$ that also becomes
a big spreader, consequentially a central node. Hence $b$ is affected by such
variation and becomes a big receiver.

The homogeneous Markov chains related to the networks of the examples share
the same state space $V=\\{a,b,c,d,e\\}$. We denote the Markov chain
associated to the network $N_{k}$ by $\mathcal{X}_{k}$, for $k=1,\dots,4$. The
transition matrix and the initial probability distribution of
$\mathcal{X}_{k}$ can be obtained by the adjacency matrix $\textbf{W}_{k}$ and
will be labelled by $\textbf{P}_{k}$ and $\pi_{0,k}$. We here list the
transition matrices.

$\textbf{P}_{1}=\left(\begin{array}[]{lllll}0&0.\bar{3}&0.\bar{3}&0.\bar{3}&0\\\
0.99&0&0&0.01&0\\\ 0&0&1&0&0\\\ 0&0.01&0&0&0.99\\\
0&0.\bar{3}&0.\bar{3}&0.\bar{3}&0\\\
\end{array}\right),\qquad\textbf{P}_{2}=\left(\begin{array}[]{lllll}0&0.01&0.98&0.01&0\\\
0.99&0&0&0.01&0\\\ 0&0&1&0&0\\\ 0&0.01&0&0&0.99\\\ 0&0.01&0.98&0.01&0\\\
\end{array}\right)$
$\textbf{P}_{3}=\left(\begin{array}[]{lllll}0&0.01&0.99&0&0\\\ 0&0&0&1&0\\\
0.5&0&0&0&0.5\\\ 0&1&0&0&0\\\ 0.5&0&0.5&0&0\\\
\end{array}\right),\qquad\textbf{P}_{4}=\left(\begin{array}[]{lllll}0&0.5&0.5&0&0\\\
0&0&0&1&0\\\ 0.5&0&0&0&0.5\\\ 0&1&0&0&0\\\ 0.5&0&0.5&0&0\\\
\end{array}\right)$

The initial probability distributions are proportional to the out-strengths
because of Equation (5). They can be written as follows:

$\displaystyle\pi_{0,1}^{\top}$ $\displaystyle=$
$\displaystyle[0.008,0.493,0,0.493,0.007]$ $\displaystyle\pi_{0,2}^{\top}$
$\displaystyle=$ $\displaystyle[0.1\bar{6},0.\bar{3},0,0.\bar{3},0.1\bar{6}]$
$\displaystyle\pi_{0,3}^{\top}$ $\displaystyle=$
$\displaystyle[0.1\bar{6},0.08\bar{3},0.25,0.1\bar{6},0.\bar{3}]$
$\displaystyle\pi_{0,4}^{\top}$ $\displaystyle=$
$\displaystyle[0.284,0.072,0.215,0.143,0.287].$

The states of the Markov chains $\mathcal{X}$ can be easily classified. $c$ is
an absorbing state for $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$ and the set of
states $\\{b,d\\}$ is an absorbing class for the Markov chains
$\mathcal{X}_{3}$ and $\mathcal{X}_{4}$. In accordance with the related
networks, the probabilities of being absorbed by $c$ and $\\{b,d\\}$ represent
the elements distinguishing $\mathcal{X}_{1}$ from $\mathcal{X}_{2}$ and
$\mathcal{X}_{3}$ from $\mathcal{X}_{4}$, respectively.

Now we will inspect the assortativity as the correlation between the out-
strength and the in-strength associated with the states of the mentioned
Markov Chains.

$h$ | $N_{1}$ | $N_{2}$ | $N_{3}$ | $N_{4}$
---|---|---|---|---
$1$ | $0.7718$ | $0.0828$ | $0.5404$ | $0.5503$
$2$ | $-0.7664$ | $-0.0739$ | $0.5204$ | $-0.0309$
$3$ | $0.6031$ | $0.0507$ | $0.3515$ | $0.2663$
$4$ | $-0.6031$ | $-0.0470$ | $0.5067$ | $-0.0130$
$5$ | $0.5086$ | $0.0361$ | $0.4034$ | $0.1688$
$6$ | $-0.5550$ | $-0.0339$ | $0.4551$ | $-0.0236$
$\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$
$9995$ | $0$ | $0$ | $0.0098$ | $0.0718$
$9996$ | $0$ | $0$ | $-0.0098$ | $-0.0718$
$9997$ | $0$ | $0$ | $0.0098$ | $0.0718$
$9998$ | $0$ | $0$ | $-0.0098$ | $-0.0718$
$9999$ | $0$ | $0$ | $0.0098$ | $0.0718$
$10000$ | $0$ | $0$ | $-0.0098$ | $-0.0718$
Table 3: First and last evaluated values of $r^{OI}_{h}$ of the four
illustrative examples. Figure 5: $r^{OI}_{h}$ for increasing values of $h$ of
the networks $N_{1}$ and $N_{2}$.

Figure 5 reports the higher order assortativity out-in, $r^{OI}_{h}$, for the
networks $N_{1}$ and $N_{2}$. In $N_{1}$ the probability of falling in the
absorbing state $c$ is small; therefore, we can consider its effect as
marginal. Nodes $a$ and $e$ have a small out-strength and high in-strength and
in contrast to $b$ and $d$, which have opposite characteristics. We identify
the pair $\\{a,e\\}$ as group 1 and the pair $\\{b,d\\}$ as group 2. At the
first step, $h=1$, the assortativity is due mainly to the correlation between
the out-strengths in group 1 with the in-strengths in group 2, and vice-versa.
It can be observed that there is concordance between them and, therefore,
positive assortativity. After two steps, $h=2$, if the Markov Chain starts
from group 1, it is highly probable that it jumps to group 2 and back to group
1, generating a disassortative behaviour. It is the same if the Markov Chain
starts from group 2. As a consequence, for increasing values of $h$, this
creates the alternation between positive and negative assortativity for odd
and even values of $h$, respectively, as confirmed in Table 3. The asymptotic
zero-assortative behaviour is caused by both the randomness and the presence
of the absorbing state $c$. In network $N_{2}$, it can be observed the same
behaviour with a smaller oscillation. This is due to the increased out-degrees
of nodes in group 1 and to the higher probability of falling in the absorbing
state $c$ where there is no centrality variability.

Observe that, even if in the sink there is no variability of strength
measures, there is not an indeterminate-form in Equation (1) because the
probability of falling in the sink is less than one. Moreover in the network
$N_{1}$, the probability of falling in the sink is high; therefore, the
asymptotic value is reached immediately. On the contrary, the asymptotic value
is reached after a long alternation between assortative and disassortative
behaviour in the network $N_{2}$.

Figure 6: $r^{OI}_{h}$ for increasing values of $h$ of the networks $N_{3}$
and $N_{4}$.

Figure 6 shows the same results for networks $N_{3}$ and $N_{4}$. Such
networks are characterized by the absorbing class $\\{b,d\\}$ and the
complementary one $\\{a,c,e\\}$ connected by the bridge represented by the
edge from $a$ to $b$. Both classes have positive assortativity. In $N_{3}$,
the bridge has a low weight; therefore, the inter-assortative behaviour
dominates. In network $N_{4}$, the weight of the bridge is very high, and
consequentially the larger oscillations are justified by the differences
between the classes.

##  6 Conclusions

This paper fills a gap in the complex network literature by introducing a
concept of higher order assortativity measure for weighted and directed
networks, where the considered nodes’ attributes are the directed version of
degree and strength centrality. Such measures have a clear interpretation in
application fields well-linked to direct and weighted networks, such as that
of the international trade, studied in the proposed empirical application.
Furthermore, we also highlight the strong connection of the introduced
assortativity measure and the autocorrelations of suitably defined Markov
chains.

The versatility of the methodological proposal makes it possible to describe
the preferential attachment of a wide set of models. Some paradigmatic
empirical instances might involve the analysis of the roots of the migration
flows (Cerqueti et al., 2019) or the structure and systemic risk profile of
the interbanking system (Bo and Capponi, 2015; Castellano et al., 2021;
Castiglionesi and Eboli, 2018; Cerqueti et al., 2022, 2021; Bargigli et al.,
2015). Moreover, the readability of the proposed assortativity measure in the
context of Markov chains might be efficiently exploited to analyse the
properties of some classes of dynamical random systems having the Markovian
properties. In this respect, it is important to notice that the decay of the
autocorrelation explains the long-term memory properties of the underlying
stochastic process (see Hurst, 1951). Such challenging themes are already in
our research agenda.

## References

  * Antonietti et al. (2022) Antonietti, R., Falbo, P., Fontini, F., Grassi, R., and Rizzini, G. (2022).  The world trade network: country centrality and the covid-19 pandemic.  Applied Network Science, 7(1):18.
  * Arcagni et al. (2017) Arcagni, A., Grassi, R., Stefani, S., and Torriero, A. (2017).  Higher order assortativity in complex networks.  European Journal of Operational Research, 262(2):708–719.
  * Arcagni et al. (2019) Arcagni, A., Grassi, R., Stefani, S., and Torriero, A. (2019).  Extending assortativity: An application to weighted social networks.  Journal of Business Research.
  * Bang-Jensen and Gutin (2008) Bang-Jensen, J. and Gutin, G. Z. (2008).  Digraphs: theory, algorithms and applications.  Springer Science & Business Media.
  * Bargigli et al. (2015) Bargigli, L., Di Iasio, G., Infante, L., Lillo, F., and Pierobon, F. (2015).  The multiplex structure of interbank networks.  Quantitative Finance, 15(4):673–691.
  * Bo and Capponi (2015) Bo, L. and Capponi, A. (2015).  Systemic risk in interbanking networks.  SIAM Journal on Financial Mathematics, 6(1):386–424.
  * Castellano et al. (2021) Castellano, R., Cerqueti, R., Clemente, G. P., and Grassi, R. (2021).  An optimization model for minimizing systemic risk.  Mathematics and Financial Economics, 15:103–129.
  * Castiglionesi and Eboli (2018) Castiglionesi, F. and Eboli, M. (2018).  Liquidity flows in interbank networks.  Review of Finance, 22(4):1291–1334.
  * Cerqueti et al. (2022) Cerqueti, R., Cinelli, M., Ferraro, G., and Iovanella, A. (2022).  Financial interbanking networks resilience under shocks propagation.  Annals of Operations Research, https://doi.org/10.1007/s10479-022-04567-w.
  * Cerqueti et al. (2019) Cerqueti, R., Clemente, G. P., and Grassi, R. (2019).  A network-based measure of the socio-economic roots of the migration flows.  Social Indicators Research, 146:187–204.
  * Cerqueti et al. (2021) Cerqueti, R., Clemente, G. P., and Grassi, R. (2021).  Systemic risk assessment through high order clustering coefficient.  Annals of Operations Research, 299:1165–1187.
  * Estrada et al. (2008) Estrada, E., Hatano, N., and Gutierrez, A. (2008).  “Clumpiness” mixing in complex networks.  Journal of Statistical Mechanics: Theory and Experiment, page P03008.
  * Foster et al. (2010) Foster, J. G., Foster, D. V., Grassberger, P., and Paczuski, M. (2010).  Edge direction and the structure of networks.  Proceedings of the National Academy of Sciences, 107(24):10815–10820.
  * Gómez et al. (2010) Gómez, S., Arenas, A., Borge-Holthoefer, J., Meloni, S., and Moreno, Y. (2010).  Discrete-time markov chain approach to contact-based disease spreading in complex networks.  EPL (Europhysics Letters), 89(3):38009.
  * Gross and Yellen (2003) Gross, J. L. and Yellen, J. (2003).  Handbook of graph theory.  CRC press.
  * Hurst (1951) Hurst, H. E. (1951).  Long-term storage capacity of reservoirs.  Transactions of the American society of civil engineers, 116(1):770–799.
  * Iannelli et al. (2017) Iannelli, F., Koher, A., Brockmann, D., Hövel, P., and Sokolov, I. M. (2017).  Effective distances for epidemics spreading on complex networks.  Physical Review E, 95(1):012313.
  * Leung and Chau (2007) Leung, C. and Chau, H. (2007).  Weighted assortative and disassortative networks model.  Physica A: Statistical Mechanics and its Applications, 378(2):591–602.
  * Mayo et al. (2015) Mayo, M., Abdelzaher, A., and Ghosh, P. (2015).  Long-range degree correlations in complex networks.  Computational Social Networks, 2(1):1–13.
  * Meghanathan (2016) Meghanathan, N. (2016).  Assortativity analysis of real-world network graphs based on centrality metrics.  Computer and Information Science, 9(3):7–25.
  * Newman (2010) Newman, M. (2010).  Networks: an introduction.  Oxford university press.
  * Newman (2002) Newman, M. E. J. (2002).  Assortative Mixing in Networks.  Physical Review Letters, 89(20):208701.
  * Norris et al. (1998) Norris, J. R., Norris, J. R., and Norris, J. R. (1998).  Markov chains.  Number 2 in Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Piraveenan et al. (2010) Piraveenan, M., Prokopenko, M., and Zomaya, A. (2010).  Assortative mixing in directed biological networks.  IEEE/ACM Transactions on computational biology and bioinformatics, 9(1):66–78.
  * UN COMTRADE (2022) UN COMTRADE (2022).  International trade statistics database.  data retrieved from International Trade Statistics Database, https://comtrade.un.org/.
  * Van Mieghem et al. (2010) Van Mieghem, P., Wang, H., Ge, X., Tang, S., and Kuipers, F. A. (2010).  Influence of assortativity and degree-preserving rewiring on the spectra of networks.  The European Physical Journal B, 76(4):643–652.
  * Yuan et al. (2021) Yuan, Y., Yan, J., and Zhang, P. (2021).  Assortativity measures for weighted and directed networks.  Journal of Complex Networks, 9(2):cnab017.


# A laser projection system for polymer transverse strain measurements in
tensile testing

Dawid Kucharski  Division of Metrology and Measurement Systems, Institute of
Mechanical Technology  Maciej Obst  Division of Strength of Materials and
Structures, Institute of Applied Mechanics  Jarosław Adamiec  Institute of
Machine Design

###### Abstract

A self-built laser projection system for optical transverse strain
measurements in tensile testing is presented. The setup based on laser diode
modules, CCD (Charge-Coupled Device) array detectors, and standard
optomechanical components is proposed to be a low-cost system for testing
polymer samples. The optical setup is mobile and can easily be mounted on a
tensile test machine. Circular cross-sectional changes in polymer diameter
were detected with ten-micrometre resolution and stress-strain curves obtained
by a tensile test machine.
The article describes signal processing using the open-source programming
language R, and the transverse deformation of the polymer sample is evaluated
in two perpendicular planes. The results obtained by the optical system may be
used to accurately describe the models of energy-based mechanical properties
of polymers for complex load conditions.

Keywords: polymers mechanical properties, projection system, optical
metrology, tensile testing, laser system, R programming.

##  1 Introduction

The stress-strain characteristics of materials supplemented with accurate
transverse deformation measurements are essential for a more precise
description of the material’s mechanical properties. A material is considered
non-classical if it violates, or at least does not fully satisfy, one or more
of the assumptions on which classical continuum mechanics are based. An
example is a deviation from Hooke’s law for polymers (elastic nonlinearity,
etc.) opposite to the metals as the classical materials. Especially modern
non-classical materials require the development of appropriate mechanical
properties and accurate approximation in a wide range of strains.
The stress-strain curves of metallic materials do not always provide a clear
mechanical description of the sample; the effect may even be more present with
polymers. As S. Tu et al. [1] pointed out for metals, the flow stress in the
post-neck region is essential to evaluate the hardening of the metal sheet.
The effect is not observed in polymers, but accurate detection of the necking
point is crucial for the mechanical evaluation of polymers during uniaxial
tensile tests, for example, for the better plastic instability creep and
relaxation evaluation process.
C. G’sell and J. Jones in [2] introduced a relatively simple circumferential
gauge specially designed to measure the neck diameters of polymer specimens.
They determined the flow curves of polyvinyl chloride (PVC) and high-density
polyethene (HDPE). More precise techniques are based on optics. F. Addiego et
al. used a video system to accurately determine axial strain using fluorescent
markers illuminated with an ultraviolet lamp [3]. A typical optical technique
is extensometry, also introduced by G’Sell et al. [4]. The method based on
displacement measurements was applied for the description of the mechanical
behaviour of PTFE (polytetrafluoroethylene) by L. Nunes et al. [5]. The method
has some severe limitations. Markers for video extensometry can modify the
mechanical response of the tested materials. S. Todros et al. introduced the
problem in biomedical polymers [6].
An exciting approach to optical methods for the mechanical evaluation of
cylindrical PMMA shells (polymethylmethacrylate) was introduced by W. Liu and
X. Zhang [7]. They interposed the samples in monochromatic light from a green
laser with power $P=50$ mW.
S. Muhammad and P. Jar postulated that in HDPEs, necking in thinner specimens
generated a more significant reduction in thickness direction but lower in
width direction [8]. The data obtained by the extensometer suggested that
resistance to neck propagation was highest for the thinnest specimens, as
their neck growth speed was lowest and the flow stress highest. The authors
used the data for the FE simulation to mimic the necking process. Our proposed
bi-axial optical setup might verify these results with direct strain
measurements.

In the Industry 4.0 age (The Fourth Industrial Revolution), where the
development of smart factories and autonomous systems is now widely observed,
optical instruments are game-changers in metrology and mechanical engineering
[9]. In addition to surface measurements [10, 11], optics also provides sub-
micrometre resolution and high precision in the evaluation of mechanical
properties. L. Zhang et al. proposed a simple laser setup to monitor the
necking process during a material tensile test [12]. The setup was not finally
constructed but presented a promising application also for polymers.
The projection technique is another group of optical methods for specimen
geometry analysis during tensile testing. R. Srinivasan et al. presented the
projection method based on the modified shadow moire technique to measure
large and small structures. They discussed a single projector, a single
exposure method for neck development detection in tensile testing [13]. H.
Hoffmann and C. Vogl used an optical projection scanner based on structured
light to measure the specimen’s geometry in the necking area [14].
An exciting study was presented by J. Ye et al. [15] on the plastic
instability of HDPE deformed under tension. Using full-field 3D-DIC (Digital
Image Correlation) and extensometry, they pointed out that a video-
extensometer based on marker tracking cannot efficiently measure the strain
rate in the specimen centre at significant strain levels when the markers are
very deformed. This limitation can be easily overcome with our proposed no-
marker setup. The most reliable way to study the rheological properties of
materials is by direct measurement of the neck profile.
An interesting optical setup for measuring edge position was described by J.
Fisher and T. Radil [16]. They compared a point light source shadow projection
system with commonly used methods. They pointed out that when the collimated
light beam illuminated the measured object, the method did not need complex
signal processing, and the measured object’s edge position was directly equal
to the position of the edge in the illumination profile. The resolution of the
presented position was approximately $1$ $\mu m$. No practical application was
shown.
A reasonably innovative method to visualise the initiation and propagation of
the neck of polymer films during tensile deformation was proposed by S. Kato
et al. [17]. They used chemical modifications of the polymer by
tetraarylsuccinonitrile (TASN) moieties and observed yellow fluorescence with
electron paramagnetic resonance (EPR) spectroscopy. They evaluated the
microforces generated in the amorphous regions during the uniaxial stretching
of the crystalline polymers.
W. Ren et al. introduced a method to measure the response to the strain of
electroactive polymers by Doppler laser interferometry [18]. They investigated
Maxwell stress actuators made of silicone and thermoplastic polyurethane.

Polymers are widely used in many applications where the correct evaluation of
mechanical properties is crucial. An exciting area is polymer-based optical
sensing [19]. A. Leal-Junior et al. [20] described a temperature sensor based
on polymer optical fibre (POF) where the fibre stress changes with temperature
variation. The authors characterised the variation of the POF mechanical
property using DMA (dynamic mechanical analysis). They proposed an analytical
model and simulation for the evaluation of sensor behaviour. Validation of the
model in temperature tests with the fibre sensor was provided. In [21] A.
Leal-Junior et al. proposed optical fibre sensing based on the dynamic
mechanical analysis (DMA) of the polymer to obtain the Young modulus
concerning the variation of strain, temperature, humidity and frequency.
Polymer optical fibre sensors can also be used for transverse force
measurements. A. Leal-Junior et al. [22] described an interesting experimental
setup for the characterisation of the strain and transverse force of fused and
non-fused POFs.
Another quite exciting polymer application was presented by A. Leal-Junior et
al. [23] about 3D printed elements to obtain low-cost and custom-made sensing
elements.
Polymer-based additive manufacturing also requires precise measurement tools
to evaluate the mechanical properties of components. The self-built optical
transverse strain measurement system proposed in the paper for the tensile
testing of polymers might be a recommended support tool.

The paper proposes a simple, low-cost optical setup for fast and precise
radius monitoring of polymer sample rods during a tensile test. The setup
based on the bi-axial laser line projection system is proposed for the first
time in the literature. The proposed construction is based on parallel beams
and linear CCD (Charge-Coupled Device) detectors, where the position of the
edge in the illumination profile is equal to the position of the measured
object’s edge. The metrology evaluation of the setup followed by the
calibration of the steel rods is introduced. A self-prepared algorithm based
on the R programming language [24] was used for signal processing and data
evaluation. This paper presents tests of the optical setup on a tensile
machine. The system is proposed as an additional measurement tool for the
experimental verification of the mechanical properties models [25].
In the paper, the following are described:

  * •

the design and construction of the new optical system for fast and precise
transverse strain monitoring of polymers in tensile testing;

  * •

the fast open-source R programme developed for signal processing. The code is
freely available on GitHub;

  * •

metrology calibration estimation of the system;

  * •

detection of the polymer sample size changes during the tensile tests on the
machine for various types of polymer rods;

  * •

different mechanical behaviour of the materials depends on the type and sample
preparation;

  * •

complex symmetrical and non-symmetrical mechanical effects observed by the
sample radius monitoring;

  * •

the importance of the direct sample strain measurements under the tensile
testing or compressing, which the proposed low-cost system can achieve;

  * •

the importance of further works with the setup to revise non-classical
material mechanical properties models.

##  2 Experimental details

###  2.1 Setup

The measurement system consists of two semiconductor lasers ($\lambda=650$ nm,
$P=5$ mW), two CCD array detectors (TCD1304), two Nucleo boards (F401RE) and a
self-built laser power supply. The optical setup is constructed with off-the-
shelf optomechanical components. The optical layout of the system is shown in
fig. 1.

Figure 1: Optical layout: $\lambda_{1,2}$ \- semiconductor lasers
($\lambda=650$ nm, $P=5$ mW), $L_{1,3}$ \- collimation lenses, $L_{2,4}$ \-
imaging lenses, $PH$ \- pinhole, $D_{1,2}$ \- CCD array detectors, $Cage$ \-
30 mm cage cube Figure 2: A perspective view of the setup

The light from the lasers is collimated into horizontal line-shaped beams by
the lenses $L_{1,3}$. A sample is irradiated from two perpendicular directions
for bi-axial radius detection. Two CCD array detectors ($D_{1}$ and $D_{2}$)
register the shadow image during a test. In this way, real-time monitoring of
the development of the sample neck can be performed. A self-built stabilised
power supply ($I_{max}=150$ $A$; $V_{max}=30$ $V$; $V_{set}=4$ $V$) powered
economic, non-stabilised lasers (no wavelength stabilisation). The lasers were
tested for light-power constancy. Long-term linear drift ($\sim 6$ hours) and
periodic oscillations were observed (see fig. 3(a) and 3(b)). A moving average
filter was used for the initial smoothing of the data. Linear and periodic
detrending (removing a trend by model fitting) methods were implemented for
optional final sample measurement results compensation if needed for time-
consuming sample testing. To eliminate the linear trend, the detrend function
of the pracma R package was applied [26]. For the elimination of periodic
oscillations, a non-linear fitting function was used from the stats R package
[24]. As shown in table 1, the power instabilities of the lasers are
approximately $0.30$ $\mu W$ and after trend reductions, less than $0.15$ $\mu
W$ (estimated as the maximum amplitude.
The power instabilities are negligible for the application presented in the
paper, for short-time measurements (approximately $10$ min. per sample), but
it may be important for more advanced measurements for a time-consuming sub-
micrometre application. The results presented in the paper with short-time
measurements also avoid the second necking appearance problem in polymers
tensile testing [27]. However, the laser parameters should be tested
experimentally prior to use. It is crucial, especially for low-cost lasers
($\leq 1\pounds$ per diode), that might also be applied for an advanced
purpose.

Table 1: Laser power instability estimations. $A$ – Amplitude, $\omega$ – frequency, $t$ – time, $\phi$ – phase, $C$ – constant. The data correspond to fig. 3 Detrending | Function |  $\Delta P_{D_{1}}$ $[\mu W]$ |  $\Delta P_{D_{2}}$ $[\mu W]$
---|---|---|---
non | - | $0.27$ | $0.31$
linear | $y=a\cdot t+b$ | $0.15$ | $0.24$
periodic | $y=A\cdot\cos(\omega\cdot t+\phi)+C$ | $0.06$ | $0.13$

(a)

(b)

(c)

(d)

Figure 3: Lasers light power stability: (a) - oscillations of the $\\#1$ laser
diode light power in time (grey dotted line: pure data, black line: average
filtered data, dashed line: linear regression), (b) - oscillations of the
$\\#2$ laser diode light power in time (grey dotted line: pure data, black
line: average filtered data, dashed line: linear regression), (c) - linearly
detrended oscillations of the $\\#1$ laser diode light power in time (black
line: average filtered data, dotted line: cosine fitting curve), (d) -
linearly detrended oscillations of the $\\#2$ laser diode light power in time
(black line: average filtered data, dotted line: cosine fitting curve)

The TCD1304 is a high-sensitivity 3648-pixel linear CCD. Two CCD detectors
were applied in the setup. Two Nucleo boards (F401RE) were used for data
registering, controlled by a command-line interface (CLI). The CCDs are
analogue devices that collect photons and convert them into charge. The
collected charge is then shifted to the output, where an analogue-to-digital
converter (ADC) digitises it. The ADC signal is cleaner than the raw output
when a capture is done with the ADC of the microcontroller. There is a linear
correlation between the input voltage and the converted value (ADC). The
STM32F401RE ADC range is $0-3.3V$. The ADC has a 12-bit resolution, so an
input of $3.3V$ corresponds to a value of $2^{12}=4096$, and $0V$ corresponds
to a value of $0$. Details about the sensors and software used can be found in
[28].

###  2.2 Signal processing with R

R is a programming language and environment for statistical computing and
graphics. Developed at Bell Laboratories, R is influential in statistics and
image processing and is highly extensible via packages [24].
An open-source R algorithm was developed to process CCD signals and to
evaluate specimen dimensions. The code is freely available on GitHub [29].
Open-source C programmes were used to control Nucleo boards and capture CCD
data [28]. The programmes were run parallel for each board using simple Unix
Bash shell commands. The signals were stored as *.dat files for further
processing. For signal denoising, the moving average filter was applied (see,
e.g., figs. 4(a), 4(b)). Next, for signal fitting, the b-spline was
implemented using the R stats package with defined degrees of freedom [30, 31,
32]. The detection of the edge of the signal is based on the first derivative
of the fitted signal. In this way, the width of the shadow is determined as
the parameter $P$ $[$Pixels$]$. The programme steps are described in detail in
appendix A. The b-spline is more described in the appendix.
Signal processing is simple, fast, and reliable. The computation time for the
$100$ files took less than $30$ s on the Unix machine based on the Intel Core
i5-5250U $1.6$ GHz (without parallel computing, Geekbench score is 5094). The
electronic noise on the graphs 4(a) and 4(b) negligibly influences the
precision of measurement, estimated as the standard deviation of
$\overline{P}$ (average $P$) under $100$ repeated detections (with speed one
detection per $10$ seconds) at the level of $2\sigma<5$ pixels.

(a)

(b)

Figure 4: The output signals from the CCDs through the analogue-to-digital
converters (ADCs) of the Nucleo microcontroller boards. Signals depicted
(solid grey lines) for both CCDs: the first channel $\\#1$ (a) and the second
$\\#2$ (b) for a measurement of a steel rod of $\phi=9.990$ mm. Solid black
lines are filtered signals (with an average rectangle moving window of 80
pixels), and dotted lines are fitted b-splines (with degrees of freedom of
27). The two-dashed vertical lines are the cutoff corresponding to the
object’s shadow detected with width $P$ $[$Pixels$]$

###  2.3 Calibration

An optical instrument needs to be calibrated. Steel rod standards with a minor
diameter step of $10$ $\mu m$ and tolerance $1$ $\mu m$ were used to evaluate
the resolution of the system. Each rod size was detected $100$ times with a
detection speed of one per $10$ seconds for the average shadow width
($\overline{P}$ $[$Pixels$]$) and standard deviation ($2\sigma$) estimations.
The linear fitting-based graphical method was used to find the mathematical
relation between the object’s size ($\phi$) and the optical output of the
system (shadow width $P$). It is evident that for $\phi=0$, the shadow width
(theoretically) should be equal to $0$). The fitting with an intercept equal
to $0$ was checked, and the lower precision appeared than the non-zero
intercept fitting, mainly affecting the fitting in the diameter range
$\phi=10$ mm (for a comparison example, see the appendix). The reason comes
from the detector noise as well. The estimated maths functions better
represent the optical output with a free non-zero intercept than equal $0$,
and the conditions were also kept because of the tested sample sizes
$\phi=7-12$ mm. Moreover, a slight difference in the calibration equation may
affect only the estimated diameter values of the sample, while the profile of
the monitored sample diameter will remain the same.

The calibration graphs of the two setup channels are depicted in fig. 5. The
$1$ $mm$ of the diameter step was chosen for the clarity of the graphs. The
linear functions are used to monitor the sample diameter during tensile tests.
An example of a sub-millimetre optical response is shown in table 2.

(a)

(b)

Figure 5: Calibration curves (linear regression $y=ax+b=\overline{P}=a\phi+b$)
for: (a) – channel $\\#1$ with fitted parameters: $a=221.96$, $b=141$, std.
errors $\Delta a=2.47$, $\Delta b=15.36$, $p<0.001$; (b) – channel $\\#2$ with
fitted parameters: $a=228.36$, $b=92.77$, std. errors $\Delta a=5.12$, $\Delta
b=31.77$, $p<0.001$ Table 2: Example of measured calibration steel rods with
$10$ $\mu m$ diameter differences. $\\#1$ – channel 1, $\\#2$ – channel 2,
$\phi$ – diameters in $[$mm$]$, $\overline{P}$ – average shadow’s width in
$[$Pixels$]$, $2\sigma$ – $2\times$standard deviation in $[$Pixels$]$.

$\overline{P}\pm 2\sigma$ $[$Pixels$]$ |
---|---
$\phi$ $[$mm$]$ | $9.99$ | $9.98$ | $8.01$ | $8.00$ |
$\\#1$ | $2334.05\pm 1.60$ | $2330.10\pm 1.20$ | $1966.61\pm 1.61$ | $1946.58\pm 1.22$ |
$\\#2$ | $2364.59\pm 1.18$ | $2352.96\pm 1.51$ | $1896.39\pm 4.70$ | $1879.00\pm 1.28$ |

$\overline{P}\pm 2\sigma$ $[$Pixels$]$
---
$\phi$ $[$mm$]$ | $7.01$ | $7.00$
$\\#1$ | $1687.00\pm 0.70$ | $1673.88\pm 1.72$
$\\#2$ | $1694.48\pm 1.47$ | $1668.59\pm 1.18$

###  2.4 Samples and experimental workflow

The experiment scheme is shown in fig. 6. PET (polyethylene terephthalate),
PVC (polyvinyl chloride) and PVDF (polyvinylidene fluoride) rods with diameter
$\phi\approx 12$ mm were used as samples. Pure and precut rods (with
approximately $\Delta\phi=2$ mm initiation of the neck region) were tested
(see fig. 7).

Engineering stress $\sigma$ is the applied load $F$ measured by the machine
sensor, divided by the original cross-sectional area of material $A_{0}$,
calculated from the initial diameter of a sample (see eq. 1).

$\sigma=\frac{F}{A_{0}}=\frac{F}{\pi r^{2}},$ (1)

where: $F$ – applied force, $A_{0}$ – original cross-sectional area, $r$ –
sample radius.
The strain is the amount of sample deformation due to an applied force. The
strain ($\epsilon$) is dimensionless, defined as the fractional change in
length (see eq.2).

$\epsilon=\frac{\Delta L}{L_{0}},$ (2)

where: $\Delta L$ – lenght change, $L_{0}$ – initial length.
In the results presented, $L_{0}$ corresponded to the initial position of the
upper mounting bracket of the machine without load and $\Delta L$ was the
displacement of the mounting bracket detected by the machine during the test.

Figure 6: The experimental scheme: $\lambda$’s - laser light sources, $CCD$s -
array detectors, $F$ \- elongation force
Figure 7: Sample dimensions: $\Phi\approx 12$ mm, $\Phi_{c}\approx 10$ mm,
$L\approx 20$ cm, $L_{c}\approx 3$ cm

The Zwick Z100 machine was used for tensile tests with the optical projector
attached to the moving traverse beam by a specially designed stiff arm (see
fig. 8). The maximum load is $10$ $kN$. The fixing point is in the centre of
the projector (the centre cube). During data acquisition, the optical setup on
the column was moved at half the elongation speed by the transmission belt.
For the first tests, it was $2.5$ mm/min. This speed kept the laser beams
close to the possible necking region. The reduced speed avoids the second
necking initiation in the sample and is chosen for the setup performance
tests. It is worth emphasising that the optical setup was attached to the
machine, replacing an extensometer, and preventing in this way, the ability to
simultaneously monitor the elongation of the sample in the area of the optical
projector detection. The approach to the correlation between displacement and
necking is shown, and a different statistical approach will be proposed in the
future to analyse the mechanical properties obtained with an extensometer.

(a)

(b)

Figure 8: The projector on the tensile machine: (a) front view, (b) side view

##  3 Results and discussion

Eight measurements in total were analysed for the first setup tests. The
results proved the remarkable reliability of the constructed setup. The
example results shown in the paper proved the complexity, symmetry, and non-
symmetry of the polymer’s mechanical behaviour. This corresponds to the two
measurement directions for the sample’s diameters (to the symmetric or non-
symmetric response of the channels). In fig. 9, the measurement results for
the precut PET rod are presented. The sample diameter profiles in both
transverse directions are similar and correspond to the detection in the neck
region, where the sample finally breaks. The b-splines with degrees of freedom
of $15$ were applied to smooth the profiles (for an additional description,
see the appendix). In fig. 9(c), the maximum of the stress-strain curve named
the ”necking point” corresponds to the fast radius reduction of the sample
after approximately $130$ seconds of elongation. The initial force $F_{max}$
corresponding to the maximum stress was approximately $5.2$ $kN$. The data
were stored every second for maximum precision. No time synchronisation
between the optics and the machine was prepared for the first tests. Research
with a specially prepared programme will begin presently.

(a)

(b)

(c)

Figure 9: Example results for a precut PET polymer rod: (a) - diameter change
registered by channel $\\#1$, (b) - diameter change registered by channel
$\\#2$, (c) - stress-strain curve obtained by the tensile machine. The initial
force $F_{max}\approx 5.2$ $kN$. The grey points correspond to registered
data; the dashed line corresponds to spline fitting

The results obtained show some asymmetry for the same polymer PET without pre-
necking preparation (precut) (see figs. 10(a) and 10(b)). The initial force
corresponding to the maximum stress was $F_{max}\approx 8.5$ $kN$. At the end
of the test, the diameter growth is presented. This group of points was
detected after the sample cracked. The rest of the sample still created the
shadow. This is related to the necking region shift initiated off the optical
axis and the sample shift after the crack. The effect will be further
investigated with time synchronisation.
Moreover, some fluctuations are observed during the first $50$ seconds of the
process. The effect corresponds to the machine pre-set procedure. Before the
direct procedure, the machine moved up and down (with a force of $0N-50N-0N$)
to avoid slippage of the polymer by increasing the hold of the clamps by
biting into the material. This affects probe radius oscillation and is
observed in the much smaller range for the probe with the necking precut (see
figs. 9(a) and 9(b)). This procedure does not significantly affect the
sample’s mechanical properties in this range. It is not detected even by the
machine sensors with such a small stress range. The effect might be considered
for the subsequent modelling properties using the proposed optical
measurements.

(a)

(b)

(c)

Figure 10: Example results for a solid PET polymer rod: (a) - diameter change
registered by channel $\\#1$, (b) - diameter change registered by channel
$\\#2$, (c) - stress-strain curve obtained by the tensile machine. The initial
force $F_{max}\approx 8.5$ $kN$. The grey points correspond to registered
data; the dashed line corresponds to spline fitting

An apparent difference between stress-strain curves shows sharper necking
initiation in the precut sample (see fig. 9(c)) than in the solid one (see
fig. 10(c)).
An interesting effect was observed during the PVDF sample measurements. The
optics detected the initiation of the necking in the precut sample (see Figs.
11(a) and 11(b)) and the absence of the necking process in the solid polymer
rod (see fig. 11(c) and 11(d)). The effect might be related to the homogeneity
of the samples and the repeatability of the measurements. This will be
considered for further investigation and metrology evaluation of the setup.
The growth of the sample diameter at the end of the test for the precut PVDF
is again related to the points detected after the sample breaks and shifts.
The initial force reached $F_{max}\approx 4.6$ $kN$ for the precut sample and
$F_{max}\approx 5.4$ $kN$ for the solid sample, respectively.

(a)

(b)

(c)

(d)

Figure 11: Example results for PVDF polymer rod: (a) - diameter change
registered by channel $\\#1$ for a precut sample, (b) - diameter change
registered by channel $\\#2$ for a precut sample. The initial force
$F_{max}\approx 4.6$ $kN$, (c) - diameter change registered by channel $\\#1$
for a solid sample, (d) - diameter change registered by channel $\\#2$ for a
solid sample. The initial force $F_{max}\approx 5.4$ $kN$

Complex results were obtained for measurements of PVC samples. For the precut
rod, neck initialisation was observed in both channels after approximately
$20s$ of elongation (see figs. 12(a) and 12(b)), while for the solid sample in
the second channel $\\#2$ periodic oscillations and diameter growth were
observed after $100s$ of elongation, opposite to the channel $\\#1$. The last
detected points for the precut sample after approximately $25s$ and $150s$ for
the solid sample correspond to the cracked rod with the shifted position on
the optical axis of the setup. The angle between the sample and the optical
axis may change with cracking and gives an asymmetric response in the
transverse channels of the setup. Furthermore, depending on the homogeneous
structure of the samples, the responses in the two channels could be
different. This proves the importance of direct bi-axial transverse
investigations using the optical setup. The behaviour of the material
corresponds to the complexity of the polymer structure and the response to
stress. Mechanical behaviour under tensile tests may not be repeatable for the
same polymer samples. The obtained results are crucial for further
investigations with specially selected, more homogeneous samples, e.g.
prepared with a particular form during manufacturing. The results show that
the samples must be precise and repeatable precuts to detect the start of the
neck to increase the measurements’ repeatability. Reducing the number of
measurement points can omit the sample vibration problem (see the appendix).

(a)

(b)

(c)

(d)

Figure 12: Example results for a PVC polymer rod: (a) - diameter change
registered by channel $\\#1$ for a precut sample, (b) - diameter change
registered by channel $\\#2$ for a precut sample. The initial force
$F_{max}\approx 3.4$ $kN$, (c) - diameter change registered by channel $\\#1$
for a solid sample, (d) - diameter change registered by channel $\\#2$ for a
solid sample. The initial force $F_{max}\approx 4.5$ $kN$

As can be seen, from the results and the calibration procedure, the proposed
setup provides $\mu m$ resolution for sample diameter detection during tensile
tests. The proposed low-cost setup may be considered the new measurement tool
for evaluating the direct mechanical properties of polymer materials. The data
obtained may be used to model the material’s properties numerically. The
repeatability and reproducibility of measurements with the tensile machine and
optical setup will begin presently.

##  4 Conclusions

The self-built laser projection system for symmetric transverse strain
measurements of polymers in tensile or compression tests has been presented.
The proposed device and measurement methodology allow for monitoring of
diameter changes in two perpendicular directions during specimen deformation.
The novel method, the constructed and calibrated setup, has been presented
with the first polymer tests using the tensile machine. The setup provided
monitor polymer diameters that can change from $12mm$ to even $1mm$ with
resolution $10$ $\mu m$. The maximum change in the diameter of the samples
presented was approximately $2mm$ before cracking a polymer rod. The response
of the sample to the load may be different in the orthogonal direction even at
the level of $50$ $\mu m$ (for example, see figs. 10(a) and 10(b)). For all of
the results, the polymer sample’s diameters decreased with elongation, which
confirmed the setup’s performance for monitoring on the tensile machine. The
noise of the data during measurements can be reduced with the number of
detected points (see the appendix). For data processing, the R programme was
prepared. The code is freely available on the Internet.
The optical setup was attached to the machine, replacing an extensometer and
preventing, in this way, the ability to simultaneously monitor the elongation
of the sample in the exact area of optical projector detection. Following this
issue, precision measurements will be conducted with well-prepared repeatable
samples, separately from the devices, and a statistical approach will be
proposed to analyse mechanical properties. The point is synchronisation
between machine and projector and repeatable conditions (homogeneous samples,
precise, repeatable geometry). In this way, the repeatability and
reproducibility of the measurements could also be evaluated.
The results have clearly shown that the proposed optical system may be a vital
tool for revising correct models of the mechanical properties of polymers
under complex loads. The low-cost setup provides fast and accurate sample
transversal dimension monitoring in tensile or compression tests, and further
work will focus on the robustness and compactness of the device. More accurate
and modern material models can be proposed based on registered data.
More comprehensive material research will begin with specially prepared
geometric and homogeneous material samples.

Notable progress in materials engineering requires more detailed research on
the mechanical properties of materials. Research on non-classical materials
requires adequate experimental results. The proposed methodology can be
applied for the comprehensive test of mechanical properties of non-classical
and linear or nonlinear behaviour materials, which can be more relevant in
material properties modelling. An additional advantage is the ability to test
samples of materials with a cylindrical or flat shape. The standard tensile or
compression test does not provide enough information on the mechanical
properties of non-classical materials such as a broad group of polymers.
Accurate registration of changes in vertical specimen dimensions is necessary
to describe modern material models. The appropriate material model is the
basis for all structural analyses. The proposed measurement system can be an
excellent tool for registering experimental data and mathematical descriptions
of the material’s mechanical properties under complex load conditions. The
proposed open-source data processing algorithms deliver the results’ clarity
and full reproducibility potential.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could appear to influence the work reported in
this document.

## Acknowledgments

The authors acknowledge Professor Szymon Wojciechowski of the Faculty of
Mechanical Engineering and Management, Poznan University of Technology, for
constructive criticism of the manuscript.
Funding: This work was supported by grant DS 061/SBAD/1547.

##  Appendix A The signal processing details

For signal processing, pracma and stats packages [26, 33] were used. Signal
processing steps are as follows:

  1. 1.

After the import of the data into R, the specially prepared moving average
function is applied for the initial data filtering. The function detects
possible NA’s when the signal has a discontinued character and prevents
filtering stops. Details can be found in the program [29]. The moving average
window with $80$ pixels was chosen for the presented result.

  2. 2.

The spline function smooth.spline() smoothed the filtered signal with the
desired equivalent number of degrees of freedom (trace of the smoother matrix)
$27$ [30, 31, 32]. In this way, the signal is prepared for an extrema
detection step.

  3. 3.

The predict function of the package car R [34] is used for the first
derivative calculation of the smoothed signal.

  4. 4.

Using the peaks function of the pracma R package, the extrema of the first
derivative are detected. The function returns the intensities as $y$s.

  5. 5.

To find the extrema positions in pixels ($x$’s) the which R function is used
for the signal edge detection (a shadow edge detection position).

  6. 6.

The estimated absolute difference between the positions of two points on the
edge of the shadow defined as the width ($P$ $[$Pixels$]$ parameter)
corresponds to the size of the measured object.

For preparing the manuscript graphs, the package tikzDevice was used [35].

##  Appendix B More result examples

The measurement results with a reduced number of captured points are presented
in figs. 13 and 14. The growths of the detected diameters at the end of the
tests correspond to the broken sample, which jumps out of the optical axis. In
this way, clear necking is observed before a break after $80$ seconds of
elongation for PET and after $150$ seconds for PVDF.

(a)

(b)

(c)

Figure 13: Example results for a solid PET rod with reduced detection points:
(a) diameter change registered by channel $\\#1$, (b) - diameter change
registered by channel $\\#2$, (c) - stress-strain curve obtained by the
tensile machine. The points correspond to registered data; the dashed line
corresponds to the b-spline fitting (degrees of freedom is $5$)

(a)

(b)

(c)

Figure 14: Example results for a solid PVDF rod with reduced points: (a) -
diameter change registered by channel $\\#1$, (b) - diameter change registered
by channel $\\#2$, (c) - stress-strain curve obtained by the tensile machine.
The points correspond to registered data; the dashed line corresponds to the
b-spline fitting (degrees of freedom is $5$)

##  Appendix C Linear fitting of the calibration graphs

For linear fitting of the models, the function $lm(formula,data)$ of R was
used. An object of class $formula$ is a symbolic description of the model to
be fitted. A $formula$ has an implicit intercept term of: $formula=y\sim x$. A
fixed zero-intercept fitting $formula=y\sim 0+x$ can be applied to remove
this. For the graph-based calibration method, the non-zero free intercept was
applied to give a better approximation to the linear function later used to
determine the diameter of samples in the $\phi=7-12$ mm range. The significant
error of approximation is clearly to be seen, especially for $\phi=10$ mm of
the zero intercept type fitting for both channels of the setup (see fig. 15)

(a)

(b)

(c)

(d)

Figure 15: Comparison of calibration curves (linear regression
$y=ax+b=\overline{P}=a\phi+b$) for two types of non-zero and zero intercepts
fitting conditions for channel $\\#1$: (a) – fitted parameters: $a=221.96$,
$b=141$, std. errors $\Delta a=2.47$, $\Delta b=15.36$, $p<0.001$; (b) –
fitted parameters: $a=241.11$, std. error $\Delta a=3.67$, $p<0.001$; and for
channel $\\#2$: (c) – fitted parameters: $a=228.36$, $b=92.77$, std. errors
$\Delta a=5.12$, $\Delta b=31.77$, $p<0.001$ and (d) – fitted parameters:
$a=241.42$, std. error $\Delta a=3.22$, $p<0.001$, respectively

##  Appendix D A note on b-spline fitting

The data’s global nonlinear fitting has several issues, the essential being
non-locality (proper function at a given value depends on data values far from
that point). In spline fitting, instead of fitting a global polynomial, the
range of $x$ is partitioned into smaller intervals, using an arbitrary number
and position of points, $K$, also called knots. Splines are created using a
polynomial with an interval between two successive knots. The spline has four
parameters in each $K+1$ region minus three constraints for each knot, giving
$K+4$ degrees of freedom. For example, a spline with three knots ($K=3$) has
seven degrees of freedom.

The R programming environment offers several packages that support b-spline
fitting. For the data processing presented in the paper, the package
smooth.spline() has been applied with specified degrees of freedom ($df$). The
lower value of the degrees of freedom meant stronger data smoothing. The
smoothing level is more dependent on the operator’s choice based on the number
of registered points rather than the error estimation procedure. For CCD
signals, it was optimised as $df=27$, influenced the estimated precision of
the system calibration $2\times\sigma<5$ pixels.

## References

  * [1] S. Tu, X. Ren, J. He, Z. Zhang, Stress–strain curves of metallic materials and post-necking strain hardening characterization: A review, Fatigue & Fracture of Engineering Materials & Structures 43 (1) (2020) 3–19.  doi:10.1111/ffe.13134.
URL https://onlinelibrary.wiley.com/doi/abs/10.1111/ffe.13134

  * [2] C. G’sell, J. J. Jonas, Determination of the plastic behaviour of solid polymers at constant true strain rate, Journal of Materials Science 14 (3) (1979) 583–591.  doi:10.1007/BF00772717.
URL http://link.springer.com/10.1007/BF00772717

  * [3] F. Addiego, A. Dahoun, C. G’Sell, J.-M. Hiver, Characterization of volume strain at large deformation under uniaxial tension in high-density polyethylene, Polymer 47 (12) (2006) 4387–4399.  doi:10.1016/j.polymer.2006.03.093.
URL https://linkinghub.elsevier.com/retrieve/pii/S003238610600406X

  * [4] C. G’Sell, J. M. Hiver, A. Dahoun, A. Souahi, Video-controlled tensile testing of polymers and metals beyond the necking point, Journal of Materials Science 27 (18) (1992) 5031–5039.  doi:10.1007/BF01105270.
URL http://link.springer.com/10.1007/BF01105270

  * [5] L. Nunes, F. Dias, H. da Costa Mattos, Mechanical behavior of polytetrafluoroethylene in tensile loading under different strain rates, Polymer Testing 30 (7) (2011) 791–796.  doi:10.1016/j.polymertesting.2011.07.004.
URL https://linkinghub.elsevier.com/retrieve/pii/S014294181100119X

  * [6] S. Todros, S. Pianigiani, N. de Cesare, P. G. Pavan, A. N. Natali, Marker Tracking for Local Strain Measurement in Mechanical Testing of Biomedical Materials, Journal of Medical and Biological Engineering 39 (5) (2019) 764–772.  doi:10.1007/s40846-018-0457-z.
URL http://link.springer.com/10.1007/s40846-018-0457-z

  * [7] W. Liu, X. Zhang, An optical method for measuring mixed-mode stress intensity factors for an arbitrarily oriented crack in cylindrical shells, Theoretical and Applied Fracture Mechanics 107 (February) (2020) 102567.  doi:10.1016/j.tafmec.2020.102567.
URL https://linkinghub.elsevier.com/retrieve/pii/S0167844220301427

  * [8] S. Muhammad, P.-Y. B. Jar, Effect of aspect ratio on large deformation and necking of polyethylene, Journal of Materials Science 46 (4) (2011) 1110–1123.  doi:10.1007/s10853-010-4880-y.
URL http://link.springer.com/10.1007/s10853-010-4880-y

  * [9] The International Bureau of Weights and Measures, The International System of Units (SI), Tech. Rep. 9 (sep 2019).
URL http://link.springer.com/10.1007/BF01449764

  * [10] R. Leach, Optical Measurement of Surface Topography, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.  doi:10.1007/978-3-642-12012-1.
URL http://link.springer.com/10.1007/978-3-642-12012-1

  * [11] D. Kucharski, H. Zdunek, A low-cost, simple optical setup for a fast scatterometry surface roughness measurements with nanometric precision, Bulletin of the Polish Academy of Sciences: Technical Sciences 68 (48) (2020) 485–490.  doi:10.24425/bpasts.2020.133364.
URL
http://journals.pan.pl/Content/116520/PDF/10_485-490_01475_Bpast.No.68-3_30.06.20_K3A_TeX.pdf

  * [12] L. Zhang, X. Song, G. Cao, Z. Yu, Optical measurement method of necking process in the material tensile test, in: S. Ye (Ed.), Automated Optical Inspection for Industry: Theory, Technology, and Applications II, Vol. 3558, International Society for Optics and Photonics, SPIE, 1998, pp. 615–617.  doi:10.1117/12.318451.
URL http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=971744

  * [13] R. Srinivasan, C. S. Hartley, B. B. Raju, J. Clave, Measurement Of Neck Development In Tensile Testing Using Projection Moire, Optical Engineering 21 (4) (1982) 655–662.  doi:10.1117/12.7972961.
URL
http://opticalengineering.spiedigitallibrary.org/article.aspx?doi=10.1117/12.7972961

  * [14] H. Hoffmann, C. Vogl, Determination of True Stress-Strain-Curves and Normal Anisotropy in Tensile Tests with Optical Strain Measurement, CIRP Annals 52 (1) (2003) 217–220.  doi:10.1016/S0007-8506(07)60569-1.
URL https://linkinghub.elsevier.com/retrieve/pii/S0007850607605691

  * [15] J. Ye, S. André, L. Farge, Kinematic study of necking in a semi-crystalline polymer through 3D Digital Image Correlation, International Journal of Solids and Structures 59 (2015) 58–72.  doi:10.1016/j.ijsolstr.2015.01.009.
URL https://linkinghub.elsevier.com/retrieve/pii/S0020768315000116

  * [16] J. Fischer, T. Radil, Simple Methods of Edge Position Measurement Using Shadow Projected on CCD Sensor, Meas. Sci. Rev. 3 (2003) 1–4.
URL https://www.measurement.sk/2003/S3/Radil.pdf

  * [17] S. Kato, D. Aoki, K. Oikawa, K. Tsuchiya, K. Numata, H. Otsuka, Visualization of the Necking Initiation and Propagation Processes during Uniaxial Tensile Deformation of Crystalline Polymer Films via the Generation of Fluorescent Radicals, ACS Macro Letters 10 (5) (2021) 623–627.  doi:10.1021/acsmacrolett.1c00185.
URL https://pubs.acs.org/doi/10.1021/acsmacrolett.1c00185

  * [18] W. Ren, G. Yang, B. K. Mukherjee, J. P. Szabo, Interferometric measurement of the transverse strain response of electroactive polymers, Smart Structures and Materials 2004: Electroactive Polymer Actuators and Devices (EAPAD) 5385 (2004) 395.  doi:10.1117/12.538531.
  * [19] A. G. Leal-Junior, A. Theodosiou, C. Marques, M. J. Pontes, K. Kalli, A. Frizera, Compensation Method for Temperature Cross-Sensitivity in Transverse Force Applications With FBG Sensors in POFs, J. Lightwave Technol. 36 (17) (2018) 3660–3665.
URL http://opg.optica.org/jlt/abstract.cfm?URI=jlt-36-17-3660

  * [20] A. Leal-Junior, A. Frizera-Neto, C. Marques, M. Pontes, A Polymer Optical Fiber Temperature Sensor Based on Material Features, Sensors 18 (2) (2018) 301.  doi:10.3390/s18010301.
URL http://www.mdpi.com/1424-8220/18/1/301

  * [21] A. Leal-Junior, A. Frizera, M. J. Pontes, A. Fasano, G. Woyessa, O. Bang, C. A. F. Marques, Dynamic mechanical characterization with respect to temperature, humidity, frequency and strain in mPOFs made of different materials, Optical Materials Express 8 (4) (2018) 804.  doi:10.1364/OME.8.000804.
URL https://opg.optica.org/abstract.cfm?URI=ome-8-4-804

  * [22] A. Leal-Junior, A. Frizera, H. Lee, Y. Mizuno, K. Nakamura, T. Paixão, C. Leitão, M. F. Domingues, N. Alberto, P. Antunes, P. André, C. Marques, M. J. Pontes, Strain, temperature, moisture, and transverse force sensing using fused polymer optical fibers, Optics Express 26 (10) (2018) 12939.  doi:10.1364/OE.26.012939.
URL https://opg.optica.org/abstract.cfm?URI=oe-26-10-12939

  * [23] A. G. Leal-Junior, C. Marques, M. R. N. Ribeiro, M. J. Pontes, A. Frizera, FBG-Embedded 3-D Printed ABS Sensing Pads: The Impact of Infill Density on Sensitivity and Dynamic Range in Force Sensors, IEEE Sensors Journal 18 (20) (2018) 8381–8388.  doi:10.1109/JSEN.2018.2866689.
URL https://ieeexplore.ieee.org/document/8444379/

  * [24] R Foundation, The R Project for Statistical Computing (2002).
URL https://www.r-project.org

  * [25] D. Kurpisz, M. Obst, The energetic and experimental based approach to description of basic material characteristics and mechanical properties of selected polymers, Journal of Theoretical and Applied Mechanics 58 (1) (2020) 183–193.  doi:10.15632/jtam-pl/115392.
URL http://www.journalssystem.com/jtam/,115392,0,2.html

  * [26] H. W. Borchers, pracma: Practical Numerical Math Functions, r package version 2.3.3 (2021).
URL https://CRAN.R-project.org/package=pracma

  * [27] E. Pink, A. Kronthaler, P. Valtingoier, Double necking in polymer tensile deformation, Journal of Materials Science 24 (1) (1989) 183–186.  doi:10.1007/BF00660951.
URL http://link.springer.com/10.1007/BF00660951

  * [28] E. R. Torbensen, Esben Rossel Torbensen website about the TCD1304 sensors, https://tcd1304.wordpress.com, Accessed: 2021-09-11.
  * [29] D. Kucharski, The R code for CCD signal processing and data evaluation used in a self-built laser projection system for polymers neck development measurements in tensile testing (2022).  doi:10.5281/zenodo.5814853.
URL https://github.com/dawidkucharski/optical_neckle_measurements

  * [30] C. de Boor, On calculating with B-splines, Journal of Approximation Theory 6 (1) (1972) 50–62.  doi:10.1016/0021-9045(72)90080-9.
URL https://linkinghub.elsevier.com/retrieve/pii/0021904572900809

  * [31] A. Perperoglou, W. Sauerbrei, M. Abrahamowicz, M. Schmid, A review of spline function procedures in R, BMC Medical Research Methodology 19 (1) (2019) 46.  doi:10.1186/s12874-019-0666-3.
URL
https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-019-0666-3

  * [32] E. Quak, About b-splines. twenty answers to one question: What is the cubic b-spline for the knots -2,-1,0,1,2?, J. Numer. Anal. Approx. Theory 45 (1) (2016) 37–83.
URL https://ictp.acad.ro/jnaat/journal/article/view/1099

  * [33] R Core Team, R: A Language and Environment for Statistical Computing (2020).
URL https://www.r-project.org/

  * [34] J. Fox, S. Weisberg, An R Companion to Applied Regression, 3rd Edition, Sage, Thousand Oaks CA, 2019\.
URL https://socialsciences.mcmaster.ca/jfox/Books/Companion/

  * [35] C. Sharpsteen, C. Bracken, tikzDevice: R Graphics Output in LaTeX Format, r package version 0.12.3.1 (2020).
URL https://CRAN.R-project.org/package=tikzDevice