# Micro-scale Mechanical Metamaterial with a Controllable Transition in the
Poisson’s Ratio and Band Gap Formation

K. K. Dudek  Institut FEMTO-ST, CNRS, Université Bourgogne Franche-Comté,
Besançon 25030, France  Institute of Physics, University of Zielona Gora, ul.
Szafrana 4a, Zielona Gora 65-069, Poland  J. A. Iglesias Martínez  Institut
FEMTO-ST, CNRS, Université Bourgogne Franche-Comté, Besançon 25030, France  G.
Ulliac  Institut FEMTO-ST, CNRS, Université Bourgogne Franche-Comté, Besançon
25030, France  L. Hirsinger  Institut FEMTO-ST, CNRS, Université Bourgogne
Franche-Comté, Besançon 25030, France  L. Wang  Institut FEMTO-ST, CNRS,
Université Bourgogne Franche-Comté, Besançon 25030, France  V. Laude  Institut
FEMTO-ST, CNRS, Université Bourgogne Franche-Comté, Besançon 25030, France  M.
Kadic  Institut FEMTO-ST, CNRS, Université Bourgogne Franche-Comté, Besançon
25030, France

###### Abstract

The ability to change significantly mechanical and wave propagation properties
of a structure without rebuilding it has been one of the main challenges in
the field of mechanical metamaterials. This stems from the enormous appeal
that, especially in the case of micro-scale systems, such tunable behavior may
offer from the perspective of applications ranging from biomedical to
protective devices. In this work, a novel micro-scale mechanical metamaterial
is proposed that can undergo a transition from one type of configuration to
another, with one configuration having a very negative Poisson’s ratio,
corresponding to strong auxeticity, and the other having a highly positive
Poisson’s ratio. The formation of phononic band gaps, at the same time, can be
controlled, which can be very useful in the design of vibration dampers and
sensors. Finally, it is shown experimentally that reconfiguration of the
system, leading to a change in its properties, can be induced and controlled
remotely through application of a magnetic field, thanks to appropriately
distributed magnetic inclusions.

##  1 Introduction

Over the years, the importance of functional materials [1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11] has been rapidly increasing in our everyday lives. Most
prominently in the span of the last two decades, different classes of these
structures have been successfully implemented in numerous industries including
biomedical [12, 13] and protective devices [14, 15, 16], sports equipment [17]
as well as electronics [18, 19]. One of the most promising directions of
studies related to such easily-applicable functional materials is that devoted
to mechanical metamaterials [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32, 33, 34, 35, 36, 37]. Mechanical metamaterials are rationally designed
structures that can exhibit counterintuitive mechanical properties based
primarily on their design rather than their chemical composition. Some of the
most commonly studied instances of such extraordinary mechanical properties
are negative Poisson’s ratio (auxetic behavior) [20, 38, 39, 40, 41], negative
stiffness [42, 43, 44, 45] and negative compressibility [42, 46, 47, 48]. Most
notably, auxetic mechanical metamaterials have been proven to exhibit superior
indentation resistance [49], energy absorption [50, 51] and wave attenuation
[14, 52] in comparison to many conventional non-auxetic materials. However,
despite the numerous advantages offered by mechanical metamaterials
characterized by a negative value of Poisson’s ratio, a vast majority of such
structures share an important limitation. Namely, once typical auxetic
metamaterials are manufactured, it is very difficult to change their
mechanical characteristic without rebuilding the system. This, in turn,
diminishes their applicability in the fields where variable response of the
material is required. Nevertheless, as reported in recent years, this
shortcoming can be addressed by the use of active mechanical metamaterials in
a form of stimuli responsive structures.

Active mechanical metamaterials are structures that can exhibit a tunable
mechanical behavior depending on the application of an external stimulus. In
recent years, it has been demonstrated that in the case of multi-material
structures, the mechanical properties of the system can be adjusted remotely
in a process that does not require reconstructing the system. In the case of
macroscopic systems, one of the most promising approaches corresponds to the
use of appropriately distributed magnet-like inclusions that respond to the
application of an external magnetic field [53, 54]. However, this task becomes
much more challenging at small scales where one of the more feasible
approaches seems to be the use of magnetic nanoparticles dispersed within a
non-magnetic host material [10, 55, 56, 57]. Another approach making it
possible to construct an active mechanical metamaterial is associated with the
use of several materials with different thermal expansion coefficients. In
such a case, the behavior of the structure can be modified by changing the
external temperature [56, 58]. This method, however, often does not allow to
obtain a very significant variation in the mechanical properties. In addition,
the corresponding reconfiguration process tends to be quite slow. Last but not
least, it is worth noting that it is possible to construct an active
metamaterial composed of only one type of residual material that, to some
extent, can undergo a change in its mechanical properties with a deformation
pattern based solely on a variation in the mechanical stimulus, e.g. the rate
of the mechanical compression. A prime example of this approach corresponds to
hierarchical metamaterials [59, 60, 61, 62, 63, 64, 65], where different
hierarchical levels can normally deform irrespective of each other.

Despite several notable achievements, studies devoted to active mechanical
metamaterials are still in their infancy and much more in-depth analysis must
be conducted to fully harness their potential. First of all, it is important
to emphasize that for a significant number of known active mechanical
metamaterials, even though it is possible to change one of their properties
without rebuilding the system, it is still very difficult to achieve control
over the other properties. This, in turn, could increase the appeal for active
functional materials from the perspective of many industries. In fact, even
though this task is not simple, it becomes even more challenging at small
scales, such as the micro-scale, where the manufacturing process is often not
only very difficult but also fairly expensive. On the other hand, despite the
aforementioned challenges, the ability to achieve tunable control over
multiple properties of a micro-scale system without the need for
reconstructing it appears to be a very appealing concept that could be
utilized in fields such as programmable robotics, vibration dampers / sensors
and effective biomedical devices. To this aim, a very promising idea seems to
be the possibility of controlling the Poisson’s ratio of a micro-scale system
[64] as well as its wave propagation properties with an emphasis on tunable
band gap formation [52, 66, 67]. In fact, in the literature [68], there are
some rare examples demonstrating the possibility of controlling to a limited
extent both Poisson’s ratio and band gap formation of micro-scale mechanical
metamaterials. However, despite their numerous advantages, such studies
typically focus on complex hierarchical structures where the deformation
process is not easily recoverable nor controlled due to the presence of
multiple independent hierarchical levels. In many cases, it is also not
possible to observe a large tunable change in the Poisson’s ratio of the
system, which would significantly increase the applicability of the considered
concept.

In this work, a novel mechanical metamaterial is proposed that in general can
be defined in the form of two distinct configurations. It is presented that
without rebuilding the system, the considered structure can undergo a
transition from one type of configuration to another during which the
Poisson’s ratio varies very significantly. Namely, the change occurs from a
structure possessing a positive Poisson’s ratio to a strongly auxetic
structure and vice versa. In addition, in terms of wave propagation
properties, it is shown that a change in the configuration of the system can
be used to achieve control over the formation of a phononic band gap.
Furthermore, it is presented that all configurations of the considered system
can be constructed at the micro-scale where they exhibit the aforementioned
tunable properties. Finally, it is shown experimentally that reconfiguration
of the system, leading to a change in its properties, can be prompted and
controlled remotely through the use of an external magnetic field.

##  2 Results and Discussion

In this work, a novel mechanical metamaterial composed of triangle-shaped
structural elements is proposed in order to assess its potential to exhibit
versatile mechanical and wave propagation properties. As shown in Fig. 1a, the
considered system can occur in two very different configurations, that in this
work, are named Type A and Type B structures. In the case of the Type A
configuration, the system consists of two types of isosceles triangles
connected to each other at vertices. Due to a specific connectivity, this
structure corresponds to three different types of apertures. In contrast, the
configuration defined as Type B has a much more regular design and consists of
only one type of isosceles triangles. In fact, it resembles a lattice of
mutually connected star-shaped elements where each of these elements consists
of six triangles. The two considered configurations were designed in such a
way that during mechanical deformation they would exhibit a very different
Poisson’s ratio. This concept is depicted schematically in Fig. 1a, where
qualitatively, type A and Type B structures exhibit strongly negative (auxetic
behavior) and positive Poisson’s ratio, respectively.

Figure 1: a) Diagrams showing the two types of structures considered in this
work. b) Experimental prototypes corresponding to specific examples of the
Type A and Type B structures that are analyzed from the point of view of their
mechanical properties. Parameters that remain the same for all structures are
the following: $a$ = 40 $\mu$m, $l_{1}$ = 4.8 $\mu$m, $l_{2}$ = 2.4 $\mu$m,
$l_{3}$ = 2.4 $\mu$m, $l_{4}$ = 3.6 $\mu$m, and $l_{5}$ = 6.7 $\mu$m. In
addition, parameters that remain the same for type A structures are $b$ = 80
$\mu$m and $c$ = 40 $\mu$m. Parameters that assume unique values both for type
A and type B systems are as follows: $d$ = 28 $\mu$m (A1), $d$ = 40 $\mu$m
(A2), $d$ = 60 $\mu$m (A3) and $b$ = 40 $\mu$m (B1), $b$ = 60 $\mu$m (B2), and
$b$ = 80 $\mu$m (B3). The out-of-plane thickness $l$ for all of the structures
is equal to 300 $\mu$m.

In order to investigate the mechanical properties of the two considered
configurations and to demonstrate that the proposed concept can be utilized in
the case of microscopic applications, experimental prototypes were constructed
at the micro-scale (see Methods section). In addition, as shown in Fig. 1b, in
order to assess the effect of the geometric parameters on the behavior of the
structure, for each configuration three different prototypes were fabricated
that are labeled A1, A2, A3 and B1, B2, B3. In fact, the only difference
between the three Type A structures is the value of the $d$ parameter.
Similarly, the only difference between the Type B prototypes is the value of
the $b$ parameter. All prototypes are composed of 5 $\times$ 3 unit-cells,
where the bottom part of the system is clamped as a result of the
manufacturing process (see Methods section).

###  2.1 Mechanical Properties

To determine the Poisson’s ratio of the considered system, the microscopic
experimental prototypes corresponding both to Type A and Type B were
compressed by a flat external indenter along the $y$-axis, with the bottom
part of each sample remaining fixed in space. Thus, to qualitatively estimate
the Poisson’s ratio of the system, one can follow the changes in the
horizontal dimension of the topmost part of the structure that is the farthest
part of each sample from the aforementioned constraint. According to optical
images of the experiment presented in Fig. 2a, one can note that from the
qualitative analysis, prototypes corresponding to Type A and Type B exhibit a
very different behavior. Namely, the $x$ dimension for all of the prototypes
associated with the Type A configuration decreases significantly during the
compression process. This, in turn, is an indication of a strong auxetic
behavior. In contrast, in the case of structures corresponding to Type B, the
horizontal dimension of the samples initially increases before returning
approximately to the initial value at large strains. Such behavior is an
indication of positive and close to zero Poisson’s ratio, respectively. This
means that Type B configurations exhibit a non-auxetic behavior unlike the
Type A samples. At this point, it should be emphasized that non-uniform
deformation of the entire system originates from the fact that the
experimental samples remain attached to the substrate after the 3D printing
process. However, if one was to consider the deformation of the detached
sample, the extent of the non-uniformity would be significantly smaller.
Nonetheless, mechanical testing could potentially become unreliable in such a
scenario due to the very small size of the system. It is also worth
emphasizing that if fabricated through the use of appropriate materials, the
considered system can maintain its properties over numerous deformation cycles
(see Fig. S1 in the Supplementary Information). This, in turn, further
enhances its appeal from the perspective of potential applications.

‘

Figure 2: a) Optical images of a compression test along the $y$-axis of the
three micro-scale experimental prototypes corresponding to Type A and Type B
configurations. b) Poisson’s ratio estimated from experimental samples. c)
Poisson’s ratio estimated with FEM (Finite Element Method) computer
simulations that was calculated based on the deformation of the central unit-
cell of the considered finite samples.

In addition to the qualitative description of the behavior of the system, it
is important to conduct a quantitative analysis in order to determine the
range of Poisson’s ratio values and to assess its usefulness in the case of
specific applications. To this aim, in order to calculate the Poisson’s ratio
for each of the structures, the central unit-cell was taken into
consideration. This stems from the fact that the mechanical behavior of this
specific unit-cell is the least affected by the factors contributing to the
uneven deformation pattern. Namely, the deformation of this unit-cell is not
significantly affected by the bottom part of the system being attached to the
substrate. In addition, its behavior is also mostly unaffected by the edge
effects on the left and right-hand sides of the structure as well as by the
frition related to the contact of the topmost part of the system with the
indenter inducing the deformation in the conducted experiments. As shown in
Fig. 2b, all of the Type A structures exhibit a strong auxetic behavior
throughout the entire deformation process. In fact, the largest extent of
auxeticity can be observed for the A3 sample corresponding to the largest
value of $d$. In this case, for a relatively large strain of around 0.07, the
Poisson’s ratio assumes the value of approximately -1.4. Furthermore, one can
note that the extent of auxeticity can be adjusted by the value of $d$. In
contrast, the behavior of the analyzed samples of type B is very different.
Namely, all type B samples tend to exhibit a positive Poisson’s ratio that in
most cases, becomes gradually less positive as strain increases. Notably, the
sample that exhibits the largest values of Poisson’s ratio throughout the
majority of the deformation process is structure B3 that possesses the largest
value of $b$ amongst the considered samples. Thus, taking all of these results
into account, there is a significant difference between the Poisson’s ratios
of type A and type B structures. In addition, the range of Poisson’s ratio for
each configuration can be adjusted by a variation in geometric parameters.

In addition to experimental studies, another interesting aspect of the
analysis corresponds to the investigation of the Poisson’s ratio of the system
through the FEM (Finite Element Method) computer simulations. In this case,
both the structure and the corresponding boundary conditions were designed to
closely match the experiment (see Methods section). As shown in Fig. 2c, the
trends observed for computational results were qualitatively the same as for
conducted experiments. Namely, in the case of these two approaches, type A
configurations exhibited a strong auxetic behavior throughout the entire
deformation process while type B configurations exhibited positive values of
this parameter. Of course, in terms of the quantitative analysis, there are
some discrepancies between the two approaches that originate from multiple
reasons. Namely, in the case of the experiment, the effect of the bottom
constraint and edge effects on the effective Poisson’s ratio is not fully
negligible. However, through the use of multiple cells in the design of the
system and the selection of only the central unit-cell for analysis purposes,
these effects are certainly significantly diminished. Thus, based on the
provided results, it is clear that the two approaches confirm a very
significant difference in the Poisson’s ratio exhibited by the two types of
analyzed configurations.

###  2.2 Effect of Reconfiguration on the Properties of the Structures

In addition to the analysis of Poisson’s ratio, another interesting direction
of study is the analysis of phononic band gap formation associated with phonon
dispersion. This stems from the fact that such analysis provides valuable
information about the ranges of frequencies for which elastic waves are not
transmitted through a given configuration of the system. Hence, in addition to
applications related explicitly to a tunable Poisson’s ratio, such information
can be used in the design of microscopic vibration dampers / sensors.

From the perspective of applications, it is very important to check what
happens to mechanical and wave propagation properties once the considered
system is subjected to a reconfiguration process. In other words, it is
essential to assess the properties of the system during a reconfiguration
process that does not involve reconstructing the structure and that could
potentially occur as a result of the application of an external stimulus. An
example of such a reconfiguration process is presented in Fig. 3a, where the
initial configuration of type A corresponds to angle $\theta_{0}=80^{\circ}$.
In this case, all triangles are the same and correspond to the $b/a$ = 1.5
aspect ratio. Furthermore, as the reconfiguration process begins, the
magnitude of the angle $\theta_{0}$ starts increasing until it reaches the
value $\theta_{0}=120^{\circ}$. At this point, the system assumes a
configuration of type B, where there are only two types of apertures. To
describe what happens to the properties of the system throughout such a
reconfiguration process, we focus on the three representative configurations
depicted schematically in Fig. 3a, where $\theta_{0}$ is equal to
$80^{\circ}$, $100^{\circ}$ and $120^{\circ}$ respectively.

‘

Figure 3: Effect of the reconfiguration of the system on phononic band gap
formation and mechanical properties. a) Diagrams showing schematically
different configurations assumed by the same system. b) Phonon dispersion
relations for the three considered configurations where the first band gap is
highlighted by means of the grey box. c) Variation of Poisson’s ratio for the
three selected structures subjected to mechanical compression along the
$y$-axis. Poisson’s ratio was calculated for individual unit-cells with
periodic boundary conditions imposed in the $x$-direction and a symmetric
prescribed displacement applied to the topmost and bottommost parts of the
system along the $y$-axis.

According to Fig. 3b, the phononic band structures for all three considered
configurations are relatively similar. Most notably, the first non-negligible
band gap occurring in the case of the system characterized by
$\theta_{0}=80^{\circ}$ has a size of 2.95 MHz. As the value of $\theta_{0}$
increases and assumes the values of $100^{\circ}$ and $120^{\circ}$, the size
of the band gap increases and is equal to 4.03 MHz and 4.13 MHz respectively.
However, even though the size of the band gap changes as a result of the
reconfiguration process, one could say that such variations are not very
significant. This, in turn, agrees with the observations made for other
complex mechanical metamaterials subjected to a reconfiguration process that
were reported in the literature [68]. At this point, it is also important to
note that whereas for Fig. 3 a structure corresponding to a specific value of
the b/a ratio was taken into account, in general, the change in width of the
first significant band gap may look different when considering the deformation
of structures characterized by different values of $b/a$. More specifically,
as described in the Supplementary Information, if one was to consider the
deformation of the structure associated with $b/a$ = 2 from the configuration
where $\theta_{0}$ = 80∘ to the configuration corresponding to $\theta_{0}$ =
120∘, the magnitude of the change in the band gap width would be very
different from the former case. Namely, in this scenario, the relative change
in the band gap width is only equal to 0.149 MHz. In contrast, for the
structure associated with $b/a$ = 1.5, the change in the band gap width equals
1.18 MHz. Thus, a relatively small variation in the aspect ratio of structural
elements may lead to a considerable difference in the magnitude of the change
in the width of the band gap during the continuous reconfiguration process
associated with the change in $\theta_{0}$. Of course, this stems from the
fact that structures corresponding to certain values of $b/a$ lead to
considerably different band gaps compared to other systems. Hence, it follows
that the change in their magnitude is also considerably different.

Contrary to the results related to band gap formation, it appears that the
reconfiguration process strongly influences the Poisson’s ratio of the system.
As shown in Fig. 3c, configurations associated with $\theta_{0}=80^{\circ}$
and $\theta_{0}=120^{\circ}$ assume very different values of Poisson’s ratio.
Namely, in the case of the system where $\theta_{0}=80^{\circ}$, the structure
exhibits a strong auxetic behavior throughout the entire deformation process.
On the other hand, for $\theta_{0}=120^{\circ}$, the value of the Poisson’s
ratio is predominantly positive. These two results are in accordance with the
expectations that one could have based on the results related to the
mechanical testing of different configurations of Type A and Type B. However,
it is very interesting to see what happens when the Type A configuration
becomes very similar to a Type B configuration. In fact, such a scenario can
be observed in the case of the intermediate configuration corresponding to
$\theta_{0}=100^{\circ}$. For this configuration, even though Poisson’s ratio
is initially positive, it quickly assumes strong negative values upon
increasing strain. All of these results are very interesting since they
indicate that a relatively simple reconfiguration process of the considered
system makes it possible to very significantly change its Poisson’s ratio.

Figure 4: Effect of the elongation of the isosceles triangles corresponding to
Type B of the system on the band gap formation. a) Diagrams showing three
different configurations of the system, where the value of $a$ is kept
constant. b) Phonon band structures corresponding to the analyzed structures.
c) The first three band gaps associated with structures corresponding to
different values of the $b/a$ ratio. The geometric parameters used for all of
the structures were set to be the following: $a$ = 200 $\mu$m, $l_{1}$ = 8
$\mu$m, $l_{2}$ = 4 $\mu$m, $l_{3}$ = 4 $\mu$m, $l_{4}$ = 4 $\mu$m, $l_{5}$ =
9 $\mu$m

Once the effect of the reconfiguration of the system on its properties has
been determined, it is also interesting to have a closer look at band gap
formation for the Type B system. As shown in Fig. 3b, the largest first
significant band gap is observed during the reconfiguration process. To this
aim, it was analyzed what happens when the aspect ratio of isosceles triangles
constituting the Type B structure changes while keeping dimension $a$
constant.

The three specific configurations of the system chosen for such analysis are
portrayed schematically in Fig.4a. In fact, the only difference between them
is the value of $b$. Hence, to refer to a specific configuration assumed by
the structure, one can use a particular value of the $b/a$ coefficient.
According to Fig.4b, the variation in the value of $b/a$ has a very
significant effect on the phonon dispersion and on the size of the first
significant band gap. More specifically, in the considered range of
frequencies, one can note that for the structure associated with $b/a$ = 0.7,
the first non-negligible band gap is relatively small. In contrast, once the
geometric dimensions of the structure are modified to correspond to $b/a$ =
1.0, the appearance of a very significant band gap can be observed.
Furthermore, upon increasing the magnitude of this coefficient up to the value
of 1.3, the size of the first significant band gap is reduced relative to the
configuration associated with $b/a$ = 1.0. This means that at relatively low
frequencies, one should expect the largest band gap to be observed for this
specific configuration. Hence, as shown in Fig. 4c, an analysis of the size of
the first three band gaps was conducted for a larger number of possible
structures corresponding to the $b/a$ coefficient ranging between 0.7 and 1.3.
Based on the results, the width of the first significant band gap gradually
increases until it assumes its maximum value at $b/a$ = 1.0. For
configurations with $b/a$ exceeding this value, the band gap gradually
decreases. Finally, it should be noted that the experimental validation of the
possibility of controlling the width of band gaps through the change in the
value of $b/a$ is provided in Fig. S2 in the Supplementary Information.

###  2.3 Active Control

Up to this point, it has been shown that the properties of the proposed
mechanical metamaterial can be controlled as a result of the reconfiguration
process. However, such reconfiguration was merely a theoretical concept and it
was not explicitly demonstrated how it could be achieved without
reconstructing the system. Thus, in this section, it is shown how such active
control over the configuration assumed by the system and the resulting
mechanical properties can be attained. Furthermore, since the concept proposed
in this work is fully scalable and all results can in general be reproduced
irrespective of the size of the system, active control over the
reconfiguration process was demonstrated through the use of a macroscopic
prototype (see Fig. 5).

Figure 5: Active reconfiguration and transition in the properties of the
system induced by the application of an external magnetic field. a) Picture of
the investigated sample placed between two parallel electromagnets that
produce an approximately uniform magnetic field. b) Diagram showing the
concept of the reconfiguration of the magneto-mechanical structure. c) A
change in the configuration of the system induced by the variation in the
magnitude of the external magnetic field. d) Mechanical deformation of
different configurations assumed by the system during the active
reconfiguration.

As shown in Fig. 5a, the system selected to demonstrate the active
reconfiguration process corresponds to a 3D-printed (see Methods section)
prototype associated with the Type B configuration that was inserted into the
uniform magnetic field induced by two large parallel electromagnets. The
response of the structure to the application of the external magnetic field
was made possible due to the use of magnetic inclusions in the form of very
small neodymium magnets (see Methods section) that were embedded within some
of the structural triangle-like elements constituting the system. The concept
of active reconfiguration corresponding to this structure is schematically
demonstrated in Fig. 5b, where for the sake of clarity, only a fragment of the
considered structure is presented. As can be seen, the structure is initially
at rest in the Type B configuration, where $\theta_{0}$ = 120∘. Once the
external magnetic field is switched on, the magnetic dipole moments of the
magnets attempt to align with the external magnetic field. Of course,
depending on the magnitude of the external field, this might not be possible
since for the triangles to rotate, the applied torques must be large enough to
overcome the resistance and the restoring forces corresponding to the elastic
non-magnetic material that the structure is made of.

According to Fig. 5c, upon being subjected to a relatively strong external
magnetic field having a magnitude of 0.24 T, the considered structure deforms
significantly and assumes a configuration corresponding to $\theta_{0}\approx
84^{\circ}$ (see Supplementary Video 1 to observe the gradual reconfiguration
process where the magnetic field was changed with increments of 0.01 T). In
other words, the structure undergoes a transition from Type B to Type A
configuration. Based on the computational results provided in Fig. 5d, one can
note that these two configurations correspond to a very different Poisson’s
ratio. Thus, the application of the external magnetic field makes it possible
to significantly modify the properties of the structure without the need to
rebuild it. In addition, after switching off the external magnetic field, the
structure returns to its initial configuration, which means that after the
transition it is possible to recover the initial properties of the system.
These results are very important as they prove that the mechanical
metamaterial proposed in this work can be used to exhibit a versatile
Poisson’s ratio that can be controlled remotely via the magnitude and
orientation of the external magnetic field. Of course, in this case, due to
the scale of the prototype, the controllable behavior was achieved thanks to
the use of magnetic inclusions in the form of neodymium magnets. However, at
much lower scales, this approach should be replaced for example by the use of
single domain superparamagnetic nanoparticles dispersed within the resin prior
to the polymerization process [69, 70]. Another approach that could be used to
achieve control over the behavior of the system through the use of an external
magnetic field, is to use a selective magnetic coating applied for example by
sputtering following the 3D printing fabrication process [71, 72]. At this
point, we emphasize that the results of this work are scale independent.
Depending on the size of the system, however, they can be better suited for
specific applications.

The results show that the proposed system can exhibit two very useful
properties at the same time that are further controllable remotely via an
applied external magnetic field without the need for reconstructing the
system. Specifically, it was shown that the considered structure can undergo a
transition from positive to highly negative Poisson’s ratio solely as a result
of a reconfiguration process that can be invoked and controlled remotely and
that does not involve rebuilding the system. In addition, after the
transition, the system maintains a stable mechanical behaviour for relatively
large strains. Such active and reversible control over mechanical properties,
combined with the observation of a very significant variation of Poisson’s
ratio, is seldom observed in mechanical metamaterials and can be very useful
for numerous applications. Furthermore, what makes these results even more
appealing, is the fact that it is demonstrated that the highly versatile
Poisson’s ratio of the considered system can be observed for samples
fabricated at the micro-scale. In addition, from the perspective of
multifunctional devices, the proposed system exhibits two highly useful
properties at once: controllable Poisson’s ratio and adjustable phononic band
gap. The latter property prevents waves in a specific range of frequencies
from propagating through the system. Mechanical metamaterials exhibiting such
band gap formation combined with fully controllable Poisson’s ratio are
extremely rare. In a few known instances where control over both aspects was
attempted to some extent [56], the studies were limited to structures at the
macro / millimetre scale and typically only one of these properties was
quantitatively assessed.

As mentioned above, the reconfiguration process leading to the transition in
the properties of the system is fully reversible and can be achieved in an
active manner that does not require rebuilding the system. This means that the
ability of the mechanical metamaterial to exhibit a tunable Poisson’s ratio
could be used at the macroscale or the millimeter scale in the design of smart
protective devices that could adjust their mechanical response depending on
the specific impacting body, e.g. active car bumpers. Similarly, at the
microscale, the controllable Poisson’s ratio can lead to the design of
materials characterized by programmable extent of their micro-indentation
resistance. In addition, should somebody use magnetic inclusions to control
the properties of the structure, then by distributing them in a non-uniform
manner, it would be possible to observe a situation where the initially
uniformly-designed structure exhibits a very different Poisson’s ratio in
different parts of the system. Such behavior, as demonstrated in another study
[64], could enable the system to exhibit complex shape morphing upon being
subjected to a mechanical deformation or through the application of an
external magnetic field. Especially at the micro-scale, and among other
applications, such shape morphing may be useful for the design of smart stents
that would provide local support to specific parts of a blood vessel, where
the location of the support could be changed throughout treatment. Finally, it
was also demonstrated that in addition to control over the Poisson’s ratio,
one can form band gaps that are almost unaffected by the reconfiguration
process. Thus, while changing the mechanical response of the sample, it would
be possible to prevent waves corresponding to a specific range of frequencies
from being transmitted through the system. Such ability could be very
interesting from the perspective of the design of multifunctional vibration-
damping protective devices. Alternatively, should the need arise, the band gap
could be adjusted by a variation in geometric parameters of the system. In
addition, for specific values of the $b/a$ parameter, the change in the width
of the band gap induced by mechanical deformation, which process does not
involve the reconstruction of the system, can be noticeable.

##  3 Conclusion

In this work, a novel micro-scale mechanical metamaterial was proposed. It was
demonstrated that depending on the assumed configuration, the system can
exhibit very different Poisson’s ratios including strongly negative and
positive values at relatively large strains. In addition, it was shown that
such control over the mechanical behavior of the system can be achieved as a
result of active reconfiguration, i.e. a process that does not require
rebuilding the system. For such reconfiguration process, the transition in the
properties of the structure is recoverable and can be controlled for example
by means of an external magnetic field. Furthermore, in addition to control
over Poisson’s ratio, it was demonstrated that during the reconfiguration
process one can approximately retain the band gap formation corresponding to
the system. Alternatively, band gap formation can be also tailored by an
appropriate variation in the geometric parameters. All of these results
indicate that the proposed mechanical metamaterial can be very useful in a
plethora of applications ranging from smart protective to biomedical devices.

##  4 Methods

###  4.1 Fabrication

The Poisson’s ratio of the considered system was assessed based on micro-scale
experimental prototypes produced through the use of a commercial 3D printer
(Photonic Professional GT+, Nanoscribe GmbH). The 3D printer used in this work
operates based on the two-photon lithography method. All samples were printed
with a resolution of around 2 $\mu$m which was possible thanks to use of the
negative tone IP-S photoresin (Nanoscribe GmbH). Prior to the printing
process, a drop of resin was deposited on an ITO-coated soda-lime glass
substrate having the following dimensions: 25 × 25 × 0.7 mm3. After
deposition, specific fragments of the drop were photopolymerized by means of a
25X-objective with a femtosecond laser that was operating at $\lambda$ = 780
nm. Furthermore, some of the standard printing parameters were set to be the
following: a galvanometric scanning speed of 100 mm/s, a laser power at the
level of 90% and slicing as well as hatching distances were set to be equal to
1 $\mu$m and 0.5 $\mu$m respectively. Once the printing process was finished,
each sample was developed for 25 min in a solution of PGMEA (Propylene Glycol
Methyl Ether acetate). This step of the procedure allows to remove the
nonpolymerized photoresist. Subsequently, the resulting structure was rinsed
for a duration of three minutes in isopropyl alcohol (IPA) in order to remove
the developer solution.

In the case of the macroscopic experimental prototype, the structure was 3D
printed by means of a commercial FDM 3D-printer (Ultimaker) by means of the
elastic TPU 95A material. Furthermore, in each of the star-shaped segments of
the structure, there were two apertures with magnetic inclusions in the form
of 3 neodymium magnets of the N42 type. Each of the magnets had a cylindrical
shape with a height of 1 mm and a base diameter equal to 1 mm. The remanent
induction of each magnet was estimated to be equal to 0.613 T. The geometric
dimensions corresponding to the initial structure were set to be the
following: $a$ = 6 mm, $b$ = 9 mm, $l_{1}$ = 0.7 mm, $l_{2}$ = 0.43 mm,
$l_{3}$ = 0.36 mm, $l_{4}$ = 0.8 mm and $l_{5}$ = 0.55 mm.

###  4.2 Mechanical Compression Testing

The experimental samples considered in this work were tested mechanically in
order to assess their Poisson’s ratio. To this aim, all samples were subjected
to the mechanical compression induced by a flat external indenter. The speed
of the indenter was constant and was set to be equal to 1 $\mu$m s-1. At the
same time, the bottom part of each sample was fixed in space due to its
connection to the substrate on top of which they were printed. Throughout the
experiment, the deformation process was recorded by means of an optical camera
with magnification factor x20. The working distance, i.e. the distance between
the lens of the camera and the sample, was approximately equal to 11 mm. Based
on the recorded pictures, the Poisson’s ratio of the system was determined
based on the central topmost unit-cell (see Supplementary Information). To
track the position of vertices that could be used to calculate the Poisson’s
ratio of the selected unit-cell, the digital image correlation method
implemented in the MATLAB software was used.

In the case of the results presented in Fig. 5, the macroscopic experimental
prototype was inserted into the uniform external magnetic field that was
induced by two parallel electromagnets. The magnitude of the applied field was
varied in a range between 0 T and 0.24 T with steps of 0.01 T in order to
induce the reconfiguration process (see Supplementary Video 1). The
deformation of the structure was recorded by means of a standard optical
camera.

###  4.3 FEM simulations

To determine the Poisson’s ratio of the considered finite structurs resembling
experimental prototypes presented in Fig. 2, FEM simulations were conducted by
means of the commercial software COMSOL Multiphysics. To this aim, for each of
the structures considered in Fig. 2, it was assumed that geometric parameters
characterizing the structure are the same as in the case of the experimental
samples. Furthermore, to match the constraint used in the experiments, it was
assumed that the bottom part of each of the systems is fixed in space. In
addition, the deformation of the structure was induced by means of contact
with a large flat indenter in order to match the deformation process observed
in the experiments as closely as possible. Mechanical properties of the
simulated material were set as $E$ = 4 GPa, $\nu$ = 0.4 and $\rho$ = 1200 kg
m-3. In addition, geometrical nonlinearities were implemented to make
simulations more realistic under large strains. Finally, it should be
emphasized that a specific method related to the calculation of the Poisson’s
ratio is provided in the Supplementary Information.

In order to determine the phonon band structure, Bloch boundary conditions
were implemented and used via an eigenvalue solver assuming real wave vectors
in the Bloch conditions and searching for real eigenfrequencies of the system.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the
corresponding author upon request.

Acknowledgements

This work was partly supported by the french RENATECH network and its FEMTO-ST
technological facility.

K.K.D. acknowledges the support of the Polish National Science Centre (NCN) in
the form of the grant awarded as a part of the SONATINA 5 program, project No.
2021/40/C/ST5/00007 under the name ”Programmable magneto-mechanical
metamaterials guided by the magnetic field”.

This research was funded by the Polish Minister of Education and Science under
the program “Regional Initiative of Excellence” in 2019-2023, project No.
003/RID/2018/19, funding amount PLN 11 936 596.10.

## References

  * [1] B. Florijn, C. Coulais, and M. van Hecke, “Programmable mechanical metamaterials,” Phys. Rev. Lett., vol. 113,, p. 175503, 2014.
  * [2] C. Coulais, E. Teomy, K. de Reus, Y. Shokef, and M. van Hecke, “Combinatorial design of textured mechanical metamaterials,” Nature, vol. 535,, pp. 529–532, 2016.
  * [3] C. Coulais, A. Sabbadini, F. Vink, and M. van Hecke, “Multi-step self-guided pathways for shape-changing metamaterials,” Nature, vol. 561,, pp. 512–515, 2018.
  * [4] R. Fleury, D. L. Sounas, C. F. Fleck, M. R. Haberman, and A. Alú, “Sound isolation and giant linear nonreciprocity in a compact acoustic circulator,” Science, vol. 343,, pp. 516–519, 2014.
  * [5] R. Fleury, D. L. Sounas, and A. Alú, “An invisible acoustic sensor based on parity-time symmetry,” Nat. Commun., vol. 6,, p. 5905, 2015.
  * [6] F. Wenz, I. Schmidt, A. Leichner, T. Lichti, S. Baumann, H. Andrae, and C. Eberl, “Designing shape morphing behavior through local programming of mechanical metamaterials,” Adv. Mater., vol. 33,, p. 2008617, 2021.
  * [7] Y. Zhu, M. Birla, K. R. Oldham, and E. T. Filipov, “Elastically and plastically foldable electrothermal micro-origami for controllable and rapid shape morphing,” Adv. Funct. Mater., vol. 30,, p. 2003741, 2020.
  * [8] T.-Y. Huang, H.-W. Huang, D. D. Jin, Q. Y. Chen, J. Y. Huang, L. Zhang, and H. L. Duan, “Four-dimensional micro-building blocks,” Sci. Adv., vol. 6,, 2020.
  * [9] H. Tao and J. Gibert, “Multifunctional mechanical metamaterials with embedded triboelectric nanogenerators,” Adv. Funct. Mater., vol. 30,, p. 2001720, 2020.
  * [10] J. A. Jackson, M. C. Messner, N. A. Dudukovic, W. L. Smith, L. Bekker, B. Moran, A. M. Golobic, A. J. Pascall, E. B. Duoss, K. J. Loh, and C. M. Spadaccini, “Field responsive mechanical metamaterials,” Sci. Adv., vol. 4,, p. eaau6419, 2018.
  * [11] J. Qi, Z. Chen, P. Jiang, W. Hu, Y. Wang, Z. Zhao, X. Cao, S. Zhang, R. Tao, Y. Li, and D. Fang, “Recent progress in active mechanical metamaterials and construction principles,” Adv. Sci., p. 2102662, 2021.
  * [12] H. M. A. Kolken, S. Janbaz, S. M. A. Leeflang, K. Lietaert, H. H. Weinans, and A. A. Zadpoor, “Rationally designed meta-implants: A combination of auxetic and conventional meta-biomaterials.,” Mater. Horiz., vol. 5,, pp. 28–35, 2018.
  * [13] T. van Manen, S. Janbaz, K. M. B. Jansen, and A. A. Zadpoor, “4d printing of reconfigurable metamaterials anddevices,” Commun. Mater., vol. 2,, p. 56, 2021.
  * [14] M. Miniaci, A. Krushynska, F. Bosia, and N. M. Pugno, “Large scale mechanical metamaterials as seismic shields,” New J. Phys., vol. 8,, p. 083041, 2016\.
  * [15] G. Imbalzano, P. Tran, T. D. Ngo, and P. V. S. Lee, “A numerical study of auxetic composite panels under blast loadings,” Compos. Struct., vol. 135,, pp. 339–352, 2016.
  * [16] N. Novak, M. Borovinsek, O. Al-Ketan, Z. Ren, and M. Vesenjak, “Impact and blast resistance of uniform and graded sandwich panels with tpms cellular structures,” Compos. Struct., vol. 300,, p. 116174, 2022.
  * [17] O. Duncan, T. Shepherd, C. Moroney, L. Foster, P. D. Venkatraman, K. Winwood, T. Allen, and A. Alderson, “Review of auxetic materials for sports applications: Expanding options in comfort and protection,” Appl. Sci., vol. 8,, p. 941, 2018.
  * [18] A. Russo, B. Y. Ahn, J. J. Adams, E. B. Duoss, J. T. Bernhard, and J. A. Lewis, “Pen-on-paper flexible electronics,” Adv. Mater., vol. 23,, pp. 3426–3430, 2011.
  * [19] S. Huang, Y. Liu, Y. Zhao, Z. Ren, and C. F. Guo, “Flexible electronics: Stretchable electrodes and their future,” Adv. Funct. Mater., vol. 29,, p. 1805924, 2019.
  * [20] R. Lakes, “Foam structures with a negative poisson’s ratio,” Science, vol. 235,, p. 1038, 1987.
  * [21] K. E. Evans and A. Alderson, “Auxetic materials: Functional materials and structures from lateral thinking!,” Adv. Mater., vol. 12,, pp. 617–628, 2000.
  * [22] A. S. Gladman, E. A. Matsumoto, R. G. Nuzzo, L. Mahadevan, and J. A. Lewis, “Biomimetic 4d printing,” Nat. Mater., vol. 15,, pp. 413–418, 2016.
  * [23] J. T. B. Overvelde, J. C. Weaver, C. Hoberman, and K. Bertoldi, “Rational design of reconfigurable prismatic architected materials,” Nature, vol. 541,, pp. 347–352, 2017.
  * [24] M. J. Mirzaali, S. Janbaz, M. Strano, L. Vergani, and A. A. Zadpoor, “Shape-matching soft mechanical metamaterials,” Sci. Rep., vol. 8,, p. 965, 2018.
  * [25] J. L. Silverberg, A. A. Evans, L. McLeod, R. C. Hayward, T. Hull, C. D. Santangelo, and I. Cohen, “Using origami design principles to fold reprogrammable mechanical metamaterials,” Science, vol. 345,, pp. 647–650, 2014.
  * [26] A. Rafsanjani, K. Bertoldi, and A. R. Studart, “Programming soft robots with flexible mechanical metamaterials,” Sci. Robot., vol. 4,, p. eaav7874, 2019\.
  * [27] P. Cai, C. Wang, H. Gao, and X. Chen, “Mechanomaterials: A rational deployment of forces and geometries in programming functional materials,” Adv. Mater., vol. 33,, p. 2007977, 2021.
  * [28] C. Song, B. Zou, Z. Cui, Z. Liang, and J. Ju, “Thermomechanically triggered reversible multi-transformability of a single material system by energy swapping and shape memory effects,” Adv. Funct. Mater., vol. 31,, p. 2101395, 2021.
  * [29] S. Li, G. Librandi, Y. Yao, A. J. Richard, A. Schneider-Yamamura, J. Aizenberg, and K. Bertoldi, “Controlling liquid crystal orientations for programmable anisotropic transformations in cellular microstructures,” Adv. Mater., vol. 33,, p. 2105024, 2021.
  * [30] K. Bertoldi, V. Vitelli, J. Christensen, and M. van Hecke, “Flexible mechanical metamaterials,” Nat. Rev. Mater., vol. 2,, p. 17066, 2017.
  * [31] J. N. Grima, R. Gatt, B. Ellul, and E. Chetcuti, “Auxetic behaviour in non-crystalline materials having star or triangular shaped perforations,” J. Non. Cryst. Solids, vol. 356,, pp. 1980–1987, 2010.
  * [32] R. M. Neville, F. Scarpa, and A. Pirrera, “Shape morphing kirigami mechanical metamaterials,” Sci. Rep., vol. 6,, p. 31067, 2016.
  * [33] X. Chen, J. Moughames, Q. Ji, J. A. I. Martínez, H. Tan, S. Adrar, N. Laforge, J.-M. Cote, S. Euphrasie, G. Ulliac, et al., “Optimal isotropic, reusable truss lattice material with near-zero poisson’s ratio,” Extreme Mechanics Letters, vol. 41, p. 101048, 2020.
  * [34] X. Chen, Q. Ji, J. A. I. Martinez, H. Tan, G. Ulliac, V. Laude, and M. Kadic, “Closed tubular mechanical metamaterial as lightweight load-bearing structure and energy absorber,” Journal of the Mechanics and Physics of Solids, vol. 167, p. 104957, 2022.
  * [35] X. Tan, J. A. I. Martínez, G. Ulliac, B. Wang, L. Wu, J. Moughames, M. Raschetti, V. Laude, and M. Kadic, “Single-step-lithography micro-stepper based on frictional contact and chiral metamaterial,” Small, vol. 18, no. 28, p. 2202128, 2022.
  * [36] X. Chen, J. Moughames, Q. Ji, J. A. I. Martínez, H. Tan, G. Ulliac, V. Laude, and M. Kadic, “3d lightweight mechanical metamaterial with nearly isotropic inelastic large deformation response,” Journal of the Mechanics and Physics of Solids, vol. 169, p. 105057, 2022.
  * [37] B. Deng, A. Zareei, X. Ding, J. C. Weaver, C. H. Rycroft, and K. Bertoldi, “Inverse design of mechanical metamaterials with target nonlinear response via a neural accelerated evolution strategy,” Adv. Mater., vol. 34,, p. 2206238, 2022.
  * [38] K. W. Wojciechowski, “Two-dimensional isotropic system with a negative poisson ratio,” Phys. Lett. A, vol. 137,, pp. 60–64, 1989.
  * [39] L. Mizzi, K. M. Azzopardi, D. Attard, J. N. Grima, and R. Gatt, “Auxetic metamaterials exhibiting giant negative poisson’s ratios,” Phys. Status Solidi RRL, vol. 9,, pp. 425–430, 2015.
  * [40] Y.-L. Wei, Q.-S. Yang, and R. Tao, “Smp-based chiral auxetic mechanical metamaterial with tunable bandgap function,” Int. J. Mech. Sci., vol. 195,, p. 106267, 2021.
  * [41] L. Wang, G. Ulliac, B. Wang, J. A. I. Martinez, K. K. Dudek, V. Laude, and M. Kadic, “3d auxetic metamaterials with elastically-stable continuous phase transition,” Adv. Sci., p. 2204721, 2022.
  * [42] J. Qu, M. Kadic, and M. Wegener, “Poroelastic metamaterials with negative effective static compressibility,” Applied Physics Letters, vol. 110, no. 17, p. 171901, 2017.
  * [43] T. A. M. Hewage, K. L. Alderson, A. Alderson, and F. Scarpa, “Double-negative mechanical metamaterials displaying simultaneous negative stiffness and negative poisson’s ratio properties,” Adv. Mater., vol. 28,, pp. 10323–10332, 2016.
  * [44] K. K. Dudek, R. Gatt, and J. N. Grima, “3d composite metamaterial with magnetic inclusions exhibiting negative stiffness and auxetic behaviour,” Mater. Des., vol. 187,, p. 108403, 2020.
  * [45] X. Tan, L. Wang, S. Zhu, S. Chen, B. Wang, and M. Kadic, “A general strategy for performance enhancement of negative stiffness mechanical metamaterials,” Eur. J. Mech. A Solids, vol. 96,, p. 104702, 2022.
  * [46] J. Qu, M. Kadic, and M. Wegener, “Three-dimensional poroelastic metamaterials with extremely negative or positive effective static volume compressibility,” Extreme Mechanics Letters, vol. 22, pp. 165–171, 2018\.
  * [47] Z. G. Nicolaou and A. E. Motter, “Mechanical metamaterials with negative compressibility transitions,” Nat. Mater., vol. 11,, pp. 608–613, 2012\.
  * [48] R. H. Baughman, S. Stafstrom, C. Cui, and S. O. Dantas, “Materials with negative compressibilities in one or more dimensions,” Science, vol. 279,, pp. 1522–1524, 1998.
  * [49] T. Li, F. Liu, and L. Wang, “Enhancing indentation and impact resistance in auxetic composite materials,” Compos. B. Eng., vol. 198,, p. 108229, 2020\.
  * [50] S. Yuan, C. K. Chua, and K. Zhou, “3d-printed mechanical metamaterials with high energy absorption,” Adv. Mater. Technol., vol. 4,, p. 1800419, 2019\.
  * [51] M. Mohsenizadeh, F. Gasbarri, M. Munther, A. Beheshti, and K. Davami, “Additively-manufactured lightweight metamaterials for energy absorption,” Mater. Des., vol. 139,, pp. 521–530, 2018.
  * [52] A. O. Krushynska, M. Miniaci, F. Bosia, and N. M. Pugno, “Coupling local resonance with bragg band gaps in single-phase mechanical metamaterials,” Extreme Mech. Lett., vol. 12,, pp. 30–36, 2017.
  * [53] J. N. Grima, R. Caruana-Gauci, M. R. Dudek, K. W. Wojciechowski, and R. Gatt, “Smart metamaterials with tunable auxetic and other properties,” Smart Mater. Struct., vol. 22,, p. 084016, 2013.
  * [54] R. Galea, K. K. Dudek, P.-S. Farrugia, L. Z. Mangion, J. N. Grima, and R. Gatt, “Reconfigurable magneto-mechanical metamaterials guided by magnetic fields,” Compos. Struct., vol. 280,, p. 114921, 2022.
  * [55] M. Pishvar and R. L. Harne, “Foundations for soft, smart matter by active mechanical metamaterials,” Adv. Sci., vol. 7,, p. 2001384, 2020.
  * [56] S. M. Montgomery, S. Wu, X. Kuang, C. D. Armstrong, C. Zemelka, Q. Ze, R. Zhang, R. Zhao, and H. H. J. Qi, “Magneto-mechanical metamaterials with widely tunable mechanical properties and acoustic bandgaps,” Adv. Funct. Mater., vol. 31,, p. 2005319, 2021.
  * [57] C. Ma, Y. Chang, S. Wu, and R. R. Zhao, “Deep learning-accelerated designs of tunable magneto-mechanical metamaterials,” ACS Appl. Mater. Interfaces, vol. 14,, pp. 33892–33902, 2022.
  * [58] Q. Ji, J. Moughames, X. Chen, G. Fang, J. J. Huaroto, V. Laude, J. A. I. Martinez, G. Ulliac, C. Clevy, P. Lutz, K. Rabenorosoa, V. Guelpa, A. Spangenberg, J. Liang, A. Mosset, and M. Kadic, “4d thermomechanical metamaterials for soft microrobotics,” Commun. Mater., vol. 2,, p. 93, 2021\.
  * [59] R. Gatt, L. Mizzi, J. I. Azzopardi, K. M. Azzopardi, D. Attard, A. Casha, J. Briffa, and J. N. Grima, “Hierarchical auxetic mechanical metamaterials,” Sci. Rep., vol. 5,, p. 8395, 2015.
  * [60] Y. Cho, J.-H. Shin, A. Costa, V. K. T. A. Kim, J. Li, S. Y. Lee, S. Yang, H. N. Han, I.-S. Choi, and D. J. Srolovitz, “Engineering the shape and structure of materials by fractal cut,” Proc. Natl. Acad. Sci., vol. 111,, pp. 17390–17395, 2014.
  * [61] Y. Tang, G. Lin, L. Han, S. Qiu, S. Yang, and J. Yin, “Design of hierarchically cut hinges for highly stretchable and reconfigurable metamaterials with enhanced strength,” Adv. Mater., vol. 27,, pp. 7181–7190, 2015.
  * [62] N. An, A. G. Domel, J. Zhou, A. Rafsanjani, and K. Bertoldi, “Programmable hierarchical kirigami,” Adv. Funct. Mater., vol. 30,, p. 1906711, 2020\.
  * [63] K. Billon, I. Zampetakis, F. Scarpa, M. Ouisse, E. Sadoulet-Reboul, M. Collet, A. Perriman, and A. Hetherington, “Mechanics and band gaps in hierarchical auxetic rectangular perforated composite metamaterials,” Compos. Struct., vol. 160,, pp. 1042–1050, 2017.
  * [64] K. K. Dudek, J. A. I. Martínez, G. Ulliac, and M. Kadic, “Micro-scale auxetic hierarchical mechanical metamaterials for shape morphing,” Adv. Mater., vol. 34,, p. 2110115, 2022.
  * [65] K. K. Dudek, J. A. I. Martínez, and M. Kadic, “Variable dual auxeticity of the hierarchical mechanical metamaterial composed of re-entrant structural motifs,” Phys. Status Solidi B, p. 2200404, 2022.
  * [66] T. Frenzel, J. Köpfler, E. Jung, M. Kadic, and M. Wegener, “Ultrasound experiments on acoustical activity in chiral mechanical metamaterials,” Nat. Commun., vol. 10,, p. 3384, 2019.
  * [67] C. Mei, L. Li, H. Tang, X. Han, X. Wang, and Y. Hu, “Broadening band gaps of shear horizontal waves of metamaterials via graded hierarchical architectures,” Compos. Struct., vol. 271,, p. 114118, 2021.
  * [68] V. Kunin, S. Yang, Y. Cho, P. Deymier, and D. J. Srolovitz, “Static and dynamic elastic properties of fractal-cut materials,” Extreme Mech. Lett., vol. 6,, pp. 103–114, 2016.
  * [69] H.-W. Huang, M. S. Sakar, A. J. Petruska, S. Pane, and B. J. Nelson, “Soft micromachines with programmable motility and morphology,” Nat. Commun., vol. 7, p. 12263, 2016.
  * [70] J. Kim, S. E. Chung, S.-E. Choi, H. Lee, J. Kim, and S. Kwon, “Programming magnetic anisotropy in polymeric microactuators,” Nature Mater., vol. 10, pp. 747–752, 2011.
  * [71] T.-Y. Huang, M. S. Sakar, A. Mao, A. J. Petruska, F. Qiu, X.-B. Chen, S. Kennedy, D. Mooney, and B. J. Nelson, “3d printed microtransporters: Compound micromachines for spatiotemporally controlled delivery of therapeutic agents,” Adv. Mater., vol. 27, pp. 6644–6650, 2015.
  * [72] J. Cui, T.-Y. Huang, Z. Luo, P. Testa, H. Gu, X.-Z. Chen, B. J. Nelson, and L. J. Heyderman, “Nanomagnetic encoding of shape-morphing micromachines,” Nature, vol. 575, pp. 164–168, 2019.


# Rao-Burbea centroids applied to the statistical characterisation of time
series and images through ordinal patterns

Diego M. Mateos ID Consejo Nacional de Investigaciones Científicas y Técnicas
(CONICET), Argentina.  Facultad de Ciencia y Tecnología. Universidad Autónoma
de Entre Ríos (UADER). Oro Verde, Entre Ríos, Argentina.  Instituto de
Matemática Aplicada del Litoral (IMAL-CONICET-UNL), CCT CONICET, Santa Fé,
Argentina.  Corresponding author: Diego M. Mateos<EMAIL_ADDRESS>Leonardo E. Riveaud  Facultad de Ingeniería, Universidad Nacional del Comahue
(FAIN, UNComa)  Pedro W. Lamberti  Consejo Nacional de Investigaciones
Científicas y Técnicas (CONICET), Argentina.  Facultad de Matemática
Astronomía, Física y Computación (FaMAF), Universidad Nacional de Córdoba.
Córdoba, Argentina.

###### Abstract

Divergences or similarity measures between probability distributions have
become a very useful tool for studying different aspects of statistical
objects such as time series, networks and images. Notably not every divergence
provides identical results when applied to the same problem. Therefore it is
convenient to have the widest possible set of divergences to be applied to the
problems under study. Besides this choice an essential step in the analysis of
every statistical object is the mapping of each one of their representing
values into an alphabet of symbols conveniently chosen. In this work we attack
both problems, that is, the choice of a family of divergences and the way to
do the map into a symbolic sequence. For advancing in the first task we work
with the family of divergences known as the Burbea-Rao centroids (BRC) and for
the second one we proceed by mapping the original object into a symbolic
sequence through the use of ordinal patterns. Finally we apply our proposals
to analyse simulated and real time series and to real textured images. The
main conclusion of our work is that the best BRC, at least in the studied
cases, is the Jensen Shannon divergence, besides the fact that it verifies
some interesting formal properties.

Keywords— Divergences, Ordinal patterns, Time series analysis, Textured Images

The notion of distinguishability is of crucial importance in probability
theory and in general in every statistical theory. Due to statistical
fluctuations it becomes difficult to distinguish between close probability
distributions. This problem involves two issues. The first one is the
definition of distance measures between probability distributions and the
second one is the introduction of distinguishability criteria. On the other
hand, many physical phenomena can be represented by means of a time series,
which could be generated by chaotic dynamics or by a random proccess. It is
necessary, in order to be able to study these dynamics, to make a
correspondence between the values of the time series under study with a
symbolic sequence. In several occasions this task is the most difficult to
perform. In this work we address the different problems indicated above.
First, we study the applicability of a family of divergences between
probability distributions, known as Burbea-Rao centroids, to the study of time
series. We then investigate the use of ordinal patterns as a way to assign a
symbolic sequence to a time series. With these two tools we study the
statistical characteristics of simulated and real world time series as well as
to textured 2D images. As the Burbea-Rao centroids depend on a concave
functional we test the behaviour of this divergences when we change it. Among
them one of the choosen funtional leads to the Jensen Shannon divergence. In
addition of demonstrating the robustness of the scheme presented it was
possible to conclude that in all the investigated examples the best results
correspond to the Jensen Shannon divergence.

##  1 Introduction

There exists a profusion of divergences between probability distributions.
Remarkably when some of them are applied to the same statistical problem, in
general they do not lead to indistinguishable results. Therefore, might be
useful to have a large set of divergences. In general the divergences have
different origins. Some are statistical, others originated in information
theory and others borrowed from the realm of pure mathematics. The Fisher’s
metric is a conspicuous example of the first kind [1]; and the Kullback-
Leibler is a well known example originated in information theory [2]. Measures
of similarity (or dissimilarity), between probability distributions have
become of great interest in physics (classical and quantum), biology and many
other areas of science and technology [3, 4, 5, 6].

In the realm of statistical physics particularly notable has been the use of
divergences in the context of out-of-equilibrium systems. The obligatory
reference on this subject is the work of G. Crooks and D. Sivak. There these
authors give a physical interpretation to different measures of divergence
when are applied to the study of conjugate ensembles of non equilibrium
trajectories. The relative entropy is related to the dissipation, the Jeffreys
divergence is the average dissipation of the forward and reverse evolution
(hysteresis), the Jensen – Shannon divergence has been proposed as a magnitude
of the time arrow, and the Chernoff divergence is the work cumulant generating
function [7].

As it was said before having a variety of divergence could be useful for both
pure theoretical studies and in the context of applications. Sometimes it is
possible to introduce families of divergences, labelling each member with a
parameter [8, 9] or by giving a general structure depending on certain
functional with adequate characteristics. An example of this last case are the
Csiszar’s divergences or $f$-divergences.

Let $\Omega=\\{\omega_{1},...,\omega_{N+1}\\}$ be a discrete sample space. The
set of probability distributions on $\Omega$ can be identified with the
simplex $\mathbb{P}^{N}=\\{x^{i}\geq
0,i=1,...,N+1,\sum_{i=1}^{N+1}x^{i}=1\\}$. Given two probability distributions
$P=\\{p_{i}\\}_{i=1}^{N+1}$ and $Q=\\{q_{i}\\}_{i=1}^{N+1}$ belonging to
$\mathbb{P}^{N}$, a Csiszar’s divergence is defined by [10]

$\mathcal{D}_{f}(P,Q)=\sum_{i=1}^{N+1}f\left(\frac{p_{i}}{q_{i}}\right)q_{i}$

where $f(x)$ is a convex function such that $f(1)=0$. This family has been
extensively studied in the context of information geometry [11]. A remarkable
result is that when $p_{i}$ and $q_{i}=p_{i}+\delta p_{i}$ are two close
probability distributions, the divergence $\mathcal{D}_{f}(P,Q)$ is
proportional to the Fisher (Riemannian) metric:

$\mathcal{D}_{f}(P,P+\delta
P)\sim\frac{f^{\prime\prime}(1)}{2}\sum_{i}\frac{\delta p_{i}^{2}}{p_{i}}$ (1)

The above mentioned Kullback-Leibler divergence corresponds to $f(t)=t~{}\log
t$ and it is not symmetric and the Jensen-Shannon divergence (JSD) it is
obtained by taking $f(t)=(t+1)\log(\frac{2}{t+1})+t~{}\log t$, which clearly
is symmetric [12].

For the study of dynamical phenomena, we need to have a sequence of
measurements related to them. These sequences are usually given in the form of
time series which allow extracting information on the underlying physical
system under study. Something similar happen with images, where we have a two-
dimensional matrix which contains the information about each pixel of the
images. The time series or images can be associated with a probability
distribution function (PDF) and, from this to have a way for applying the
divergences. To this aim it is necessary to map the time series or images in a
finite alphabet. There exist many methods to discretize a continuous time
series, such as binarization or the wavelet analysis. However these methods do
not take into account the relation between the value for a given time and the
neighbouring values. In 2002 C. Bandt and B. Pompe introduced the ordinal
patterns for time series discretization [13]. This approach has been
extensively used to different problems with very relevant results. Its main
advantage is that reveals the underlying dynamics in the process that
generates the time series. For a review of these properties we suggest the
work of M. Zanin and F. Olivares in which the ordinal pattern approach is
studied in detail [14].

Here we use the Burbea-Rao centroids (BRC) as the dissimilarity measures and
we build a symbolic alphabet formed by ordinal patterns, in which the original
time series or the 2D images are mapped. Then we apply an segmentation
procedure meaning this a way to detect dynamical changes in the corresponding
symbolic sequence. All this scheme is used to the study of the statistical
properties of simulated time series and to real 2D images.

This manuscript is organised as follows. In section 2 we review some
properties of the BRC. We show that the BRC can be thought as a deformations
of the Euclidean metric and we indicate a particular character of the JSD as a
BRC. In section 3 we give a brief explanation of the ordinal patterns and how
them can be evaluated in some concrete cases. In section 4 we combine the use
of the BRC, really a particular case that we called the $\gamma$-divergences,
with the ordinal patterns to four examples, three related to time series
analysis and one to images analysis. Finally, in section 5 we discussed our
results as well as some other possible frameworks of applicability of our
scheme.

##  2 Burbea-Rao centroids

In 1982 J. Burbea and C.R. Rao introduced a family of divergences now known as
the Burbea-Rao centroids (BRC). They started from a concave $\Phi$ functional.
A generic member of this family is given by:

$\mathcal{J}_{\Phi}(P,Q)=\frac{1}{2}[\Phi(P)+\Phi(Q)]-\Phi\left(\frac{P+Q}{2}\right)$
(2)

where $P$ and $Q$ are two discrete probability distributions belongings to
$\mathbb{P}^{N}$.

Incidentally it is worth mentioning that a centroid can be interpreted as a
deformation of the square of the Euclidean metric. Indeed, let us write the
square of the Euclidean distance between $P$ and $Q$ in the form:

$\mathcal{E}(P,Q)=\mathop{\sum}_{i}(q_{i}-p_{i})^{2}=\mathop{\sum}_{i}\left(2~{}q_{i}^{2}+2~{}p_{i}^{2}-(p_{i}+q_{i})^{2}\right)$
(3)

which can be rewritten as a BRC with $\Phi=\sum_{i}x_{i}g(x_{i})$ and
$g(x)=x$. In this sense we can think that a BRC is a “deformation” of the
Euclidean metric through the map $\sum_{i}x_{i}^{2}\rightarrow\Phi$.
Analogously the Jensen-Shannon divergence (JSD) mentioned above can be written
as a BRC with $\Phi=\sum_{i}x_{i}g(x_{i})$ and $g(x)=\frac{1}{2}\log x$. The
explicit expression of the JSD is:

$\mathcal{D}_{JS}(P,Q)=\sum_{i}\frac{1}{2}p_{i}\log p_{i}+\frac{1}{2}q_{i}\log
q_{i}-\frac{1}{2}(p_{i}+q_{i})\log\left(\frac{p_{i}+q_{i}}{2}\right)$ (4)

An important characteristic of the JSD is that it is the only BRC that is a
Csiszar divergence. This can be easily proved by doing a Taylor expansion of
(2) up to second order. Not every BRC is a metric in the sense that verifies
the triangle inequality. However it has been shown that the square root of the
JSD is a true metric for the simplex $\mathbb{P}^{N}$ as well as the Euclidean
metric (3). In general the metric character of a BRC will depend on the
functional $\Phi$.

One advantage of a BRC that is particularly interesting for us is that by
modifying the function $\Phi$, we can determine the divergence that more
adequately allows to characterise the statistical object under study.

###  2.1 Weighted BRC

In 1991 J. Lin introduced a generalised version of the JSD by assigning
different weights to each probability distribution, $P$ and $Q$. Let $\pi_{P}$
and $\pi_{Q}$ two non negative numbers such that $\pi_{P}+\pi_{Q}=1$. Then the
weighted JSD is given by:

$\mathcal{D}_{JS}(P,Q;\pi_{P},\pi_{Q})=H_{S}(\pi_{P}P+\pi_{Q}Q)-\pi_{P}H_{S}(P)-\pi_{Q}H_{S}(Q)$
(5)

with $H_{S}(P)\equiv-\sum_{i}p_{i}\log(p_{i})$ the Shannon entropy. We can
proceed in the same way with a generic BRC and define the weighted BRC as
follows:

$\mathcal{J}_{\Phi}(P,Q,\pi_{P},\pi_{Q})\equiv\pi_{P}\Phi(P)+\pi_{Q}\Phi(Q)-\Phi(\pi_{P}P+\pi_{Q}Q)$
(6)

Due to the functional $\Phi$ it is assumed to be concave, as a consequence of
the Jensen inequality, the quantity (6) is always non negative.

Another possible generalisation of the BRC is to evaluate it among several
probability distributions each one with different weight. To fix ideas let us
suppose $K$ probability distributions $P^{(1)},...,P^{[K)}$ and weights
$\pi_{1},...\pi_{K}$. We can define a weighted RBC among these probability
distributions in the form:

$\mathcal{J}_{\Phi}(P^{(1)},...,P^{(K)};\pi_{1},...,\pi_{N})\equiv\sum_{k=1}^{K}\pi_{k}\Phi(P^{(k)})-\Phi\left(\sum_{k=1}^{K}\pi_{k}P^{(k)}\right)$
(7)

Again due to the concavity of $\Phi$ and the Jensen inequality this last
quantity is non negative.

##  3 Ordinal Patterns and Burbea Rao centroids

As was remarked in the introduction for studying the statistical properties of
a time series, and with the purpose of assigning probability distributions to
it, it is indispensable to map the original data set, that is, the values of
the time series, into symbolic sequences. There are different methods in the
literature for doing this, such as the discretization, binarization and
wavelet’s histograms, among others. An alternative method proposed by Bandt
and Pompe (BP) [13] consists in transforming the time series, via a non
parametric transformation, into a sequence of patterns and then making
inference over these patterns. With this, the analysis gains robustness and
becomes apt to detect relevant causal information related to the unobserved
variables that control the underlying dynamics of the system. The BP approach
is based on the computation of the Shannon entropy from the histogram of
causal patterns. Later the ideas of BP were extended in different ways, being
one particularly relevant for us the the application to the analysis of two
dimensional images [15, 16].

For a time series $\mathcal{X}(t)$ the ordinal patterns procedure is based on
the relative values of the neighbours belonging to the series, and in
consequence takes into account the time structure or causality of the process
that generated the sequence. To understand this idea, let us consider a real-
valued discrete-time series $\mathcal{X}(t)=\\{x_{t}\in\mathbb{R}\\}$, and the
parameters $d\geq 2$ ( embedding dimension) and $\tau\geq 1$ (the time delay).
From these we construct a $d$-dimensional vector $\mathbf{Y}_{t}$:

$\mathbf{Y}_{t}=(x_{t-(d-1)\tau},\dots,x_{t-\tau},x_{t})~{}\qquad~{}{\mathrm{with}}~{}\qquad~{}t\geq(d-1)\tau\
.$ (8)

The dynamical properties of the overall system is preserved by the vector
$\mathbf{Y}_{t}$. We sorted in ascending order the components of the phase
space trajectory $\mathbf{Y}_{t}$. Then, we can define a permutation vector,
$\mathbf{\Pi}_{t}$, with components given by the original position of the
sorted values. Each one of these vectors represents a pattern (or motif) with
$N_{\Pi}=d!$ possible patterns. For example let us consider the time series
$\mathcal{X}(t)=\\{0.42,~{}2.7,~{}4.2,~{}0.35,~{}1.5\\}$ and take the
parameters $d=3$ and $\tau=1$. We define the embedding vectors
$\mathbf{Y}_{t}$ as $\mathbf{Y}_{1}=(0.42,~{}2.7,~{}4.2)$;
$\mathbf{Y}_{2}=(2.7,~{}4.2,~{}0.35)$; $\mathbf{Y}_{3}=(4.2,~{}0.35,~{}1.5)$,
and the respective permutation vectors as $\mathbf{\Pi}_{1}=(0,~{}1,~{}2)$,
$\mathbf{\Pi}_{2}=(1,~{}2,~{}0)$ and $\mathbf{\Pi}_{3}=(2,~{}0,~{}1)$.
Regarding the selection of the parameters, Bandt and Pompe [13] suggested
working with $3\leq d\leq 6$ and specifically consider an embedding delay
$\tau=1$. Nevertheless, other values of $\tau$ could provide additional
information. It has been recently shown that this parameter is strongly
related to the intrinsic time scales of the system under analysis [17, 18]. In
addition, it is important to take into account that for a reliable estimate of
the ordinal probability distribution we need that, the length of the sequence
$N_{\mathcal{X}}$ must to be greater than the number of possible ordinal
patterns: $N_{\Pi}<<N_{\mathcal{X}}$ [13].

Ribeiro and Zunnino extended the permutation entropy framework to two-
dimensional data [15, 16]. We can consider a two-dimensional data array
$\\{x_{u}^{v}\\}$ with $u=1,...,N_{x}$ and $v=1,...,N_{y}$ whose elements
could be consider as pixels of an image. We define the embedding dimensions of
the horizontal ($d_{x}$) and vertical directions ($d_{y}$), and the
corresponding $\tau_{x}$ and $\tau_{y}$. Using a similar scheme that in the
time series case, we slice the data array in partitions of size $d_{x}\times
d_{y}$ defined by:

$\mathbf{Y}_{i}^{j}=\begin{bmatrix}x_{i}^{j}&x_{i}^{j+\tau_{y}}&\ldots&x_{i}^{j+(d_{y}-1)\tau_{y}}\\\
x_{i+\tau_{x}}^{j}&x_{i+\tau_{x}}^{j+\tau_{y}}&\ldots&x_{i+\tau_{x}}^{j+(d_{y}-1)\tau_{y}}\\\
\vdots&\vdots&\ddots&\\\
x_{i+(d_{x}-1)\tau_{x}}^{j}&x_{i+(d_{x}-1)\tau_{x}}^{j+\tau_{y}}&\ldots&x_{i+(d_{x}-1)\tau_{x}}^{j+(d_{y}-1)\tau_{y}}\end{bmatrix}$
(9)

where $i=1,\ldots,n_{x}$ and $j=1,\ldots,n_{y}$ , with
$n_{x}=N_{x}-(d_{x}-1)\tau_{x}$ and $n_{y}=N_{y}-(d_{y}-1)\tau_{y}$.

To associate a permutation symbol with each two-dimensional partition we
concatenate line by line the partition $y_{i}^{j}$

$\begin{split}\mathbf{Y}_{i}^{j}=&(x_{i}^{j},x_{i}^{j+\tau_{y}},\ldots,x_{i}^{j+(d_{y}-1)\tau_{y}},\\\
&x_{i+\tau_{x}}^{j},x_{i+\tau_{x}}^{j+\tau_{y}},\ldots,x_{i+\tau_{x}}^{j+(d_{y}-1)\tau_{y}}\\\
&x_{i+(d_{x}-1)\tau_{x}}^{j},x_{i+(d_{x}-1)\tau_{x}}^{j+\tau_{y}},\ldots,x_{i+(d_{x}-1)\tau_{x}}^{j+(d_{y}-1)\tau_{y}}.)\end{split}$
(10)

Then we evaluate the permutation symbol associated with each data partition as
in the one-dimensional case to define the symbolic array
${\mathbf{\Pi}_{i}^{j}}$ for $i=1,\ldots,n_{x}$ and $j=1,\ldots,n_{y}$ related
to the data set. In this case the number of possible ordinal patterns are
$N_{\Pi}=(d_{x}d_{y})!$ [15]. For a better comprehension of the method we will
present a simple example. Let us suppose we have the following $3\times 3$
array:

$A=\begin{bmatrix}2&3&7\\\ 4&5&6\\\ 1&7&8\end{bmatrix}$

For this case we take $d_{x}=d_{y}=2$ and $\tau_{x}=\tau_{y}=1$ obtaining four
partitions: $\mathbf{Y}_{1}=\left[\begin{smallmatrix}2&3\\\
4&5\end{smallmatrix}\right]=(2,3,4,5)$ with ordinal pattern
$\mathbf{\Pi}_{1}=(0,1,2,3)$; $\mathbf{Y}_{2}=\left[\begin{smallmatrix}3&7\\\
5&6\end{smallmatrix}\right]=(3,7,5,6)$ with ordinal pattern
$\mathbf{\Pi}_{2}=(0,3,1,2)$; $\mathbf{Y}_{3}=\left[\begin{smallmatrix}4&5\\\
1&7\end{smallmatrix}\right]=(4,5,1,7)$ with ordinal pattern
$\mathbf{\Pi}_{3}=(1,2,0,3)$ and
$\mathbf{Y}_{4}=\left[\begin{smallmatrix}5&6\\\
7&8\end{smallmatrix}\right]=(5,6,7,8)$ with ordinal pattern
$\mathbf{\Pi}_{4}=(0,1,2,3)$.

Once we have all the ordinal patters $\mathbf{\Pi}$ corresponding to a signal
or a image we compute the probability distribution $P_{op}=P(\Pi_{i})$ with
$i=1,\ldots,d!$ for one-dimensional sequence or $i=1,\dots,(d_{x}d_{y})!$ for
two-dimensional arrays.

We will compare two signals or images by evaluating a divergence between the
corresponding probability distributions for the associated ordinal patterns,
$P_{op}$ and $Q_{op}$. We used the BRC, that we name the $\gamma$-divergence,
explicitly given by:

$\mathcal{D}_{\gamma}(P_{op}||Q_{op})=2\mathop{\sum}_{i=1}^{N_{\Pi}}\gamma_{g}(p(\Pi_{i}),q(\Pi_{i}))$
(11)

with

$\gamma_{g}(p(\Pi_{i}),q(\Pi_{i}))=\frac{1}{2}p(\Pi_{i})~{}g(p(\Pi_{i}))+\frac{1}{2}q(\Pi_{i})~{}g(q(\Pi_{i}))-\frac{1}{2}(p(\Pi_{i})+q(\Pi_{i}))~{}g\left(\frac{p(\Pi_{i})+q(\Pi_{i})}{2}\right)$
(12)

The only requirement on the function $g(x)$ is that the sum
$\sum_{i}x_{i}g(x_{i})$ results concave.

##  4 Applications

To investigate the behaviour of the $\gamma$-family we use different
$\gamma$-functions, $xg(x)$, by changing the $g$ function. We apply them to
four cases. In the first example we analyse the coupling of the Henon-Henon
system; in the second one we detect the change in the dynamics between chaotic
and stochastic signals; in the third example we apply our scheme to detect the
changes in the stages of an EEG signal and finally we evaluate the distances
between textured 2D images. In each one of this examples, we work with four
$g(x)$ function; i) $g(x)=e^{x}$; $g(x)=~{}log(x)$; $g(x)=~{}\sqrt[]{x}$ and
$g(x)=~{}sinh(x)$.

###  4.1 Henon-Henon coupled system

The coupled Henon-Henon system is described by the following set of equations
[19, 20]:

$\displaystyle y_{1}(n+1)$ $\displaystyle=1.4-y_{1}^{2}(n)+b~{}y_{2}(n)$
$\displaystyle y_{2}(n+1)$ $\displaystyle=y_{1}(n)$ $\displaystyle x_{1}(n+1)$
$\displaystyle=1.4-(\epsilon~{}y_{1}(n)+(1-\epsilon)~{}x_{1}(n)+b~{}x_{2}(n))$
$\displaystyle x_{2}(n+1)$ $\displaystyle=x_{1}(n)$

where, in our choice $b=0.3$. For the simulation the coupling parameter
$\epsilon$ was varied from zero to $1$ in steps of $0.1$. For each $\epsilon$,
$200$ realisations were computed using random initial conditions. The signals
length were $N=100000$ and the ordinal patterns were evaluated with the
parameters $d=3,4,5$ and $\tau=1$.

The Figure 1 shows the results for the four divergences as a function of the
coupling parameter $\epsilon$. The median value of estimator
$\mathcal{D}_{\gamma}$ is maximum for $\epsilon=0$ in all the cases. This is
expected because there is no information sharing between the two systems. The
$\mathcal{D}_{\gamma}$ values decrease with the coupling parameter up to the
value $\epsilon=0.6$. In contrast, for $\epsilon\geq 0.7$ the
$\mathcal{D}_{\gamma}$ values are zero. For these values of the coupling
parameter, the Henon-Henon system is synchronised in such a way that both
systems are statistically indistinguishable. Analogous results were obtained
using the parameter $d=3,5$.

Similar behaviour has been already observed on other information quantifiers
such as the transfer entropy [19, 20, 21]. It is important to note that for
the function $g(x)=log(x)$ (that is, the JSD) the values between the maximum
and zero is $\sim 0.3$ almost an order of magnitude higher than for the other
functions $g(x)$.

Figure 1:  Henon-Henon coupled system. Boxplot $\gamma$-divergence as a
function of the coupling parameter $\epsilon$. $\mathcal{D}_{\gamma}$ was
calculated with function $g(x)$ A) $e^{x}$, B) $log(x)$, C) $\sqrt[]{x}$, D)
$sinh(x)$. The ordinal patterns parameter used were $d=4$ and $\tau=1$.

###  4.2 Dynamical changes detection

In this example we use the $\gamma$-divergence to detect changes in the
dynamics of a signal. For this, we generated a group of mixed signals:
stochastic-chaotic and chaotic-chaotic. For the chaotic signal we use the
logistic map [22] given by the equation
$x_{n+1}=r~{}x_{n}\left(1-x_{n}\right)$ for $r=4$ and the cubic map [23]
$x_{n+1}=A~{}x_{n}\left(1-x_{n}^{2}\right)$ for $A=3$. For stochastic signal
we use a white noise. We generate $N=100$ mixed signals
$S_{s-c}^{i}=S_{s}+S_{c}=2000+2000=4000$ and
$S_{c_{1}-c_{2}}^{i}=S_{c_{1}}+S_{c_{2}}=2000+2000=4000$ with $i=1,...,100$
(Figure 2 top). For each signal we used a pointer method to study the
divergence. This method consists of a sliding pointer $p$ which move point to
point over the signal in the range $2d!\leq i\leq L-2d!$ splitting the signals
in two subsequences $S_{l}[x(1),\ldots,x(p)]$ and $S_{r}[x(p+1),\ldots,x(L)]$.
In this case we apply the weighted divergence define in eq.7 with weight
$w_{l}=L_{S_{l}}/L$ and $w_{r}=L_{S_{r}}/L$. We compute the
$\mathcal{D}_{\gamma}^{w_{l},w_{r}}(S^{l}||S^{r})$ for all $i$-signal and
finally measure the mean value $\mu=<\mathcal{D}_{\gamma}>_{i}$ and the
standard deviation $\sigma$. In figure 2 are shown the results for the two
studied cases: white noise-logistic map signal (A) and cubic-logistic maps
signal (B). The figure shows the $\mu$ (straight line) and $\sigma$ (shadow
bar) for the four above indicated $g(x)$ functions. I both cases the maximum
value ($\mathcal{D}_{\gamma}^{max}$) of the divergence is reached in the
transition between the two signals. Is important to remark that for every
function $g(x)$ the method detects the change in the dynamics but the highest
values of $\mathcal{D}_{\gamma}^{max}$ corresponds to $g(x)=log(x)$.

Figure 2: Dynamical changes detection in mixed signal using
$\gamma$-divergence for different function $g(x)$. A) Mixed signal white noise
- logistic map, B) Mixed signal cubic map - logistic map. The ordinal patterns
parameter used were $d=4$ and $\tau=1$.

###  4.3 Detection of different states of sleep

In this example we use the proposed scheme to identify the transition between
sleep states from an electroencephalogram (EEG) record.

The problems associated with sleep have serious repercussions on people’s
health. That is why their study is of great clinical relevance. [24, 25, 26].

Usually two primary sleep stages are identified: rapid eye movement (REM)
sleep and non-REM sleep (NREM). REM is defined as an active sleeping period
showing an intense brain activity. Brain waves remain fast and desynchronised
similar to those in the waking state. In NREM sleep stage the physiological
activity decrease, the brain waves get slower and have greater amplitude,
breathing, and heart rate slows down and blood pressure drops. The NREM phase
is composed by three stages: N1, N2, and N3. The N1 stage is characterised by
perceived drowsiness or transition from being awake to falling asleep observed
by slowing down the brain waves and muscle activity. Stage N2 is a period of
light sleep during which eye movement stops. Brain waves become slower (Theta
waves (4-7 Hz)) with occasional bursts of rapid waves (12-14 Hz) called sleep
spindles, coupled with spontaneous periods of muscle tone mixed. Lastly, stage
N3 is characterised by the presence of Delta (0.5-4 Hz) slow waves,
interspersed with smaller, faster waves [27]. The N3 stage is a deep sleep
stage, without eye movement, and a decrease of muscle activity, resembling a
coma state. Usually, sleepers pass through these four stages (REM, N1, N2, and
N3) cyclically. A complete sleep cycle takes an average of 90 to 110 minutes,
with each one lasting between 5 to 15 minutes. These cycles must be maintained
for healthy body function in awake state [28]. Developing tools that can
detect the changes in dynamic sleep stages over EEG signal are highly
essential for studying patients with sleep disorders [29]. Here we use our
$\gamma$-divergence to detect changes in the dynamics of the EEG to detect the
transition from one state to another.

The data were taken from the Physionet database: The Sleep-EDF Database
[Expanded] [30, 31], and are freely available at [32]. The EEG were recorded
(Fpz-Cz) thorough bipolar channels and the sampling frequency was $100$ Hz.
Initially, we extracted segments from the original signal belonging to the
five different sleep states: Awake, REM, N1, N2, N3 111We used the notes
provided by the database. Each segments had $6000$ points (corresponding to
$60$ sec of recording) and were joined into a single signal 222This was done
for a better visualisation of the results, because the time of each sleep
state can vary from seconds to several minutes. (Figure 3A.) The signal was
preprocessed with a band-pass filter between $0.5-60$ Hz. The ordinal pattern
were computed using the parameter $d=4$ and $\tau=1$.

Figure 3B shows that our method could detect the transition between different
sleep states for all the used functions $g(x)$. The $\log(x)$ function shows
the highest values and the best differentiation between stages. The transition
between Awake-REM and N2-N3 are more remarkable than in REM-N1 and N1-N2.
There is a significant differentiation between N2-N3 because of N3 is the
deepest sleep stage. In this stage the body becomes more insensitive to
outside stimuli. In this state, the EEG signal is mostly composed by slow
waves (Delta and Theta) causing the brain dynamics to be sharply different
from the other states. Lower $\gamma$-divergences values were found between
N1-N2 states showing that both states share similar characteristics in their
dynamics. Particularly, N1 is characterised by drowsiness slowing down the
brain waves and muscle activity. N1-N2) are very similar stages being the main
difference that in N1 occasional bursts of rapid waves (12-14 Hz), called
sleep spindles appear.

Similar result can be found between REM and N1. REM and N1 also present lower
values of the gamma divergence. This is something to be expected considering
that the EEG during the REM stage contains frequencies present in the “awake”
state and in the lighter stages of sleep N1 [33]. Despite the similarities of
the REM and N1, there are still enough differences between them such that
statistically different values are obtained. The presence of $11-16$ Hz
activity (sleep spindles) in N1, and more abundant alpha activity ($8-13$Hz)
in REM sleep means that these two stages present activity at an overlapping
frequency range, which explains the proximity of the divergence values
obtained. Difficulty in detecting N1 and REM sleep has also been found using
other measures [34, mateos2021using].

Figure 3:  Application of the $\gamma$-divergence over a sleep EEG signal
using a running windows. A) The EEG signal is composed by $5$ sub-signal
belonging to different sleep stages (Awake, REM, N1, N2, N3). Each state has
$6000$ point corresponding to $60$ sec recording (dashes horizontal lines). B)
Results corresponding to the running windows method (see text for details).
The signal was quantified with the permutation vectors with parameter $d=4$
and $\tau=1$. For all functions, the maximum values
$\mathcal{D}_{\gamma_{max}}$ were reached in the exact point where a
transition between sleep states exists.

###  4.4 Distance between textured images

As a last application we use the $\gamma$-divergence to measure the distance
between 2D textured images. The images were taken from the Normalized Brodatz
Texture (NBT) album avalible in
https://multibandtexture.recherche.usherbrooke.ca. This normalised database
has 112 texture images an improvement regarding the original Brodatz texture
database since that grayscale background effects have been removed [35].

Because of this, it is impossible to discriminate between textures from this
normalised database using only ﬁrst order statistics. Six different samples of
the NBT database are illustrated in Figure 4A. The images of the NBT album
have dimensions of $640\times 640$ pixels and 8 bits/pixel, which provides 256
grayscale levels.

For the 2D quantization, the ordinal parameters used were $d_{x}=d_{y}=2$ and
$\tau_{x}=\tau_{y}=1$. Figure 4B shows the divergence matrices calculated for
the different functions $g(x)$. We can see all the divergence values allow us
to distinguish different textures from each other. The divergence is higher
when the texture patterns are more different, for example between D15 and D47,
and lower when the patterns have a similar geometry for example between D1 and
D71. The four divergences show similar relative values between the images
however, the values for $g(x)=log(x)$ have an order of magnitude higher than
the other functions.

Figure 4:  A) Six samples of the NBT album with their corresponding labels are
illustrated. B) Location of the 112 texture images from the NBT album in the
CECP. $\gamma$-divergence matrix between textured images showing in A using 2D
ordinal patterns with $d_{x}=d_{y}=2$ and $\tau_{x}=\tau_{y}=1$.

##  5 Discussion

We showed some particular aspects of the JSD when thought of as a BRC. In
particular that is the only BRC that is a Csiszar divergence besides being one
of the few that has a metric character. Then we defined the weighted BRC,
which can be explored in the context of Bayesian inference. We restrict
ourselves to the application of these weighted version to the analysis of the
dynamical behaviour of time series. Afterwards, we showed that the ordinal
patterns method is a natural way to apply divergence in signal and image
analysis. The ordinal patterns allow a natural quantification for the use of
one special kind of BRC that we called $\gamma$-divergence.

We applied the $\gamma$-divergence to four examples, three for signals
analysis and one for image analysis. Among the signal analyses, we studied the
coupling between two signals generated by the Hennon-Hennon system. We
observed that divergence can quantify the coupling of the system depending on
the coupling parameter. In the second example, we use divergence to detect
changes in the dynamics of combined signals. We noted that divergence can
accurately detect the transition point between two chaotic-stochastic and
chaotic-chaotic systems. Later, we analysed EEG signal from sleep patient. We
could detect the points where the signal changes its dynamics because of the
change in the state of sleep. Showing it can be an alternative tool for
detecting different sleep states. Finally, we used this family of divergences
to study distances between textured images. We could realised divergence
values were smaller for similar textures and larger for different textures.

One interesting result is for all the examples studied, the largest divergence
values occur for the generating function $g(x)=log(x)$ which is precisely the
Jensen-Shannon divergence.

It is also important to mention that when using divergences, we had consider
the significance of the value obtained. This significance is what allows us to
say if the distance found is really a real and not just a measure of the
statistical fluctuations. In this work we used the standard deviation
calculated over $N$ realizations as a measure of significance. However, many
times we do not have more than one realization being this method useless.
Therefore, a theoretical study of the field is necessary. Some studies on this
topic have already been addressed by Grosse et al. for the Jensen-Shannon
divergence [36]. Similar studies should be carried out for this family of
$\gamma$-divergence. However, this exceeds this work but will be a addressed
in the near future.

##  6 Aknowledgements

We are grateful to Secretaria de Ciencia y Tecnica de la Universidad Nacional
de Córdoba, Argentina, for financial assistance. This work was partially
supported by Grant PIP 519/2019 from CONICET, Argentina and by the Grant PICT
2019-01750 from Angencia Nacional de Promoción de la Investigación, el
Desarrollo Tecnológico y la Inovación, Argentina.

## References

  * [1] R. A. Fisher.  Theory of statistical estimation.  In Mathematical Proceedings of the Cambridge Philosophical Society, volume 22, pages 700–725. Cambridge University Press, 1925.
  * [2] T. M. Cover and J. A. Thomas.  Elements of information theory.  John Wiley & Sons, 2012.
  * [3] M. Parry and E. Fischbach.  Probability distribution of distance in a uniform ellipsoid: Theory and applications to physics.  Journal of Mathematical Physics, 41(4):2417–2433, 2000.
  * [4] S. N. Ondimu and H. Murase.  Effect of probability-distance based markovian texture extraction on discrimination in biological imaging.  Computers and Electronics in Agriculture, 63(1):2–12, 2008.
  * [5] A. Kostin.  Probability distribution of distance between pairs of nearest stations in wireless network.  Electronics Letters, 46(18):1299–1300, 2010.
  * [6] A. P. Majtey, P. W. Lamberti, and D. P. Prato.  Jensen-shannon divergence as a measure of distinguishability between mixed quantum states.  Physical Review A, 72(5):052310, 2005.
  * [7] G. Crooks and D. Sivak.  Measures of trajectory ensemble disparity in nonequilibrium statistical dynamics.  Journal of Statistical Mechanics: Theory of Experiments, P06003, 2011.
  * [8] T. M Osán, D. G. Bussandri, and P. W. Lamberti.  Monoparametric family of metrics derived from classical jensen–shannon divergence.  Physica A: Statistical Mechanics and its Applications, 495:336–344, 2018.
  * [9] A. P. Majtey, P. W. Lamberti, and A. Plastino.  A monoparametric family of metrics for statistical mechanics.  Physica A: Statistical Mechanics and its Applications, 344(3-4):547–553, 2004.
  * [10] I. Csiszár.  Information measures: A critical survey.  In Transactions of the Seventh Prague Conference on Information Theory, Statistical Decision Functions, Random Processes, pages 73–86, 1974\.
  * [11] J. Burbea and C. R. Rao.  Entropy differential metric, distance and divergence measures in probability spaces: A unified approach.  Journal of Multivariate Analysis, 12(4):575–596, 1982.
  * [12] S. M. Ali and S. D. Silvey.  A general class of coefficients of divergence of one distribution from another.  Journal of the Royal Statistical Society: Series B (Methodological), 28(1):131–142, 1966.
  * [13] C. Bandt and B. Pompe.  Permutation entropy: a natural complexity measure for time series.  Physical Review Letters, 88(17):174102, 2002.
  * [14] Massimiliano Zanin and Felipe Olivares.  Ordinal patterns-based methodologies for distinguishing chaos from noise in discrete time series.  Communications Physics, 4(1):1–14, 2021.
  * [15] H. V. Ribeiro, L. Zunino, E. K. Lenzi, P. A. Santoro, and R. S. Mendes.  Complexity-entropy causality plane as a complexity measure for two-dimensional patterns.  2012\.
  * [16] L. Zunino and H. V. Ribeiro.  Discriminating image textures with the multiscale two-dimensional complexity-entropy causality plane.  Chaos, Solitons & Fractals, 91:679–688, 2016.
  * [17] L. Zunino, M. C. Soriano, I Fischer, O. A. Rosso, and C. R. Mirasso.  Permutation-information-theory approach to unveil delay dynamics from time-series analysis.  Physical Review E, 82(4):46212, 2010.
  * [18] L. Zunino, M. C. Soriano, and O. A. Rosso.  Distinguishing chaotic and stochastic dynamics from time series by using a multiscale symbolic approach.  Physical Review E, 86(4):46210, 2012.
  * [19] V. Paluš, M.and Komárek, Z. Hrnčíř, and K. Štěrbová.  Synchronization as adjustment of information rates: Detection from bivariate time series.  Physical Review E, 63(4):046211, 2001.
  * [20] A. Krakovská, J. Jakubík, H. Budáčová, and M. Holecyová.  Causality studied in reconstructed state space. examples of uni-directionally connected chaotic systems.  arXiv preprint arXiv:1511.00505, 2015.
  * [21] J. F. Restrepo, D. M. Mateos, and G. Schlotthauer.  Transfer entropy rate through lempel-ziv complexity.  Physical Review E, 101(5):052117, 2020.
  * [22] R. M. May.  Simple mathematical models with very complicated dynamics.  In The Theory of Chaotic Attractors, pages 85–93. Springer, 2004\.
  * [23] Z. Wan-zhen, D. Ming-zhou, and L. Jia-nan.  Symbolic description of periodic windows in the antisymmetric cubic map.  Chinese Physics Letters, 2(7):293, 1985.
  * [24] R. M. Benca, M. Okawa, M. Uchiyama, S. Ozaki, T. Nakajima, K. Shibui, and W. H. Obermeyer.  Sleep and mood disorders.  Sleep medicine reviews, 1(1):45–56, 1997.
  * [25] R. Stickgold.  Sleep-dependent memory consolidation.  Nature, 437(7063):1272–1278, 2005.
  * [26] J. J. Pilcher and A. S. Walters.  How sleep deprivation affects psychological variables related to college students’ cognitive performance.  Journal of American College Health, 46(3):121–126, 1997.
  * [27] R. Boostani, F. Karimzadeh, and M. Nami.  A comparative review on sleep stage classification methods in patients and healthy individuals.  Computer methods and programs in biomedicine, 140:77–91, 2017.
  * [28] R. D. Ogilvie.  The process of falling asleep.  Sleep Medicine Reviews, 5(3):247–270, 2001.
  * [29] G. Tononi and C. Cirelli.  Sleep function and synaptic homeostasis.  Sleep medicine reviews, 10(1):49–62, 2006.
  * [30] A. L. Goldberger, L. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C. K. Peng, and H. E. Stanley.  Physiobank, physiotoolkit, and physionet: Circulation.  Discovery, 101(23):1, 1997.
  * [31] B. Kemp, A. H. Zwinderman, B. Tuk, H. A. C. Kamphuisen, and J. J. L. Oberye.  Analysis of a sleep-dependent neuronal feedback loop: the slow-wave microcontinuity of the eeg.  IEEE Transactions on Biomedical Engineering, 47(9):1185–1194, 2000\.
  * [32] Physionet databanck.  Sleep-EDF Database Expanded.  $\backslash$url{https://physionet.org/content/sleep-edfx/1.0.0/}.
  * [33] N. Nicolaou and J. Georgiou.  The use of permutation entropy to characterize sleep electroencephalograms.  Clinical EEG and Neuroscience, 42(1):24–28, 2011.
  * [34] Q. Noirhomme, M. Boly, V. Bonhomme, P. Boveroux, C. Phillips, P. Peigneux, and Others.  Bispectral index correlate with regional cerebral blood flow during sleep.  Archives Italiennes de Biologie, 147(1/2):51–57, 2009.
  * [35] S. Abdelmounaime and H. Dong-Chen.  New brodatz-based image databases for grayscale color and multiband texture analysis.  International Scholarly Research Notices, 2013, 2013.
  * [36] I. Grosse, P. Bernaola-Galván, P. Carpena, R. Román-Roldán, J. Oliver, and H.E. Stanley.  Analysis of symbolic sequences using the jensen-shannon divergence.  Physical Review E, 65(4):041905, 2002.