# Density of states in locally ordered amorphous organic semiconductors:
emergence of the exponential tails

S.V. Novikov<EMAIL_ADDRESS>A.N. Frumkin Institute of Physical
Chemistry and Electrochemistry, Leninsky prosp. 31, Moscow 119071, Russia
National Research University Higher School of Economics, Myasnitskaya Ulitsa
20, Moscow 101000, Russia

###### Abstract

We present a simple model of the local order in amorphous organic
semiconductors which naturally produces a spatially correlated exponential
density of states (DOS). The dominant contribution to the random energy
landscape is provided by electrostatic contributions from dipoles or
quadrupoles. An assumption of the preferable parallel orientation of neighbor
quadrupoles or antiparallel orientation of dipoles directly leads to the
formation of the exponential tails of the DOS even for a moderate size of the
ordered domains. The insensitivity of the exponential tail formation to the
details of the microstructure of the material suggests that this mechanism is
rather common in amorphous organic semiconductors.

##  I Introduction

The energy density of states (DOS) is one of the most fundamental
characteristics of the amorphous semiconductors. It describes the distribution
density of the random carrier energy $E$ and defines or greatly influences all
charge transport properties of the materials. In organic semiconductors the
exact DOS $P(E)$ depends on the particular features of the material in
question and typically has a complex shape. Taking into account that for the
quasi-equilibrium charge transport a behavior of the low energy tail is
relevant, simple DOS models emphasizing the particular shape of the tail are
widely used in analytical or computer simulation analysis of the charge
transport properties.

There are numerous possible sources of energetic disorder in amorphous organic
semiconductors; citing just a few of them we could mention electrostatic
disorder, originated from randomly located and oriented permanent dipoles and
quadrupoles, [1, 2, 3] conformational disorder, [4] or mesoscopic disorder
associated with various irregularities at the boundaries between microcrystals
(i.e., ordered domains of organic material).[5, 6, 7, 8, 9] In this paper we
limit our consideration to the case where the electrostatic contribution to
the DOS is dominant. Detailed microscopic computer simulation suggests that
this is possible in some typical organic semiconductors.[4]

A popular model of the DOS in amorphous organic materials is the Gaussian
DOS.[10] Not surprisingly, the Gaussian DOS quite naturally emerges in the
models where the electrostatic contribution is dominating because the long
range nature of the unscreened electrostatic potential implies the
applicability of the Central Limit Theorem.[1, 2, 3] If we consider the
simplest case of the low carrier density and weak applied electric field, then
for the thick transport layers having long transit time the quasi-equilibrium
distribution of occupied states develops

$P_{eq}(E)\propto
P(E)\exp(-E/kT)\propto\exp\left[-\frac{\left(E-E_{eq}\right)^{2}}{2\sigma^{2}}\right],\hskip
10.0ptE_{eq}=-\frac{\sigma^{2}}{kT},$ (1)

and the functional form of the DOS at $E\simeq E_{eq}$ determines properties
of the quasi-equilibrium charge transport. For typical rms disorder
$\sigma\simeq 0.1$ eV and room temperature $\sigma/kT\approx 4-5$,
$P(E_{eq})/P(0)\simeq 10^{-6}-10^{-4}$ and relevant transport states belong to
the tail of the DOS. At the same time, tails of the distribution are typically
affected by the subtle details of the material structure and, quite
frequently, the shape of the distribution in this region notably deviates from
the Gaussian one.

Most popular electrostatic models use the basic assumption of the absolute
independent orientations of dipoles or quadrupoles, [1, 2, 3] yet it is
evident that this assumption cannot be valid at least for neighbor molecules.
Quite large and asymmetric molecules of organic semiconductors just cannot be
oriented in a true random and independent way at the short distance.
Obviously, short range correlation should affect the density of states in
amorphous materials as well as the correlation properties of the random energy
landscape and, hence, the charge transport behavior of organic materials. Such
effects were indeed observed.[11, 12]

Additional argument for the local order comes from the different and rather
unexpected side. Typically, rms electrostatic disorder $\sigma$ was estimated
using the macroscopic value of static dielectric constant $\varepsilon$. This
approximation is unreasonable for the short range distance comparable to the
intermolecular distance. Qualitatively, more proper description could be
achieved by setting $\varepsilon\approx 1$ for short distances,[13] leading to
the greater $\sigma$ value. Indeed, it was demonstrated that the more accurate
microscopic calculation provide approximately twice greater disorder, i.e.
instead of $\sigma\simeq 0.1$ eV we should rather have $\sigma\simeq 0.2$
eV.[14] Such value is unreasonably high and disagree with the experimentally
obtained $\sigma$.[15, 16] Again, the plausible explanation for the agreement
is the effect of the local order diminishing the total energetic rms disorder.

The tail of the DOS in amorphous organic materials is frequently described by
the simple exponential form[17, 18, 19, 20]

$P(E)\propto\exp(E/E_{e}),\hskip 10.0ptE\rightarrow-\infty,$ (2)

and this type of DOS leads to a drastically different transport behavior for
low temperature $kT<E_{e}$, where transport becomes dispersive with no well-
defined time-independent average carrier velocity. Several models of the
development of the exponential or near-exponential tails in amorphous
materials have been suggested and some of then are hardly applicable to
organic materials while others required the presence of charged dopants.[21,
22, 23, 24] Until now the exponential tails in amorphous organic
semiconductors are mostly attributed to the effect of charged dopants.

In this paper we suggest a new and probably quite general model of the
development of the exponential tails of the DOS in locally ordered organic
glasses with dominating electrostatic disorder without presence of charged
dopants. The model is directly applied to quadrupolar materials and may be
applied to dipolar materials with the particular type of the local ordering.

##  II Model of the local order in quadrupole glass

The quadrupole moment of a molecule is described by the symmetric traceless
tensor Q, $\sum\limits_{i}Q_{ii}=0$. Interaction energy for the charge and
point quadrupole separated by the distance $\vec{r}$ is

$E_{q}(\vec{r})=\frac{e\sum\limits_{i,j}Q_{ij}r_{i}r_{j}}{2\varepsilon
r^{5}}.$ (3)

Quadrupole tensor may be transformed to the diagonal form where any quadrupole
could be considered as a linear combination of axial and planar quadrupoles
$\textbf{Q}=\textbf{Q}_{a}+\textbf{Q}_{p}$, where
$\textbf{Q}_{a}=\textrm{diag}(-Q_{a}/2,-Q_{a}/2,Q_{a})$ and
$\textbf{Q}_{p}=\textrm{diag}(Q_{p},-Q_{p},0)$. Even for absolutely random
non-correlated distribution of quadrupoles the cases of pure axial and planar
quadrupoles are different: for planar quadrupoles the DOS is symmetric around
$E=0$ while for the axial quadrupoles the tails are asymmetric. [3] For planar
quadrupoles the inversion $Q\rightarrow-Q$ that transforms the positive
carrier energy to negative one and vice versa is equivalent to the spatial
rotation of the quadrupole. All spatial configurations of the quadrupole have
the equal weight, hence, the DOS is symmetric around $E=0$. For axial
quadrupoles the inversion is not equivalent to any possible rotation and the
DOS is not symmetric.

Naturally, we should also expect such difference for locally ordered materials
without any preferred direction and consider the cases of axial quadrupole
glass (QG) and planar QG separately.

In the simplest QG model quadrupoles occupy sites of the simple cubic lattice
having lattice scale $a$ and all lengths are measured in the units of $a$. For
the true amorphous organic material with random position of molecules $a$ is
the distance between nearest neighbors and typically $a\simeq 1$ nm. All
energies will be measured in the units of characteristic energy
$E_{0}=eQ/2\varepsilon a^{3}$, where $Q$ is the axial or planar quadrupole
moment.

Again, we consider the simplest possible model of the local order assuming
that the amorphous organic material is built from the blocks of quadrupoles
and in any particular block they have the identical orientations while the
corresponding orientations in different blocks are independent and absolutely
random. The blocks are built by the Voronoi tesselation.[25] We set $N_{s}$
randomly located seeding sites in the cubic basic sample having linear size of
$L$ sites (with proper periodic boundary conditions) and then subdivide the
basic sample into $N_{s}$ cells where sites of each cell are more close to one
of the seeding sites than to all others. Average volume of the cell is
$L^{3}/N_{s}$ sites and the typical linear size of a cell may be estimated as
$l_{a}=L/N_{s}^{1/3}$ sites. For the majority of simulations we use $L=50$ and
comparison with some particular simulations carried out for $L=100$ shows no
significant difference.

Boundaries of each cell are planes, thus providing sharp planar interfaces
between cells. In order to check an effect of the sharp interfaces on the
shape of the DOS we consider also the smoother version of the boundary by
applying the following procedure: sweeping over the lattice and replacing
every quadrupole with the average over nearest quadrupoles (with proper
rescaling to preserve the magnitude of the quadrupole moment). In the depth of
the Voronoi cell this procedure does not make any difference, while at the
boundary it provides the gradual transition from one cell orientation to
another (see Fig. 1). With the increase of the number of sweeps $m$ the
resulting distribution of quadrupoles orientations becomes smoother and
smoother.

a)

b)

Figure 1: Distribution of the axes of quadrupoles for $15\times 15\times 15$
chunk of the axial Voronoi QG with $l_{a}=5$ for $m=0$ (a) and $m=3$ (b). Note
that due to projection effect for some sites several axes seem to belong to
the single site; in fact, they are located at the different sites which are
seen overlapping in projection. Figure 2: Shape of the central region of the
DOS peak for $l_{a}=6.3$ for the axial ($\blacktriangle$) and planar
($\vartriangle$) locally ordered QG ($L=50$). Lines show the best fit for the
Gaussian function. Inset shows the dependence of $\sigma^{2}$ on $l_{a}$ for
planar ($\vartriangle$) and axial ($\blacktriangle$) QG.

a)

b)

Figure 3: Emergence of the exponential tails of the DOS in the axial (a) and
planar (b) QG for $l_{a}=2.9$ ($+$) and $l_{a}=10$ ($\triangledown$), $L=50$.
Central region of the DOS is approximately Gaussian with the same width
irrespectively of $l_{a}$. Straight lines show the best exponential fits for
the tails (plotted only for a well developed exponential asymptotics).

##  III Results and Discussion

###  III.1 Central peak of the DOS

For the absolutely random model of the QG it was already found that the DOS is
mostly Gaussian, though for the axial QG the low energy side shows deviation
from the Gaussian shape for rather moderate $E$ (decays faster than the
Gaussian one), while for the planar QG the DOS is symmetric and retains the
Gaussian shape in a more broad region.[3] We describe here the behavior of the
low and high energy sides of the DOS assuming $eQ>0$, for the opposite sign of
the carrier the sides are interchanged.

If we consider the main body of the DOS for the locally ordered material, then
the variance should be approximately unaffected by $l_{a}$. Indeed, for the
totally random QG [3]

$\sigma\propto\frac{eQ}{\varepsilon a^{3}}.$ (4)

Correlated glass could be approximately considered as the random lattice
having scale $al_{a}$ and built by quadrupoles with the quadrupole moment
$Ql_{a}^{3}$, thus $\sigma$ remains approximately constant. This consideration
should be valid for for moderate and large $l_{a}$ and is indeed supported by
the simulation data (Fig. 2, inset).

For the random glass the variance $\sigma^{2}$ is proportional to the only
quadratic invariant $\sum\limits_{i,j}Q^{2}_{ij}$, so for the same value of
$Q$ there is a simple relation between variance for axial and planar QG

$\sigma^{2}_{a}=\frac{3}{4}\sigma^{2}_{p},$ (5)

and the previous consideration hints that it should be valid for locally
ordered glasses, too. Inset of Fig. 2 shows that this is indeed so.
Universality means also that the main body of the distribution should have
approximately the same shape, and it indeed retains the Gaussian shape (Fig.
2).

###  III.2 Emergence of exponential tails

The most prominent feature of the DOS in locally ordered QG is the development
of the exponential tails for moderate degree of order (see Fig 3). Axial QGs
demonstrate the high energy tail which is rather easily transforming to the
exponential one and the low energy tail which is much more stable and retains
its Gaussian shape even for rather large $l_{a}$ (Fig. 3a). Characteristic
decay parameter of the tail grows with the increase of $l_{a}$ (Fig. 4). For
low $l_{a}$ (e.g, for $l_{a}\leq 3$) there is no indication of the development
of the exponential tails in the reachable region of $E$.

Figure 4: Characteristic decay parameters (equivalents of $E_{e}$ from Eq. (2)
for the corresponding tails) of the exponential tails (for
$E\rightarrow\pm\infty$) for planar ($\vartriangle$, here $E_{+}\approx
E_{-}$) and axial ($\blacktriangle$, only $E_{+}$ is shown) QG. Low energy
tail for the axial glass is approximately Gaussian for rather high $l_{a}$.
Values for $l_{a}\simeq 3-4$ are not very reliable due to poor statistics.

Probably, more convincing evidence for the gradual formation of the
exponential tails is the DOS plot as $\ln P(E)$ vs $E^{2}$ (Fig. 5, axial QG).
Low $E$ tail is almost perfectly Gaussian as a result of the competition of
two factors acting in the opposite directions: initially, in the random axial
QG the low $E$ tail decays faster than the Gaussian one, but the local order
leads to a more slow decay, so the resulting tail becomes more close to the
true Gaussian shape. High $E$ tail demonstrates a clear development of the
exponential tail.

Figure 5: DOS for the axial QG with $l_{a}=5$: low($\circ$) and high
($\bullet$) energy tails, correspondingly. Solid line shows the Gaussian fit
for the central peak and the broken line shows the exponential fit.

###  III.3 Spatial distribution of the deep states

It turns out that the more simple model of the local order where the
microstructure is provided by the small cubic chunks of quadrupoles with the
linear size of $l_{a}$ having the same orientation gives approximstely the
same DOS (see, for example, the comparison in Fig. 6). This analogy gives us a
possibility to carry out a simple analysis of the spacial location of the deep
or, equivalently, very high states. Such states are mostly located in the
outer regions of the small cubes (Fig. 7).

This analysis shows that the spatial distribution of deep states is very
different in the locally ordered organic semiconductors in comparison to the
totally disordered materials. In disordered semiconductors with electrostatic
DOS (dipolar or quadrupolar) there is a cluster structure in the distribution
of random energies, and in the cluster the absolute value of $E$ increases
from the outer region into the inner part of the cluster. [26, 27] In
semiconductors with the structural local order the deepest states are located
at the interfaces between ordered granules, and these states form clusters of
sites with close energies. Clusters are spreading along the interfaces (Fig.
8). Thus, charge transport occurs over 2D random manifold of interfaces with
occasional incursions into the bulk of the granules.

Figure 6: Comparison of the DOS for the axial QG and $l_{a}=7$ for the Voronoi
($\vartriangle$) and small cubes ($+$) models. Figure 7: Partial DOS for the
small cubes model of the axial QG with $l_{c}=7$ for states with various
distance to the cube surface: 0 ($+$), 1 ($\circ$), 2 ($\bullet$), and 3
($\vartriangle$) lattice scales, correspondingly. Figure 8: Distribution of
the deepest (light color) and highest (dark color) states ($|E|\geq
0.5E_{max}$, $E_{max}$ is the maximal energy for this particular sample) in
the $L=50$ cubic sample of the axial Voronoi QG for $l_{a}=17$, high value of
$l_{a}$ was chosen for the better visibility of the structure. Radius of the
ball is proportional to $|E|$. Concentration of the deepest/highest states at
the interfaces (outer shells of the Voronoi polyhedra) is clearly visible.

To test the importance of the sharp interfaces between Voronoi cells for the
emergence of the exponential tails we performed a limited test of the DOS
features for the smoothed model of the interface (Fig. 9). Spatial
orientations of the quadrupole axes for the sharp and smooth model of the
interface is shown in Fig. 1. We see that the sharp interface is not a
necessary requirement of the formation of the exponential tail; in fact, for
the smooth interface the tails become even more pronounced. Additional
analysis shows that for the smooth version of the small cube model the deepest
states again are located at the interfaces. We may conclude that the necessary
condition for the exponential tails formation is the parallel orientation of
quadrupoles at short distance. Fig. 9 provides an additional support for the
sufficiency of $L=50$ for the proper evaluation of the DOS.

Figure 9: DOS for the axial Voronoi QG ($l_{a}=5$) for $L=50$, $m=0$ ($+$),
$L=50$, $m=3$ ($\triangledown$), $L=100$, $m=0$ ($\vartriangle$), and $L=100$,
$m=3$ ($\blacktriangle$), correspondingly.

###  III.4 Locally ordered dipole glasses

Until now we considered locally ordered quadrupolar glasses, and it turned out
that the local order for moderate or large $l_{a}$ leads to the development of
the exponential tails of the DOS. Each quadrupole could be considered as a
pair of dipoles with equal and oppositely oriented dipole moments. We may
expect that the dipole glass with anti-ferroelectric local order could mimic
the QG and develop exponential tails as well. One of the simplest models
demonstrating such behavior could be built in two steps. First, we set up the
Voronoi model with dipoles having the same orientation in each cell but
orientations in different cells are totally uncorrelated and random (the exact
analog of the Voronoi QG). Next, we iteratively flip orientation of each
dipole $\vec{p}\rightleftarrows-\vec{p}$ in order to minimize the scalar
product $\vec{p}\cdot\sum\vec{p}_{i}$, where $\vec{p}_{i}$ are dipoles in the
nearest neighborhood of $\vec{p}$. We stop flipping when the stationary state
is reached. Fig. 10 shows the DOS for the resulting dipole glass.

Figure 10: DOS for the Voronoi dipole glass for $l_{a}=5$ ($\vartriangle$) and
$l_{a}=8.5$ ($+$, in both cases $L=50$). Straight lines show the best
exponential fit for the tails.

As in the case of QG glass, the central peak weakly depends on $l_{a}$, but
the rms disorder is strongly reduced in comparison to the case of totally
random dipole glass, where $\sigma\approx 2.35ep/\varepsilon a^{2}$, while for
Fig. 10 $\sigma\approx 0.8ep/\varepsilon a^{2}$. This is a good illustration
of the possible reduction of $\sigma$ due to local ordering. The DOS is
symmetric in close analogy with the planar QG.

We are not going to study the DOS in the dipole glasses in details, just to
demonstrate that the development of the exponential tails is possible, albeit
for the more strict limitations on the local structure of the glass (anti-
ferroelectric correlations).

###  III.5 Transition region between the Gaussian peak and exponential tail

A principal difference between the Gaussian and exponential DOS is the
energence of the dispersive transport regime for the latter at low temperature
$kT<E_{e}$, where $E_{e}$ is the characteristic decay parameter for the tail.
At the same time, if the exponential tail develops for too low energy, then
for the moderately thick transport layer the true dispersive regime may be not
observed due to the rather short transit time.

Let $E_{t}$ be the estimation for the transition region from the Gaussian
dependence to the exponential one. Assuming the smooth transition between the
Gaussian and exponential regions, $E_{t}$ may be estimated as

$\frac{E_{t}}{\sigma^{2}}\simeq\frac{1}{E_{e}}.$ (6)

For time-of-flight experiments in the quasi-equilibrium regime the exponential
tail becomes relevant if

$|E_{eq}|=\frac{\sigma^{2}}{kT}\simeq E_{t}\simeq\frac{\sigma^{2}}{E_{e}}.$
(7)

This relation indicates that the exponential tail becomes relevant if
$kT\simeq E_{e}$, i.e. in the vicinity of the transition to the dispersive
transport regime. Observation of the nondispersive transport regime which is
governed by the exponential DOS should be difficult. This conclusion should be
considered with great care because all estimates here are very approximate.

Assuming typical $\sigma\simeq 0.1$ eV, at the room temperature $kT_{r}\approx
0.025$ eV $\simeq\sigma/4$ where is no indication of the dispersive transport
(hence, $kT_{r}>E_{e}$) if

$E_{e}\simeq\frac{\sigma^{2}}{E_{t}}<kT_{r},\hskip 10.0ptE_{t}\gtrsim
4\sigma.$ (8)

For a wide range of $l_{a}$ the rms disorder can be estimated as $\sigma\simeq
1.5-2E_{0}$, so $E_{t}\gtrsim 6-8\hskip 2.0ptE_{0}$. Comparison with the DOS
figures suggests that for the room temperature the absence of the dispersive
regime implies $l_{a}\leq 3-4$.

This brief consideration shows that the description and analysis of
experimental time-of-flight curves is, probably, even more difficult and
ambiguous than suggested previously; [16, 18] emergence of the exponential
tails makes transport behavior much more complicated.

An important open problem is the continuation of the exponential tail to an
arbitrary deep energy or a possible emergence of another asymptotics for very
deep states. Computer simulation cannot solve this problem due to enormous
simulation time necessary for the accumulation of the reliable tail
statistics. Unfortunately, analytic approach to the estimation of the tail
asymptotics is not developed for locally ordered organic glasses.

##  IV Conclusions

In conclusion, we consider the simple model of quadrupolar glasses having
local order and show that for the linear spatial size of the ordered regions
$\gtrsim 5$ molecules well defined exponential tails are developed. There is a
principal difference between axial and planar quadrupoles: for axial QG the
DOS is asymmetric and the low energy tail is less prone to the development of
the exponential asymptotics, while for the planar QG the DOS is symmetric. For
this reason charge transport of electrons and holes in the axial QG may be
very different due to principal difference of the shape of the tails. Deep
sites, where carriers spend the most time, are organized in clusters located
at the interfaces between ordered granules; hence, in most cases charge
carrier transport occurs over the random 2D manifold of interfaces.

Comparison of various models suggests that the details of the ordering are not
important: what is important is the existence of the domains with preferable
parallel orientation of quadrupoles. For dipolar materials the tail is
developed for the anti-ferroelectric local ordering. Primary potential
candidates for the amorphous organic semiconductors with significant local
order are liquid crystalline materials [28, 29, 30, 31, 32] or macrocyclic
aromatic hydrocarbons, [33] but other materials could develop the sizable
local order as well, such as molecularly doped polymers [34] or low molecular
weight organic glasses. [12] Quite probably, the dispersive charge transport
observed in some organic materials [18] may be related to the development of
the exponential tails due to effects of the local order.

## Acknowledgements

Financial support from the Ministry of Science and Higher Education of the
Russian Federation (A.N. Frumkin Institute) and Program of Basic Research of
the National Research University Higher School of Economics is gratefully
acknowledged.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Dieckmann, Bässler, and Borsenberger [1993] A. Dieckmann, H. Bässler, and P. M. Borsenberger, J. Chem. Phys. 99, 8136 (1993).
  * Novikov and Vannikov [1994] S. Novikov and A. Vannikov, J. Exp. Theor. Phys. 79, 482 (1994).
  * Novikov [2015] S. V. Novikov, J. Chem. Phys. 143, 164510 (2015).
  * Massé  _et al._ [2017] A. Massé, P. Friederich, F. Symalla, F. Liu, V. Meded, R. Coehoorn, W. Wenzel, and P. A. Bobbert, Phys. Rev. B 95, 115204 (2017).
  * Rivnay  _et al._ [2009] J. Rivnay, L. H. Jimison, J. E. Northrup, M. F. Toney, R. Noriega, S. Lu, T. J. Marks, A. Facchetti, and A. Salleo, Nat. Mater. 8, 952 (2009).
  * Kalihari  _et al._ [2008] V. Kalihari, E. B. Tadmor, G. Haugstad, and C. D. Frisbie, Adv. Mater. 20, 4033 (2008).
  * Tang, Carter, and Cannon [2006] M. Tang, W. C. Carter, and R. M. Cannon, Phys. Rev. B 73, 024102 (2006).
  * Rivnay  _et al._ [2011] J. Rivnay, R. Noriega, R. J. Kline, A. Salleo, and M. F. Toney, Phys. Rev. B 84, 045203 (2011).
  * Jimison  _et al._ [2009] L. H. Jimison, M. F. Toney, I. McCulloch, M. Heeney, and A. Salleo, Adv. Mater. 21, 1568 (2009).
  * Bässler [1993] H. Bässler, Phys. Status Solidi B 175, 15 (1993).
  * Inoue  _et al._ [2005] M. Inoue, H. Monobe, M. Ukon, V. Petrov, T. Watanabe, A. Kumano, and Y. Schimizu, Opto-Electron. Rev. 13, 303 (2005).
  * Thurzo  _et al._ [2006] I. Thurzo, H. Méndez, C. Iacoviţa, and D. Zahn, Synth. Met. 156, 1108 (2006).
  * Novikov [2017] S. V. Novikov, Russ. J. Electrochem. 53, 275 (2017).
  * Madigan and Bulovic [2006] C. Madigan and V. Bulovic, Phys. Rev. Lett. 97, 216402 (2006).
  * Borsenberger and Weiss [1998] P. M. Borsenberger and D. S. Weiss,  _Organic Photoreceptors for Xerography_ (CRC Press, Boca Raton, 1998).
  * Schein and Tyutnev [2008] L. Schein and A. Tyutnev, J. Phys. Chem. C 112, 7295 (2008).
  * Vissenberg and Matters [1998] M. C. J. M. Vissenberg and M. Matters, Phys. Rev. B 57, 12964 (1998).
  * Dunlap  _et al._ [2010] D. H. Dunlap, L. B. Schein, A. Tyutnev, V. Saenko, E. D. Pozhidaev, P. E. Parris, and D. S. Weiss, J. Phys. Chem. C 114, 9076 (2010).
  * Abdalla, Zuo, and Kemerink [2017] H. Abdalla, G. Zuo, and M. Kemerink, Phys. Rev. B 96, 241202(R) (2017).
  * Xiao  _et al._ [2020] B. Xiao, P. Calado, R. C. I. MacKenzie, T. Kirchartz, J. Yan, and J. Nelson, Phys. Rev. Applied 14, 024034 (2020).
  * Silver, Pautmeier, and Bässler [1989] M. Silver, L. Pautmeier, and H. Bässler, Solid State Commun. 72, 177 (1989).
  * Arkhipov  _et al._ [2005] V. I. Arkhipov, P. Heremans, E. V. Emelianova, and H. Bässler, Phys. Rev. B 71, 045214 (2005).
  * Economou  _et al._ [1985] E. N. Economou, C. M. Soukoulis, M. H. Cohen, and A. D. Zdetsis, Phys. Rev. B 31, 6172 (1985).
  * Bacalis, Economou, and Cohen [1988] N. Bacalis, E. N. Economou, and M. H. Cohen, Phys. Rev. B 37, 2714 (1988).
  * Møller [1994] J. Møller,  _Lectures on random Voronoi tessellations_ (Springer, New York, 1994).
  * Novikov and Vannikov [1995] S. V. Novikov and A. V. Vannikov, J. Phys. Chem. 99, 14573 (1995).
  * Novikov and Van der Auweraer [2009] S. V. Novikov and M. Van der Auweraer, Phys. Rev. E 79, 041139 (2009).
  * Garnier  _et al._ [1998] F. Garnier, R. Hajlaoui, A. El Kassmi, G. Horowitz, L. Laigre, W. Porzio, M. Armanini, and F. Provasoli, Chem. Mater. 10, 3334 (1998).
  * McCulloch  _et al._ [2006] I. McCulloch, M. Heeney, C. Bailey, K. Genevicius, I. MacDonald, M. Shkunov, D. Sparrowe, S. Tierney, R. Wagner, W. Zhang, M. L. Chabinyc, R. J. Kline, M. D. McGehee, and M. F. Toney, Nat. Mater. 5, 328 (2006).
  * Bushby, Kelly, and O’Neill [2012] R. J. Bushby, S. M. Kelly, and M. O’Neill, eds.,  _Liquid Crystalline Semiconductors: Materials, Properties and Applications_ (Springer, Dordrecht, The Netherlands, 2012).
  * Eichhorn [2000] H. Eichhorn, J. Porphyrins Phthalocyanines 4, 88 (2000).
  * Woon  _et al._ [2006] K. L. Woon, M. P. Aldred, P. Vlachos, G. H. Mehl, T. Stirner, S. M. Kelly, and M. O’Neill, Chem. Mater. 18, 2311 (2006).
  * Nakanishi  _et al._ [2011] W. Nakanishi, T. Yoshioka, H. Taka, J. Y. Xue, H. Kita, and H. Isobe, Angew. Chem. 123, 5435 (2011).
  * Parashchuk  _et al._ [2013] O. D. Parashchuk, S. Grigorian, E. E. Levin, V. V. Bruevich, K. Bukunov, I. V. Golovnin, T. Dittrich, K. A. Dembo, V. V. Volkov, and D. Y. Paraschuk, J. Phys. Chem. Lett. 4, 1298 (2013).


# Micro-scale Graded Mechanical Metamaterials Exhibiting Versatile Poisson’s
Ratio

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  L. Mizzi  Department of Engineering
Sciences and Methods, University of Modena and Reggio Emilia, Reggio Emilia,
Italy  J. A. Iglesias Martínez  Institut FEMTO-ST, CNRS, Université Bourgogne
Franche-Comté, Besançon 25030, France  A. Spaggiari  Department of Engineering
Sciences and Methods, University of Modena and Reggio Emilia, Reggio Emilia,
Italy  G. Ulliac  Institut FEMTO-ST, CNRS, Université Bourgogne Franche-Comté,
Besançon 25030, France  R. Gatt  Metamaterials Unit, Faculty of Science,
University of Malta, Msida MSD 2080, Malta  Centre for Molecular Medicine and
Biobanking, University of Malta, Msida, MSD 2080, Malta  J. N. Grima
Metamaterials Unit, Faculty of Science, University of Malta, Msida MSD 2080,
Malta  Department of Chemistry, Faculty of Science, University of Malta, Msida
MSD 2080, Malta  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 control Poisson’s ratio of functional materials has been one of
the main objectives of researchers attempting to develop structures efficient
from the perspective of protective, biomedical and soundproofing devices. This
task becomes even more challenging at small scales, such as the microscale,
where the possibility to control mechanical properties of functional materials
is very significant, like in the case of flexible electronics. In this work,
we propose novel microscopic 2D and 3D functionally-graded mechanical
metamaterials capable of exhibiting a broad range of Poisson’s ratio depending
on their composition. More specifically, we show that upon adjusting the
number of structural elements corresponding to one type of the substructure at
the expense of another, it is possible to change the resultant Poisson’s ratio
of the entire system from highly positive to highly negative values as well as
to achieve arbitrary intermediate values. Finally, in addition to static
properties, we also analyze dynamic properties of the considered structures.
Namely, we show that the variation in the composition of considered mechanical
metamaterials affects the velocity of a wave propagating through the system.
This, in turn, could be essential in the case of applications utilizing
localized wave attenuation or sensors.

##  1 Introduction

In recent years, functional materials have been the subject of intense studies
due to their ability to exhibit versatile mechanical as well as other types of
behavior. Particularly important in this regard is the specific class of
functional materials called mechanical metamaterials Florijn et al. (2014);
Neville et al. (2016); Coulais et al. (2016, 2018); Wenz et al. (2021);
Christensen et al. (2015); Kadic et al. (2019); Chen et al. (2020); Frenzel et
al. (2017); Jiang et al. (2020); Yao et al. (2022); Wang et al. (2014); Gao et
al. (2019); Placidi et al. (2014); Wang et al. (2020); Zheng et al. (2022); Wu
et al. (2020); Xu et al. (2015); Ning et al. (2018); Zhang et al. (2021a);
Hamzehei et al. (2020); Sun et al. (2020) which stems from their continuously
growing applicability in numerous industries. In general, mechanical
metamaterials are rationally designed structures that can exhibit unusual
mechanical properties such as auxetic behavior (negative Poisson’s ratio)
Wojciechowski (1989); Lakes (1987); Evans and Alderson (2000); Babaee et al.
(2013); Novak et al. (2021a); Lim (2019); Jiang et al. (2022); Grima and Evans
(2000); Grima et al. (2011); Albertini et al. (2021); Novak et al. (2021b);
Ruzzene et al. (2002); Hou et al. (2014), negative stiffness Hewage et al.
(2016); Tan et al. (2022); Chen et al. (2021a); Bertoldi et al. (2017) as well
as negative compressibility Nicolaou and Motter (2012); Baughman et al.
(1998); Dudek et al. (2016); Lakes and Wojciechowski (2008). As demonstrated
over the last thirty years, these counterintuitive mechanical properties can
be utilized in a plethora of applications ranging from protective Miniaci et
al. (2016); Imbalzano et al. (2018) to biomedical Kolken et al. (2018a); van
Manen et al. (2021); Wang et al. (2016); Xue et al. (2021) and soundproofing
devices Quadrelli et al. (2021); Guenneau et al. (2007); Fleury et al. (2014,
2016). Recently, it has been also possible to observe a rapid increase in the
commercial appeal of mechanical metamaterials capable of exhibiting such
properties in the case of easily accessible products that are not necessarily
devoted to high-technology industries. A prime example of such an approach are
auxetic mechanical metamaterials used in the design of sports equipment Duncan
et al. (2018).

Despite numerous advantages of standard mechanical metamaterials, it is
necessary to note that they often share several limitations. Mechanical
metamaterials are typically designed in a way where the entire structure is
composed of one type of geometry corresponding to a specific mechanical
characteristic Grima and Evans (2000). As a result, mechanical properties of
such structures cannot normally be significantly modified which in some cases
reduces the applicability of the system. In order to overcome this limitation,
a number of new approaches were proposed that allow to significantly change
mechanical properties of the system even after it is manufactured Rafsanjani
et al. (2018). A prime example of such approach corresponds to the use of
active composite metamaterials that can have their behavior controlled via the
magnitude of the external stimulus for example in the form of the change in
the temperature or the electric and magnetic field Galea et al. (2022); Korpas
et al. (2021); Xin et al. (2020); Dudek et al. (2023). Another very
interesting example corresponds to hierarchical mechanical metamaterials Lakes
(1993); Gatt et al. (2015); Cho et al. (2014); Tang et al. (2015); An et al.
(2020); Oftadeh et al. (2014); Dudek et al. (2022b); Zhang et al. (2021b);
Dudek et al. (2022a), i.e. systems composed of substructures that can deform
irrespective from the rest of the structure. As reported in recent years,
hierarchical mechanical metamaterials can exhibit very versatile mechanical
behavior that can be controlled throughout the deformation process An et al.
(2020); Dudek et al. (2022b). Nevertheless, despite their undeniable appeal,
the manufacturing process of such active mechanical metamaterials is often
challenging and expensive relative to the standard mechanical metamaterials.
In addition, to fully showcase their tunable properties, the aforementioned
active mechanical metamaterials often require access to equipment that is not
accessible outside research facilities Galea et al. (2022). Hence, it is
relatively difficult to make applications related to these structures to be
easily available to the general public. Nevertheless, it is also possible to
observe a broad range of the mechanical behavior in the case of standard
mechanical metamaterials Mizzi et al. (2018). One of the possible ways of
achieving it corresponds to the use of several different substructures in the
design of the system Liu et al. (2020); Guo et al. (2021a); Kolken et al.
(2018b); Mirzaali et al. (2020); Dudek et al. (2019).

In the literature, there are relatively few examples of graded mechanical
metamaterials composed of several substructures Mirzaali et al. (2020); Guo et
al. (2021b); Han and Wei (2022); Niknam and Akbarzadeh (2020); Shao et al.
(2021); Pagliocca et al. (2022); Yao et al. (2021). However, despite this
direction of studies being still at its infancy, there are already several
interesting papers related to this concept. More specifically, in several
relatively recent studies Liu et al. (2020); Guo et al. (2021a); Kolken et al.
(2018b); Boldrin et al. (2016), it has been shown that the resultant
mechanical property of a graded structure (e.g. Poisson’s ratio) can be
different from the properties of its constituent substructures and typically
depends on the interplay between them. In addition, it has also been
demonstrated that multi-substructure mechanical metamaterials can exhibit
shape morphing Mirzaali et al. (2018). Nevertheless, some of the studies
related to this concept correspond to graded metamaterials where the resultant
structures are relatively unstable and can lose their ability to exhibit the
desired mechanical properties if one deforms the structure in an arbitrary
direction (especially in the case of graded composites). This in turn
significantly reduces their applicability. In addition, it should be noted
that most of the reported multi-substructure systems are only two-dimensional
which does not allow to replicate their behavior in 3D. Finally, it should be
emphasized that multi-substructure mechanical metamaterials known in the
literature are typically constructed at the macroscale. On the other hand, in
the case of numerous applications including amongst others flexible
electronics and biomedical devices, it is essential to use mechanical
metamaterials capable of exhibiting a broad range of mechanical properties at
much lower scales including the micro-scale.

In this work, we propose a novel 2D and 3D functionally-graded mechanical
metamaterial composed of multiple substructures that can exhibit a broad range
of the Poisson’s ratio. Due to the specific selection of substructures as well
as their connectivity, the deformation of the analyzed system depends
primarily on one parameter. This leads to a high stability of the structure
and allows to ensure that its deformation will not be significantly affected
by the choice of the direction in which it is being deformed. We also show
that the considered system can be constructed at the micro-scale in order to
prove its suitability for applications where the ability to exhibit specific
mechanical properties at small scales is essential.

##  2 Analyzed structures

Figure 1: a) The graphical representation of the cross-section of an example
of the considered system. The highlighted part of the structure corresponds to
the unit-cell of the system. b) SEM (Scanning Electron Microscopy) images of
experimental prototypes of the three examples (structure 1: $N_{\rm{S,x}}$ =
1, $N_{\rm{T,x}}$ = 10, structure 2: $N_{\rm{S,x}}$ = 2, $N_{\rm{T,x}}$ = 6,
structure 3: $N_{\rm{S,x}}$ = 5, $N_{\rm{T,x}}$ = 1) of quasi-2D structures
analyzed in this work. The overall $x$-dimension for the three considered
structures is equal to 1.09 mm, 0.95 mm and 1.13 mm respectively. For all
structures, their $y$ and $z$ dimensions are the same and are equal to 0.47 mm
and 0.4 mm.

In this work, we propose a novel mechanical metamaterial composed of two
different substructures corresponding to very different values of the
Poisson’s ratio. Thanks to the implementation of multiple substructures
associated with a very different mechanical profile, we aim to show that
depending on the ratio of a number of structural elements corresponding to
each of the substructures, it is possible to significantly adjust the
resultant properties of the graded metamaterial. The two mechanical
metamaterials selected in this work as substructures correspond to the
rotating squares system Grima and Evans (2000) and the rotating triangles
system with four triangles forming its rectangular unit-cell Milton (2013);
Dudek et al. (2016) (in the literature there are also other types of
mechanical metamaterials based on rotating triangles Grima et al. (2011)). The
first of these 2D substructures is known to exhibit auxetic behavior
corresponding to the Poisson’s ratio equal to -1 irrespective of the loading
direction in its idealized form. On the other hand, as shown in the literature
Dudek et al. (2016), the later system based on isosceles triangle motifs
always exhibits a positive Poisson’s ratio along the principal axes. This
means that in general, one should expect that the mechanical metamaterial
composed primarily of the rotating square-based elements should exhibit the
negative Poisson’s ratio while the graded metamaterial consisting mainly of
the rotating triangles will exhibit much larger values of the Poisson’s that
may result in the positive values of this parameter.

As shown in Fig. 1(a), the quasi-2D version of the considered mechanical
metamaterial is defined in a way where two identical auxetic substructures are
separated by the positive Poisson’s ratio substructure corresponding to the
rotating triangles. In order to quantitatively describe the amount of each
substructure within the system, we introduce variables $N_{\rm{S,x}}$ and
$N_{\rm{T,x}}$, where $N_{\rm{S,x}}$ stands for the number of square-based
unit-cells ($2\times 2$ squares) on one side of the system (measured along the
$x$-axis). Similarly, $N_{\rm{T,x}}$ corresponds to the number of triangle-
based unit-cells (2 $\times$ 2 triangles) measured along the $x$-direction.
Thus, in the case of the example presented in Fig. 1(a), one can write down
the following: $N_{\rm{S,x}}$ = 2 and $N_{\rm{T,x}}$ = 6.

In terms of the geometric parameters, for all types of structures investigated
in this study, the rotating squares used in the auxetic substructure have the
side length of $a$ = 40 $\mu$m while the isosceles triangles have dimensions
$b$ = 30 $\mu$m and $c$ = 56 $\mu$m (see Fig. 1(a)). As shown in Fig. 1(a),
the structural elements are connected to each other at vertices by means of
hinges having the non-zero thickness that varies depending on the position
within the graded metamaterial. It is also important to note that since both
types of substructures are ‘unimode’ systems, i.e. their deformation can be
described solely by means of one parameter corresponding to the angle between
structural elements, the resultant structures analyzed in this study are also
‘unimode’ systems. Thus, the deformation of the entire graded structure like
the one shown in Fig. 1(a) can be described through the change in only one
parameter. In this study, this parameter was selected as the angle between the
adjacent rotating square elements and is denoted as $\theta$. Using this
variable, one can derive analytical expressions corresponding to the linear
dimensions of the unit-cell of the entire system. This unit-cell is
highlighted by means of the yellow color in Fig. 1(a). Furthermore, upon
assuming an idealized scenario where respective squares and triangles are
connected to each other at vertices in a pin-jointed manner (i.e. freely
rotating structural elements with no flexure of the material), the linear
dimensions of such unit-cell can be defined as follows:

$L_{x}=4N_{\rm{S,x}}a\left[\sin\left(\frac{\theta}{2}\right)+\cos\left(\frac{\theta}{2}\right)\right]+2N_{\rm{T,x}}\left[c\cos\left(\frac{\gamma}{2}\right)+b\cos\left(\frac{\theta_{1}}{2}\right)\right]$
(1)

$L_{y}=2a\left[\sin\left(\frac{\theta}{2}\right)+\cos\left(\frac{\theta}{2}\right)\right]$
(2)

where,
$\gamma=\arccos\left(\left(2c^{2}-L_{y}^{2}\right)/\left(2c^{2}\right)\right)$,
$\theta_{1}$ = $\pi$ \- $\gamma$ \+ $\beta$ and $\beta$ = $\arccos$
$\left(\left(2c^{2}-b^{2}\right)/\left(2c^{2}\right)\right)$. In addition, it
is possible to determine the maximum range of the $\theta$ angle corresponding
to the analyzed geometry (i.e. the range of angles associated with
configurations that can be assumed by the system without the structural
elements overlapping):

$\theta_{min}=0^{\circ}\,\,\mathrm{or}\,\,\theta_{min}=\arcsin\left[\left(\frac{2}{ac}\right)^{2}s(s-c)^{2}(s-b)-1\right]$
(3)

$\theta_{max}=\arcsin\left(\left(\frac{c}{a}\right)^{2}-1\right)\,\,\rm{or}\,\,\theta_{max}=90^{\circ}$
(4)

where, $s=(2c+b)/2$. At this point, it is important to note that the above
analytical model is only applicable if geometric parameters allow the system
to assume the geometry with structural elements connected as shown in Fig.
1(a). It is also important to mention that all of the quasi-2D structures
considered in this work have a non-zero thickness equal to 400 $\mu$m that is
denoted by means of the $l$ variable in Fig. 1(b). In addition, it should be
emphasized that, in the latter part of this work, we also investigate the
behavior of three-dimensional structures inspired by quasi-2D structures
presented in Fig. 1(b).

##  3 Materials and Methods

###  3.1 Fabrication

In order to produce experimental micro-scale prototypes investigated in this
work, the commercial 3D printer (Photonic Professional GT+, Nanoscribe GmbH)
operating based on the two-photon lithography method was used. The negative
tone IP-S photoresin (Nanoscribe GmbH) was used to manufacture the samples
that enabled us to print structural elements with connections not smaller than
2 $\mu$m. A drop of IP-S resin was deposited on an ITO-coated soda-lime glass
substrate (dimensions: 25 × 25 × 0.7 mm3) and then photopolymerized through a
25X-objective with a femtosecond laser operating at $\lambda$ = 780 nm. The
standard printing parameters are slicing and hatching distances equal to 1
$\mu$m and 0.5 $\mu$m respectively, a laser power of 90% and a galvanometric
scanning speed of 100 mm/s. During the printing process, respective unit-cells
were polymerized individually while ensuring that there is a small overlap
between them. After the printing process, the sample was developed for the
duration of 25 min in the PGMEA solution (Propylene Glycol Methyl Ether
acetate). This procedure was conducted to remove the unexposed photoresist.
Finally, the sample was rinsed for 3 min in isopropyl alcohol (IPA) in order
to clean it from the aforementioned developer solution. At this point, it is
important to note that changing values of the above parameters such as the
laser power may significantly affect the quality of the printed samples. This
is particularly important in terms of the hinges connecting the adjacent
structural elements where the change in one of these parameters may
significantly affect their thickness. Hence, the performance of the structure
could be modified relative to computer simulations where the desired thickness
of hinges would be used.

###  3.2 Experiment

In this work, to test the mechanical properties of micro-scale prototypes
corresponding to the considered graded laminate-like mechanical metamaterials,
the samples were compressed along the $y$-direction (see Fig. 2) at the
constant speed of 1 $\mu$m s-1. To this aim, an external indenter in the form
of a polished screw was used that had a circular cross-section having a
diameter of 3 mm. It is important to note that this dimension is significantly
larger than the longest dimension of the quasi-2D structures (1.09 mm, 0.95
mm, and 1.13 mm along the x-axis for structure 1, structure 2, and structure 3
respectively). Hence, all elements constituting the top part of a given system
were compressed simultaneously during the deformation. In addition to quasi-2D
samples, as explained in the further part of this work, a 3D version of the
considered structures was also investigated. In this case, 3D samples were
subjected to the same compressive conditions as was the case for quasi-2D
prototypes. Namely, through the use of a large indenter (3 mm in diameter)
that covers the entire surface of a sample, 3D structures were compressed
along the y-axis at the constant rate of 1 $\mu$m s-1. At this point, it is
important to emphasize that throughout the compression process of both
quasi-2D and 3D prototypes, the bottom part of all of the analyzed samples was
constrained from any movement since samples were not removed from the
substrate after the printing process.

To analyze Poisson’s ratio of finite quasi-2D structures presented in Fig.
1(b), only a specific part of the sample was taken into consideration, i.e.
the highlighted unit-cell in Fig. 1(a). For these structures, Poisson’s ratio
was calculated based on the average outline of the considered cell throughout
its deformation process. More specifically, positions of the outmost vertices
constituting the initial undeformed unit-cell were tracked throughout the
deformation process to determine the average vertical and horizontal bounds of
the rectangular region associated with the unit-cell. The variation in the
shape of such rectangular region was later used to determine the strains as
well as the corresponding Poisson’s ratio. A detailed explanation of the
selection of the specific points within the unit-cell used for such analysis
is provided in the Supplementary Information. It should be also emphasized
that the selection of the specific unit-cell presented in Fig. 1(b) was not
coincidental. This stems from the fact that this specific unit-cell was
relatively far away from the constrained bottom of the system so that its
behavior was not be significantly affected by it. At the same time, one could
have not selected the topmost unit-cell to conduct the quantitative analysis
since such unit-cell was expected to exhibit significant edge effects and
would have been in direct contact with the indenter. Of course, this approach
does not resolve the issue of free edges of the considered unit-cell on the
left and right hand side of the investigated region. To avoid it, it would be
normally necessary to consider the structure composed of several unit-cells in
each row of the structure. However, given the size and complexity of the
analysed unit-cell, this would be extremely expensive and very difficult to
achieve during the manufacturing process at the micro-scale. Because of this,
in order to provide information about mechanical properties of the system
unaffected by the edge effects, in the later part of this work, we also
provide results corresponding to unit-cells with periodic boundary conditions.

In the case of 3D structures, Poisson’s ratio was calculated based on the
outline of the entire structure rather than the specific unit-cell. This stems
from the fact that the considered 3D structures have a much more complex shape
than their quasi-2D counterparts which makes it very difficult to track the
deformation of a specific part of the system. More specifically, since the
analyzed samples were fabricated at the micro-scale, cameras must operate at a
very small working distance and can capture a focused image of only those
elements that are located approximately at the same distance away from its
lens. In the case of considered 3D structures, specific parts of a sample are
located at very different distances from the lens. Thus, the best approach to
estimate Poisson’s ratio is to focus the camera image on the outmost edges of
the entire system while neglecting the behavior within the structure. As a
result, the edges of 3D samples at the micro-scale can be easily distinguished
from the background. It should be also noted that to record the deformation
process, the optical camera corresponding to the 10$\times$ (3D prototypes)
and 20$\times$ (quasi-2D structures) magnification factor was used.

###  3.3 FEM simulations

In this work, to verify experimental results and conduct other types of
simulations, the Finite Element Method (FEM) was used. First of all, the
behavior of the quasi-2D models matching exactly experimental prototypes (see
Fig. 1(b)) were simulated through the use of the commercial software COMSOL
Multiphysics. To this aim, the boundary conditions were defined in a way so
that the resultant behavior of the model would be as similar to the experiment
as possible. More specifically, the bottommost parts of the considered models
were fixed in space while the other elements of the system remained
unconstrained. Then, to induce the deformation process, a rigid, undeformable
plate was used to compress the structure. Contact elements were applied to the
surface edge nodes of the system as well as to internal surfaces, where
necessary, in order to replicate as accurately as possible the experimental
compression test. Furthermore, in the conducted simulations, the use of
geometric nonlinearities was assumed. It was also assumed that the material
corresponding to the investigated samples was isotropic and its linear
mechanical properties, in the form of the Young’s modulus and Poisson’s ratio,
were defined as 4 GPa and 0.4 respectively, i.e. mechanical properties
matching the polymerized resin used to produce the experimental prototypes. It
should be also mentioned that the element type used in the conducted
simulations was CPS4R four-node plane stress element. In addition, the average
dimension of an element was equal to 3 $\mu$m with the total number of
elements per unit-cell being equal to 34875, 60396, and 71676 respectively for
the three considered types of structures. The solver used in the conducted
simulations utilized the Quasi-Newton Method (static Implicit analysis). It
should be also noted that in the case of the anlysis of the dynamic properties
of the system, approximately the same average size of elements was used. Most
importantly, the size of elements at the edge of the unit-cell had a size of
around 0.3 $\mu$m.

In addition to FEM simulations related to finite systems, the behavior of
individual unit-cells was also investigated through FEM computer simulations
to compare it with the analytical model derived in this work and gain a better
understanding of how cells within much larger systems than those considered in
this study could behave. To this aim, the same material properties as those
defined for the finite models were implemented. However, in this case, instead
of constraining the bottom part of the system, periodic boundary conditions
were implemented following the method described in Mizzi et al. (2021) (see
Supplementary Information). Furthermore, to obtain realistic results, the
geometric nonlinearity was used in our simulations. Finally, it is worth
mentioning that for all of the FEM simulations conducted in this study, the
convergence test was conducted to determine the appropriate mesh size.

###  3.4 Parameters

In addition to the already defined parameters, it is important to mention that
the thicknesses of hinges defined in Fig. 1(a) corresponded to initial values
shown in Table 1. For all structures considered in this work, the initial
value of the angle $\theta$ was set to be equal to 45∘.

| $d_{\rm{1}}$ | $d_{\rm{2}}$ | $d_{\rm{3}}$ | $d_{\rm{4}}$ | $d_{\rm{5}}$ | $d_{\rm{6}}$ | $d_{\rm{7}}$ | $d_{\rm{8}}$ | $d_{\rm{9}}$
---|---|---|---|---|---|---|---|---|---
$\rm{S_{1}}$ | 4.2 | 4.2 | 6.0 | 6.0 | 4.0 | 3.2 | 3.8 | 8.0 | 12.0
$\rm{S_{2}}$ | 4.4 | 4.4 | 7.4 | 6.5 | 5.7 | 4.8 | 4.7 | 11.4 | 15.8
$\rm{S_{3}}$ | 3.6 | 3.6 | 7.5 | 7.0 | 5.5 | 6.8 | 6.6 | 11.0 | 12.0
Table 1: Thickness of hinges expressed in $\mu$m corresponding to the three
considered types of structures denoted by means of $\rm{S_{1}}$, $\rm{S_{2}}$
and $\rm{S_{3}}$ abbreviations.

##  4 Results

In order to assess the possibility of changing the Poisson’s ratio of the
considered graded mechanical metamaterial, three different quasi-2D structures
presented in Fig. 1(b) are analyzed. The main difference between them is the
number of structural elements corresponding to substructures characterized by
the negative and positive Poisson’s ratio respectively. Thus, it is expected
that upon changing the ratio of structural units associated with these two
substructures, it may be possible to adjust the resultant Poisson’s ratio of
the entire structure.

###  4.1 Quasi-2D mechanical metamaterials

Figure 2: Behavior of finite quasi-2D structures subjected to the compression
along the $y$-axis. a) Mechanical deformation of Structure 1 corresponding to
the $N_{\rm{S,x}}$=1 and $N_{\rm{T,x}}$=10 parameters. b) Behavior of
Structure 2 associated with $N_{\rm{S,x}}$=2 and $N_{\rm{T,x}}$=6 parameters.
c) Compression of Structure 3 corresponding to $N_{\rm{S,x}}$=5 and
$N_{\rm{T,x}}$=1. d) Engineering Poisson’s ratio corresponding to the
compression of the considered finite structures (see Fig. 1(b)) in the $xy$
plane along the $y$ axis. e) Engineering Poisson’s ratio corresponding to
individual unit-cells with periodic boundary conditions determined through the
use of the analytical model and FEM simulations. In this case, the dashed and
solid lines correspond to the FEM simulations and the analytical model
respectively.

As shown in Fig. 2(a), in the case of the first of the analyzed structures,
i.e. Structure 1, the top part of the system which is unconstrained expands
significantly with respect to its initial dimension upon being compressed.
This in turn is an indication of the positive Poisson’s ratio. However, as
mentioned in the Experiment subsection, in order for the results not to be
strongly influenced by the edge effects at the top of the compressed system,
the second topmost unit-cell is taken into consideration to assess mechanical
properties of the structure (see highlighted unit-cell in Fig. 1(a)). As can
be seen based on the provided pictures and images extracted from FEM
simulations, the horizontal dimension of this specific unit-cell still expands
significantly with respect to its initial length. This, in addition to the
shrinkage of the unit-cell along the $y$-axis, results in the positive
Poisson’s ratio behavior. Furthermore, it should be noted that the specific
way how the Poisson’s ratio was calculated for finite systems analyzed in this
work is provided in the Supplementary Information.

According to Fig. 2(b), the behavior of the system is significantly changed in
comparison to what was observed for Structure 1. More specifically, it appears
that the change in the horizontal dimension of the structure is very small.
This is indicative of the resultant Poisson’s ratio being in the vicinity of
zero.

As presented in Fig. 2(c), the last of the analyzed finite systems, i.e.
Structure 3, exhibits a drastically different mechanical response than
structures 1 and 2. More specifically, upon having a closer look at the
generated results, it is possible to note that the horizontal dimension
(measured along the $x$-axis) of the analyzed part of the system was
significantly decreased as a result of the compression along the $y$-axis.
This in turn is indicative of the highly negative Poisson’s ratio exhibited by
the structure. Nevertheless, to gain a better understanding of the effect that
the change in a number of structural elements of individual substructures has
on the overall behavior of the resultant graded system, it is necessary to
analyze quantitative results.

According to Fig. 2(d), it is clearly visible that depending on the number of
structural elements corresponding to the two types of substructures
constituting the system, it is possible to very significantly change the
resultant Poisson’s ratio of the entire system. In fact, it is even possible
to change it from highly positive to highly negative values which is certainly
very appealing from the perspective of the future design of complex multi-
substructure graded mechanical metamaterials. Furthermore, it should be
emphasized that all observations made during the above qualitative analysis
are confirmed by the quantitative results. In addition, both experimental
results and results extracted from FEM simulations appear to be in a very good
agreement which confirms their validity. Upon having a closer look at the
generated results, it is possible to gain a better understanding of how the
individual structures behave during the compression process. In the case of
results associated with Structure 1, it is possible to observe that the
Poisson’s ratio assumes positive values throughout the entire deformation
process. More specifically, this quantity initially assumes values around
0.75. However, upon exceeding the strain of magnitude 0.02, the Poisson’s
ratio seems to maintain the value from the interval between 0.3 and 0.5.
Furthermore, in the case of Structure 2, it is possible to note that the
Poisson’s ratio initially assumes a small positive value in order to later
exhibit the Poisson’s ratio around -0.1 upon reaching larger strains. This is
a great indication that in general, one can design the considered graded
mechanical metamaterial in a way where the resultant Poisson’s ratio will
remain in the vicinity of zero which can be important in the case of
applications utilizing mechanical materials that are supposed to retain their
lateral dimension upon being deformed. In the case of the last of the analyzed
systems, i.e. Structure 3, one can note that the value of the Poisson’s ratio
is maintained in the vicinity of -0.5 throughout the entire deformation
process. This, in turn, is a very interesting result from the point of view of
devices that utilize materials that are not supposed to lose their auxeticity
throughout the deformation (e.g. protective equipment).

Even though the experimental and computational results presented in Fig. 1(d)
are very promising, it is essential to check how the considered systems would
behave if one was to form a very large structure composed of unit-cells
proposed in this work instead of using relatively small finite structures. To
this aim, the behavior of the analytical model corresponding to the highly
idealized pin-jointed hinges was compared to the FEM simulations corresponding
to the single unit-cell with realistic hinges (see Fig. 1(a)) and periodic
boundary conditions being imposed. As shown in Fig. 2(e), similarly to the
results generated for finite structures, it is possible to observe a very
large variation in the resultant Poisson’s ratio depending on the type of the
structure. However, it is also possible to note that in this case, Poisson’s
ratios associated with Structure 1 and Structure 2 assume more positive values
than was formerly the case. Similarly, the Poisson’s ratio corresponding to
Structure 3 assumes more negative values throughout the deformation process
than was observed for the formerly analyzed finite system. Finally, it is
worth noting that both analytical results associated with idealized structures
and FEM simulations lead to almost identical results. This in turn is an
indication that relatively small changes in the thicknesses of hinges used in
this work are not expected to influence the results. In addition, one can note
that all of the results presented in the main text correspond to the
compression of the analyzed samples. This stems from the fact that it is
extremely difficult to stretch micro-scale structures. However, one can
conduct such an analysis in the case of the FEM simulations. As shown in Fig.
S3 in the Supplementary Information, trends corresponding to Poisson’s ratio
exhibited by the considered systems are unaffected by the change in the
loading direction. This originates from the fact that all of the considered
samples are unimode structures, i.e. their deformation can be described
through the change in only one parameter.

###  4.2 3D structures

The concept reported in this work is not limited to 2D or quasi-2D systems and
in general can be extended to 3D structures capable of exhibiting the desired
mechanical properties in multiple directions. In order to show it, in this
work, we designed two new structures that are based on quasi-2D systems.

Figure 3: The behavior of finite 3D structures subjected to the compression
along the $y$-axis. a) SEM images of the two types of 3D structures analyzed
in this work. b) A conceptual diagram showing how the considered 3D structures
could become a part of a larger hexagonal lattice. c) Pictures taken by means
of the optical camera showing the change in the cross-section of the analyzed
structures in the $xy$ plane upon being compressed along the $y$-axis. d) The
graph showing the engineering Poisson’s ratio corresponding to the analyzed 3D
structures. In this case, the horizontal dimension used in the calculations of
the Poisson’s ratio corresponds to the topmost part of the finite system.

The considered 3D structures are called Structure A and Structure B and are
presented in Fig. 3(a). These mechanical metamaterials are based on quasi-2D
structures named Structure 1 and Structure 3 respectively that were already
proven to result in very different Poisson’s ratios. Upon taking Structure A
as an example, it should be emphasized that in order to construct it, three
indentical elements in the form of the quasi-2D Structure 1 were used that
were then spatially arranged to form a configuration resembling an equilateral
triangle (in terms of its cross-section) in the $xz$ plane. As shown in Fig.
3(a), the corners of such three blocks were also connected together by means
of additional layers of a material to ensure that they move together during
the deformation process of the entire system. Similarly, Structure B was
designed based on three block-like elements corresponding to Structure 3.

As shown in Fig. 3(b-c), as Structure A is being compressed, its horizontal
dimension increases which results in the positive Poisson’s ratio. It should
be noted that in this case, the horizontal and vertical dimensions used in the
calculation of the Poisson’s ratio correspond to the horizontal dimension of
the topmost part of the cross-section shown in Fig. 3(c) and the overall
vertical dimension respectively. Conversely, in the case of Structure B, the
system exhibits auxetic behavior throughout the deformation process. Both of
these results should be expected based on properties of quasi-2D samples
referred to as Structure 1 and Structure 3. However, one can note that for 3D
structures, the magnitude of the Poisson’s ratio is significantly smaller than
was formerly the case. This stems from the fact that the horizontal dimension
of the cross-section of the system incorporates additional non-deformable
elements that were used to connect the adjacent structural blocks together. At
this point, it is also worth noting that even though this is not the subject
of investigation in this work, the regions of maximum stress accumulation
within the structures can be ascertained based on the already conducted FEM
simulations of the behavior of 2D structures. More specifically, the maximum
stress can be observed at the hinges connecting mutually-rotating structural
elements. On the other hand, the localized stresses within the central regions
of the square-shaped or triangle-shaped structural elements are negligible.

###  4.3 Wave propagation

In addition to mechanical properties, one can also analyze the effect that the
variation in the composition of considered mechanical metamaterials has on the
wave propagation through these structures. The phonon dispersion relation was
calculated using COMSOL Multiphysics software for the three types of unit-
cells presented in Fig. 2. To generate the results, in the case of all of the
analyzed examples, Floquet periodic boundary conditions were implemented in
the $x$\- and $y$-directions assuming two-dimensional unit-cells.

The main objective of this part of the work is to determine the group
velocities associated with longitudinal and transverse waves propagating in
two principal directions, i.e. along the $x$ and $y$ axes. This will give an
idea of how does the speed of the wave propagation change upon considering
systems composed of a different number of auxetic and positive Poisson’s ratio
structural elements. This, in turn, can prove to be interesting for example in
the case of sensors capable of detecting the configuration assumed by the
system in the case when it is not easy to determine it via a direct optical
observation. Such possibility is of particular significance for micro-scale
structures where such procedure could prove to be challenging. In addition, it
should be mentioned that for each of the structures, the results were
generated for the 20 first eigenmodes. However, according to the generated
data, no band gaps were observed. Instead, as shown in Fig. 4, one can note a
significant hybridization originating from the complexity of the resultant
structures and their general design based on 2 effectively different
metamaterials. Finally, one can note that the observed results correspond to
the Megahertz range which stems from the fact that the analyzed structures
were constructed at the micro-scale. However, in general, together with the
size of the structure itself, the frequency range can be also adjusted to
match the desired wavelength. At this point, one should also note that in our
analysis, we are only focusing on the behavior of the 2D structures. This
stems from the fact that primitive cells associated with the 3D structures
portrayed in Fig. 3 have to be defined in 3D and are very complex. Thus, the
corresponding numerical simulations would be very time-consuming and
difficult. At the same time, we do not expect for the conclusions related to
the wave propagation in 3D structures to be significantly different from those
associated with 2D systems. In addition, it should be noted that the presented
results are valid for specific structures investigated in this work. However,
should someone fabricate similar samples at a different scale or through the
use of a different material, the band structure associated with the system
would be expected to be rescaled. Thus, to obtain more general results, one
can normalize the presented frequencies with respect to the lattice constant
and the P wave velocity. These results are provided in the Supplementary
Information in a form of Fig. S2 and correspond to the same trends as shown in
Fig. 4.

Figure 4: Phonon dispersion derived for unit-cells corresponding to the
considered structures. a) Unit-cells of considered structures (in real space).
b) Phonon dispersion graphs corresponding to the $x$ direction. c) Phonon
dispersion graphs associated with the $y$ directions. For all graphs, there is
always one band highlighted by means of the red and blue colors. These bands
were used in the quantitative analysis to determine the group velocities for
each of the structures. The auxiliary diagram of the first Brillouin zone is
not up to scale.

In order to determine group velocities, only the first two branches were taken
into consideration for each of the graphs presented in Fig. 4. These branches
are graphically highlighted by means of the blue and red colors that
correspond to transverse and longitudinal waves respectively that propagate in
a given direction. In addition, since the group velocity is in general defined
as $\upsilon_{g}=\frac{d\omega}{dk}$ ($\omega$ \- angular frequency, $k$ \-
wave vector), only the initially linear fragments of these branches were
considered in order to calculate respective group velocities. At this point,
it should also be noted that overall, phonon dispersion graphs for the wave
propagation in the $x$ and $y$ directions is very different. In fact, this
stems from the high anisotropy of the considered structures.

| S1 | S2 | S3
---|---|---|---
$x$ direction |  $\upsilon_{g}^{L}$ [m/s] | 596.4 | 620.0 | 702.9
$\upsilon_{g}^{T}$ [m/s] | 520.6 | 514.6 | 460.2
$y$ direction |  $\upsilon_{g}^{L}$ [m/s] | 781.9 | 567.9 | 574.3
$\upsilon_{g}^{T}$ [m/s] | 515.0 | 520.6 | 461.8
Table 2: Values of group velocities $\upsilon_{g}^{L}$ and $\upsilon_{g}^{T}$
corresponding to longitudinal and transverse waves respectively. Abbreviations
S1, S2 and S3 stand for Structure 1, Structure 2 and Structure 3.

As shown in Table 2, in the case of waves propagating in the $x$-direction,
the group velocity associated with longitudinal waves ($\upsilon_{g}^{L}$)
seems to be the smallest for Structure 1, in order to reach the highest value
of 702.9 m/s in the case of Structure 3. This means that in this specific
direction, it is possible to observe a significant change in
$\upsilon_{g}^{L}$ upon changing the composition of the structure. Conversely,
in the case of transverse waves propagating in the $x$-direction, the group
velocity $\upsilon_{g}^{T}$ gradually decreases with a relatively small change
in the magnitude. Furthermore, it should be emphasized that for waves
propagating in the $y$-direction, it is hard to observe any specific trends
related to the change in the magnitude of the group velocities corresponding
either to longitudinal or transverse waves. In addition, as shown in Fig.
4(b), the branches form a rather perturbed picture in comparison to what could
have been observed for waves propagating in the $x$-direction. This can be
explained by the fact that in the $y$-direction, the substructures
constituting the considered systems form parallel stripes corresponding to a
very different topology. Hence, it is possible to expect a strong interference
that influences the results. Furthermore, it should be noted that trends
presented in Fig. 4 and Table 2 can also be observed if one was to replace the
two types of substructures considered in this work with bulk materials having
the same external dimensions and effective mechanical properties (see
Supplementary Information). This indicates that the presented results are
reliable and that similar types of the complex phonon dispersion relations
could potentially be observed for other multi-substructure graded structures.
Last but not least, we note that both longitudinal and transverse bands for
Structure 1 show some very important non-monotonic behavior where its form was
recently reported as the roton effect in metamaterials. The origin of this
effect was attributed to the third-order interactions in equivalent mass-
springs models Chen et al. (2021b); Iglesias Martínez et al. (2021). Finally,
is should be noted that further details related to the analysis of the
observed results are provided in the Supplementary Information. Amongst
others, the phonon dispersion relation corresponding solely to geometries
related to rotating squares and rotating triangles mechanical metamaterials
are presented. This, in turn provides, an information about the wave
propagation properties of materials constituting individual part of the more
complex graded structure considered in this work.

##  5 Discussion

The novel graded mechanical metamaterial proposed in this work offers numerous
advantages over a majority of known mechanical metamaterials and other graded
structures. The considered graded system incorporates two structures that can
have their deformation fully described by a change in one parameter assuming
relative rigidity of their structural elements. As a result, the evolution of
the entire system investigated in this work can be also described by a
variation in one geometric parameter which distinguishes it from many other
graded structures Mirzaali et al. (2020); Guo et al. (2021b); Han and Wei
(2022); Niknam and Akbarzadeh (2020). Such unique design, also results in the
very high stability of the system which is often not observed for other graded
metamaterials composed for example of different re-entrant structural motifs.
In addition, the two types of substructures present within the system are
characterized by drastically different Poisson’s ratio. Hence, depending on
the composition of the considered mechanical metamaterial, it can exhibit
either very positive or very negative values of this quantity. This, in turn,
as demonstrated in another study Dudek et al. (2022b), can enable the
structure to exhibit a complex shape-morphing behavior. Finally, it should be
mentioned that in terms of the manufacturing process, the proposed design
makes it possible to 3D print different complex components of the structure at
the micro-scale with a very high precision. However, In the case of many other
designs known in the literature, this task would become significantly more
challenging due to the presence of overhanging elements.

All of the results presented in this work suggest that the considered
mechanical metamaterials can exhibit a very broad range of the Poisson’s ratio
including both highly negative and positive values depending solely on the
ratio of the number of structural elements corresponding to the negative and
positive Poisson’s ratio part of the system. It means that within certain
limits, one can achieve an arbitrary value of this parameter based solely on
two types of substructures. This, in turn, is of great significance as it may
enable engineers to produce relatively simple mechanical materials that can
exhibit desired mechanical properties without the need of utilizing very
complex and unique designs for each type of application. It is also worth
mentioning that the latter approach is relatively expensive and may lead to
malfunctions due to the unnecessary complexity of structural designs. In
addition, in this study, it was demonstrated that the considered results can
be observed in the case of microscopic systems. This means that in addition to
applications typical for many auxetic mechanical metamaterials, e.g.
protective devices or functional materials utilizing an enhanced fracture
toughness, indentation resistance as well as strain-dependent fully-reversible
porosity and deformation, this concept may also be implemented in the case of
novel biomedical devices and flexible electronics. Particularly important in
this regard is the potential ability of structures presented in this work to
exhibit complex shape morphing upon connecting multiple elements similar to
considered structures in order to form larger lattices. As a result, one can
obtain a graded metamaterial capable of exhibiting very different local values
of the Poisson’s ratio which, as presented in the recent studies devoted to
micro-scale shape-morphing mechanical metamaterials Dudek et al. (2022b), may
lead to the reconfiguration of the initial system into a complex predefined
shape. Such ability could for example prove to be useful in the design of
efficient stents where upon being subjected to the external stimulus, devices
based on the proposed concept could provide a better support to specific parts
of the blood vessel.

##  6 Conclusion

In this work, novel 2D and 3D functionally graded mechanical metamaterials
were proposed that consist of two substructures corresponding to the negative
and positive Poisson’s ratio. It was shown that depending on the number of
structural elements corresponding to each type of the substructure, it is
possible to significantly change the resultant value of Poisson’s ratio to
reach its desired magnitude. It was also shown that structures investigated in
this work can be constructed at the microscale in addition to larger scales
that are typical for mechanical metamaterials. This means that the proposed
concept could be used for example in the design of the biomedical devices such
as stents where the proposed solutions could enhance the efficiency of the
already existing devices by locally supporting specific parts of the blood
vessel at scales that normally cannot be achieved in the case of mechanical
metamaterials. In addition, results reported in this study could also be
utilized at a larger scale in the design of efficient protective devices or
sports equipment corresponding to a specific value of the Poisson’s ratio that
is beneficial from the perspective of the specific application. Finally, it
was demonstrated that the change in the composition of considered structures
affects the wave propagation through the system. This in turn may be
interesting from the point of view of artificial materials capable of
exhibiting a local wave attenuation.

## Declaration of competing interest

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

## 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. aa

## Data Availability

The data presented in this work are available from the corresponding authors
uopn reasonable request.

## References

  * Albertini et al. [2021] Albertini, F., Dirrenberger, J., Sollogoub, C., Maconachie, T., Leary, M., Molotnikov, A., 2021\.  Experimental and computational analysis of the mechanical properties of composite auxetic lattice structures.  Addit. Manuf. 47, 102351\.
  * An et al. [2020] An, N., Domel, A.G., Zhou, J., Rafsanjani, A., Bertoldi, K., 2020\.  Programmable hierarchical kirigami.  Adv. Funct. Mater. 30, 1906711\.
  * Babaee et al. [2013] Babaee, S., Shim, J., Weaver, J.C., Chen, E.R., Patel, N., Bertoldi, K., 2013\.  3d soft metamaterials with negative poisson’s ratio.  Adv. Mater. 25, 5044–5049.
  * Baughman et al. [1998] Baughman, R.H., Stafstrom, S., Cui, C., Dantas, S.O., 1998\.  Materials with negative compressibilities in one or more dimensions.  Science 279, 1522–1524.
  * Bertoldi et al. [2017] Bertoldi, K., Vitelli, V., Christensen, J., van Hecke, M., 2017\.  Flexible mechanical metamaterials.  Nat. Rev. Mater. 2, 17066\.
  * Boldrin et al. [2016] Boldrin, L., Hummel, S., Scarpa, F., Di Maio, D., Lira, C., Ruzzene, M., Remillat, C., Lim, T., Rajasekaran, R., Patsias, S., 2016\.  Dynamic behaviour of auxetic gradient composite hexagonal honeycombs.  Compos. Struct. 149, 114–124.  doi:https://doi.org/10.1016/j.compstruct.2016.03.044.
  * Chen et al. [2021a] Chen, S., Tan, X., Hu, J., Zhu, S., Wang, B., Wang, L., Jin, Y., Wu, L., 2021a.  A novel gradient negative stiffness honeycomb for recoverable energy absorption.  Compos. B. Eng. 215, 108745\.
  * Chen et al. [2020] Chen, Y., Frenzel, T., Guenneau, S., Kadic, M., Wegener, M., 2020.  Mapping acoustical activity in 3d chiral mechanical metamaterials onto micropolar continuum elasticity.  J. Mech. Phys. Solids 137, 103877\.
  * Chen et al. [2021b] Chen, Y., Kadic, M., Wegener, M., 2021b.  Roton-like acoustical dispersion relations in 3d metamaterials.  Nat. Commun. 12, 1–8.
  * Cho et al. [2014] Cho, Y., Shin, J.H., Costa, A., T. A. Kim, V.K., Li, J., Lee, S.Y., Yang, S., Han, H.N., Choi, I.S., Srolovitz, D.J., 2014.  Engineering the shape and structure of materials by fractal cut.  Proc. Natl. Acad. Sci. 111, 17390–17395.
  * Christensen et al. [2015] Christensen, J., Kadic, M., Kraft, O., Wegener, M., 2015\.  Vibrant times for mechanical metamaterials.  Mrs Commun. 5, 453–462.
  * Coulais et al. [2018] Coulais, C., Sabbadini, A., Vink, F., van Hecke, M., 2018\.  Multi-step self-guided pathways for shape-changing metamaterials.  Nature 561, 512–515.
  * Coulais et al. [2016] Coulais, C., Teomy, E., de Reus, K., Shokef, Y., van Hecke, M., 2016.  Combinatorial design of textured mechanical metamaterials.  Nature 535, 529–532.
  * Dudek et al. [2016] Dudek, K.K., Attard, D., Caruana-Gauci, R., Wojciechowski, K.W., Grima, J.N., 2016.  Unimode metamaterials exhibiting negative linear compressibility and negative thermal expansion.  Smart Mater. Struct. 25, 025009\.
  * Dudek et al. [2022a] Dudek, K.K., Iglesias Martínez, J.A., Kadic, M., 2022a.  Variable dual auxeticity of the hierarchical mechanical metamaterial composed of re-entrant structural motifs.  Phys. Status Solidi B , 2200404.
  * Dudek et al. [2022b] Dudek, K.K., Iglesias Martínez, J.A., Ulliac, G., Kadic, M., 2022b.  Micro-scale auxetic hierarchical mechanical metamaterials for shape morphing.  Adv. Mater. 34, 2110115\.
  * Dudek et al. [2023] Dudek, K.K., Iglesias Martínez, J.A., Ulliac, G., Hirsinger, L., Wang, L., Laude, V., Kadic, M., 2023.  Micro-scale mechanical metamaterial with a controllable transition in the poisson’s ratio and band gap formation.  Adv. Mater. n/a, 2210993\.  doi:https://doi.org/10.1002/adma.202210993.
  * Dudek et al. [2019] Dudek, K.K., Wolak, W., Gatt, R., Grima, J.N., 2019\.  Impact resistance of composite magnetic metamaterials.  Sci. Rep. 9, 3963\.
  * Duncan et al. [2018] Duncan, O., Shepherd, T., Moroney, C., Foster, L., Venkatraman, P.D., Winwood, K., Allen, T., Alderson, A., 2018\.  Review of auxetic materials for sports applications: Expanding options in comfort and protection.  Appl. Sci. 8, 941\.
  * Evans and Alderson [2000] Evans, K.E., Alderson, A., 2000\.  Auxetic materials: Functional materials and structures from lateral thinking!  Adv. Mater. 12, 617–628.
  * Fleury et al. [2016] Fleury, R., Khanikaev, A.B., Alù, A., 2016.  Floquet topological insulators for sound.  Nat. Commun. 7, 11744\.
  * Fleury et al. [2014] Fleury, R., Sounas, D.L., Sieck, C.F., Haberman, M.R., Alù, A., 2014.  Sound isolation and giant linear nonreciprocity in a compact acoustic circulator.  Science 343, 516–519.
  * Florijn et al. [2014] Florijn, B., Coulais, C., van Hecke, M., 2014.  Programmable mechanical metamaterials.  Phys. Rev. Lett. 113, 175503\.
  * Frenzel et al. [2017] Frenzel, T., Kadic, M., Wegener, M., 2017.  Three-dimensional mechanical metamaterials with a twist.  Science 358, 1072–1074.
  * Galea et al. [2022] Galea, R., Dudek, K.K., Farrugia, P.S., Mangion, L.Z., Grima, J.N., Gatt, R., 2022\.  Reconfigurable magneto-mechanical metamaterials guided by magnetic fields.  Compos. Struct. 280, 114921\.
  * Gao et al. [2019] Gao, J., Xue, H., Gao, L., Luo, Z., 2019.  Topology optimization for auxetic metamaterials based on isogeometric analysis.  Comput. Methods. Appl. Mech. Eng. 352, 211–236.
  * Gatt et al. [2015] Gatt, R., Mizzi, L., Azzopardi, J.I., Azzopardi, K.M., Attard, D., Casha, A., Briffa, J., Grima, J.N., 2015\.  Hierarchical auxetic mechanical metamaterials.  Sci. Rep. 5, 8395\.
  * Grima et al. [2011] Grima, J.N., Chetcuti, E., Manicaro, E., Attard, D., Camilleri, M., Gatt, R., Evans, K.E., 2011.  On the auxetic properties of generic rotating rigid triangles.  Proc. R. Soc. A. 468, 810–830.
  * Grima and Evans [2000] Grima, J.N., Evans, K.E., 2000\.  Auxetic behavior from rotating squares.  J. Mater. Sci. Lett. 19, 1563–1565.
  * Guenneau et al. [2007] Guenneau, S., Movchan, A., Petursson, G., Ramakrishna, S.A., 2007\.  Acoustic metamaterials for sound focusing and confinement.  New J. Phys. 9, 399\.
  * Guo et al. [2021a] Guo, M.F., Yang, H., Zhou, Y.M., Ma, L., 2021a.  Mechanical properties of 3d hybrid double arrow-head structure with tunable poisson’s ratio.  Aerosp. Sci. Technol. 119, 107177\.
  * Guo et al. [2021b] Guo, W., Huang, Y., Ritchie, R.O., Yin, S., 2021b.  Dissipative dual-phase mechanical metamaterial composites via architectural design.  Extreme Mech. Lett. 48, 101442\.
  * Hamzehei et al. [2020] Hamzehei, R., Rezaei, S., Kadkhodapour, J., Anaraki, A.P., Mahmoudi, A., 2020.  2d triangular anti-trichiral structures and auxetic stents with symmetric shrinkage behavior and high energy absorption.  Mech. Mater. 142, 103291\.  doi:https://doi.org/10.1016/j.mechmat.2019.103291.
  * Han and Wei [2022] Han, Z., Wei, K., 2022.  Multi-material topology optimization and additive manufacturing for metamaterials incorporating double negative indexes of poisson’s ratio and thermal expansion.  Addit. Manuf. 54, 102742\.
  * Hewage et al. [2016] Hewage, T.A.M., Alderson, K.L., Alderson, A., Scarpa, F., 2016\.  Double-negative mechanical metamaterials displaying simultaneous negative stiffness and negative poisson’s ratio properties.  Adv. Mater. 28, 10323–10332.
  * Hou et al. [2014] Hou, Y., Neville, R., Scarpa, F., Remillat, C., Gu, B., Ruzzene, M., 2014\.  Graded conventional-auxetic kirigami sandwich structures: Flatwise compression and edgewise loading.  Compos. B. Eng. 59, 33–42.
  * Iglesias Martínez et al. [2021] Iglesias Martínez, J.A., Groß, M.F., Chen, Y., Frenzel, T., Laude, V., Kadic, M., Wegener, M., 2021.  Experimental observation of roton-like dispersion relations in metamaterials.  Sci. Adv. 7, eabm2189.
  * Imbalzano et al. [2018] Imbalzano, G., Linforth, S., Ngo, T.D., Lee, P.V.S., Tran, P., 2018.  Blast resistance of auxetic and honeycomb sandwich panels: Comparisons and parametric designs.  Compos. Struct. 183, 242–261.
  * Jiang et al. [2020] Jiang, H., Barbenchon, L.L., Bednarcyk, B.A., Scarpa, F., Chen, Y., 2020.  Bioinspired multilayered cellular composites with enhanced energy absorption and shape recovery.  Addit. Manuf. 36, 101430\.
  * Jiang et al. [2022] Jiang, W., Ren, X., Wang, S.L., Zhang, X.G., Zhang, X.Y., Luo, C., Xie, Y.M., Scarpa, F., Alderson, A., Evans, K.E., 2022.  Manufacturing, characteristics and applications of auxetic foams: A state-of-the-art review.  Compos. B. Eng., 235, 109733\.
  * Kadic et al. [2019] Kadic, M., Milton, G.W., van Hecke, M., Wegener, M., 2019\.  3d metamaterials.  Nat. Rev. Phys. 1, 198–210.
  * Kolken et al. [2018a] Kolken, H.M.A., Janbaz, S., Leeflang, S.M.A., Lietaert, K., Weinans, H.H., Zadpoor, A.A., 2018a.  Rationally designed meta-implants: A combination of auxetic and conventional meta-biomaterials.  Mater. Horiz. 5, 28–35.
  * Kolken et al. [2018b] Kolken, H.M.A., Janbaz, S., Leeflang, S.M.A., Lietaert, K., Weinans, H.H., Zadpoor, A.A., 2018b.  Rationally designed meta-implants: a combination of auxetic and conventional meta-biomaterials.  Mater. Horiz. 5, 28–35.
  * Korpas et al. [2021] Korpas, L.M., Yin, R., Yasuda, H., Raney, J.R., 2021\.  Temperature-responsive multistable metamaterials.  ACS Appl. Mater. Interfaces 13, 31163–31170.
  * Lakes [1987] Lakes, R., 1987.  Foam structures with a negative poisson’s ratio.  Science 235, 1038\.
  * Lakes [1993] Lakes, R., 1993.  Materials with structural hierarchy.  Nature 361, 511–515.
  * Lakes and Wojciechowski [2008] Lakes, R., Wojciechowski, K.W., 2008\.  Negative compressibility, negative poisson’s ratio, and stability.  Phys. Status Solidi B 245, 545–551.
  * Lim [2019] Lim, T.C., 2019.  Composite microstructures with poisson’s ratio sign switching upon stress reversal.  Compos. Struct. 209, 34–44.
  * Liu et al. [2020] Liu, K., Han, L., Hu, W., Ji, L., Zhu, S., Wan, Z., Yang, X., Wei, Y., Dai, Z., Zhao, Z., Li, Z., Wang, P., Tao, R., 2020\.  4d printed zero poisson’s ratio metamaterial with switching function of mechanical and vibration isolation performance.  Mater. Des. 196, 109153\.
  * van Manen et al. [2021] van Manen, T., Janbaz, S., Jansen, K.M.B., Zadpoor, A.A., 2021\.  4d printing of reconfigurable metamaterials anddevices.  Commun. Mater. 2, 56\.
  * Milton [2013] Milton, G.W., 2013.  Complete characterization of the macroscopic deformations of periodic unimode metamaterials of rigid bars and pivots.  J. Mech. Phys. 61, 1543–1560.
  * Miniaci et al. [2016] Miniaci, M., Krushynska, A., Bosia, F., Pugno, N.M., 2016\.  Large scale mechanical metamaterials as seismic shields.  New J. Phys. 18, 083041\.
  * Mirzaali et al. [2018] Mirzaali, M.J., Janbaz, S., Strano, M., Vergani, L., Zadpoor, A.A., 2018.  Shape-matching soft mechanical metamaterials.  Sci. Rep. 8, 965\.
  * Mirzaali et al. [2020] Mirzaali, M.J., Pahlavani, H., Yarali, E., Zadpoor, A.A., 2020\.  Non-affinity in multi-material mechanical metamaterials.  Sci. Rep. 10, 11488\.
  * Mizzi et al. [2021] Mizzi, L., Attard, D., Gatt, R., Dudek, K.K., Ellul, B., Grima, J.N., 2021\.  Implementation of periodic boundary conditions for loading of mechanical metamaterials and other complex geometric microstructures using finite element analysis.  Eng. Comput. 37, 1765–1779.
  * Mizzi et al. [2018] Mizzi, L., Mahdi, E., Titov, K., Gatt, R., Attard, D., Evans, K.E., Grima, J.N., Tan, J.C., 2018\.  Mechanical metamaterials with star-shaped pores exhibiting negative and zero poisson’s ratio.  Mater. Des. 146, 28–37.
  * Neville et al. [2016] Neville, R.M., Scarpa, F., Pirrera, A., 2016.  Shape morphing kirigami mechanical metamaterials.  Sci. Rep. 6, 31067\.
  * Nicolaou and Motter [2012] Nicolaou, Z.G., Motter, A.E., 2012\.  Mechanical metamaterials with negative compressibility transitions.  Nat. Mater. 11, 608–613.
  * Niknam and Akbarzadeh [2020] Niknam, H., Akbarzadeh, A., 2020\.  Graded lattice structures: Simultaneous enhancement in stiffness and energy absorption.  Mater. Des. 196, 109129\.  doi:https://doi.org/10.1016/j.matdes.2020.109129.
  * Ning et al. [2018] Ning, X., Yu, X., Wang, H., Sun, R., Corman, R.E., Li, H., Lee, C.M., Xue, Y., Chempakasseril, A., Yao, Y., Zhang, Z., Luan, H., Wang, Z., Xia, W., Feng, X., Ewoldt, R.H., Huang, Y., Zhang, Y., Rogers, J.A., 2018\.  Mechanically active materials in three-dimensional mesostructures.  Sci. Adv. 4, eaat8313.  doi:10.1126/sciadv.aat8313.
  * Novak et al. [2021a] Novak, N., Biasetto, L., Rebesan, P., Zanini, F., Carmignato, S., Krstulović-Opara, L., Vesenjak, M., Ren, Z., 2021a.  Experimental and computational evaluation of tensile properties of additively manufactured hexa- and tetrachiral auxetic cellular structures.  Addit. Manuf. 45, 102022\.
  * Novak et al. [2021b] Novak, N., Duncan, O., Allen, T., Alderson, A., Vesenjak, M., Ren, Z., 2021b.  Shear modulus of conventional and auxetic open-cell foam.  Mech. Mater. 157, 103818\.
  * Oftadeh et al. [2014] Oftadeh, R., Haghpanah, B., Vella, D., Boudaoud, A., Vaziri, A., 2014.  Optimal fractal-like hierarchical honeycombs.  Phys. Rev. Lett. 113, 104301\.
  * Pagliocca et al. [2022] Pagliocca, N., Uddin, K.Z., Anni, I.A., Shen, C., Youssef, G., Koohbor, B., 2022\.  Flexible planar metamaterials with tunable poisson’s ratios.  Mater. Des. 215, 110446\.  doi:https://doi.org/10.1016/j.matdes.2022.110446.
  * Placidi et al. [2014] Placidi, L., Rosi, G., Giorgio, I., Madeo, A., 2014\.  Reflection and transmission of plane waves at surfaces carrying material properties and embedded in second-gradient materials.  Math. Mech. Solids 19, 555–578.
  * Quadrelli et al. [2021] Quadrelli, D.E., Craster, R., Kadic, M., Braghin, F., 2021\.  Elastic wave near-cloaking.  Extreme Mech. Lett. 44, 101262\.
  * Rafsanjani et al. [2018] Rafsanjani, A., Zhang, Y., Liu, B., Rubinstein, S.M., Bertoldi, K., 2018.  Kirigami skins make a simple soft actuator crawl.  Sci. Robot. 3, eaar7555.
  * Ruzzene et al. [2002] Ruzzene, M., Mazzarella, L., Tsopelas, P., Scarpa, F., 2002\.  Wave propagation in sandwich plates with periodic auxetic core.  J. Intell. Mater. Syst. Struct. 13, 587–597.
  * Shao et al. [2021] Shao, Y., Meng, J., Ma, G., Ren, S., Fang, L., Cao, X., Liu, L., Li, H., Wu, W., Xiao, D., 2021\.  Insight into the negative poisson’s ratio effect of the gradient auxetic reentrant honeycombs.  Compos. Struct. 274, 114366\.  doi:https://doi.org/10.1016/j.compstruct.2021.114366.
  * Sun et al. [2020] Sun, R., Li, Q., Yao, J., Scarpa, F., Rossiter, J., 2020\.  Tunable, multi-modal, and multi-directional vibration energy harvester based on three-dimensional architected metastructures.  Appl. Energy 264, 114615\.  URL: https://www.sciencedirect.com/science/article/pii/S0306261920301276, doi:https://doi.org/10.1016/j.apenergy.2020.114615.
  * Tan et al. [2022] Tan, X., Wang, B., L. Wang, S.Z., Chen, S., Yao, K., Xu, P., 2022.  Effect of beam configuration on its multistable and negative stiffness properties.  Compos. Struct. 286, 115308\.
  * Tang et al. [2015] Tang, Y., Lin, G., Han, L., Qiu, S., Yang, S., Yin, J., 2015.  Design of hierarchically cut hinges for highly stretchable and reconfigurable metamaterials with enhanced strength.  Adv. Mater. 27, 7181–7190.
  * Wang et al. [2020] Wang, H., Zhao, D., Jin, Y., Wang, M., Mukhopadhyay, T., You, Z., 2020.  Modulation of multi-directional auxeticity in hybrid origami metamaterials.  Appl. Mater. Today 20, 100715\.
  * Wang et al. [2016] Wang, K., Wu, C., Qian, Z., Zhang, C., Wang, B., Vannan, M.A., 2016.  Dual-material 3d printed metamaterials with tunable mechanical properties for patient-specific tissue-mimicking phantoms.  Addit. Manuf. 12, 31–37.
  * Wang et al. [2014] Wang, Y., Luo, Z., Zhang, N., Kang, Z., 2014.  Topological shape optimization of microstructural metamaterials using a level set method.  Comput. Mater. Sci. 87, 178–186.
  * Wenz et al. [2021] Wenz, F., Schmidt, I., Leichner, A., Lichti, T., Baumann, S., Andrae, H., Eberl, C., 2021.  Designing shape morphing behavior through local programming of mechanical metamaterials.  Adv. Mater. 33, 2008617\.
  * Wojciechowski [1989] Wojciechowski, K.W., 1989.  Two-dimensional isotropic system with a negative poisson ratio.  Phys. Lett. A 137, 60–64.
  * Wu et al. [2020] Wu, L., Wang, Y., Zhai, Z., Yang, Y., Krishnaraju, D., Lu, J., Wu, F., Wang, Q., Jiang, H., 2020.  Mechanical metamaterials for full-band mechanical wave shielding.  Appl. Mater. Today 20, 100671\.
  * Xin et al. [2020] Xin, X., Liu, L., Liu, Y., Leng, J., 2020.  4d printing auxetic metamaterials with tunable, programmable, and reconfigurable mechanical properties.  Adv. Funct. Mater. 30, 2004226\.
  * Xu et al. [2015] Xu, S., Yan, Z., Jang, K.I., Huang, W., Fu, H., Kim, J., Wei, Z., Flavin, M., McCracken, J., Wang, R., Badea, A., Liu, Y., Xiao, D., Zhou, G., Lee, J., Chung, H.U., Cheng, H., Ren, W., Banks, A., Li, X., Paik, U., Nuzzo, R.G., Huang, Y., Zhang, Y., Rogers, J.A., 2015\.  Assembly of micro/nanomaterials into complex, three-dimensional architectures by compressive buckling.  Science 347, 154–159.  doi:10.1126/science.1260960.
  * Xue et al. [2021] Xue, H., Saha, S.C., Beier, S., Jepson, N., Luo, Z., 2021.  Topological optimization of auxetic coronary stents considering hemodynamics.  Front Bioeng Biotechnol. 9, 728914\.
  * Yao et al. [2022] Yao, J., Su, Y., Scarpa, F., Li, Y., 2022.  An optimization approach to design deformation patterns in perforated mechanical metamaterials using distributions of poisson’s ratio-based unit cells.  Compos. Struct. 281, 115015\.
  * Yao et al. [2021] Yao, J., Sun, R., Scarpa, F., Remillat, C., Gao, Y., Su, Y., 2021.  Two-dimensional graded metamaterials with auxetic rectangular perforations.  Compos. Struct. 261, 113313\.  doi:https://doi.org/10.1016/j.compstruct.2020.113313.
  * Zhang et al. [2021a] Zhang, W., Zhao, S., Sun, R., Scarpa, F., Wang, J., 2021a.  Mechanical properties of a hybrid auxetic metamaterial and metastructure system.  Journal of Reinforced Plastics and Composites 40, 785–799.  doi:10.1177/07316844211009599.
  * Zhang et al. [2021b] Zhang, Z., Scarpa, F., Bednarcyk, B.A., Chen, Y., 2021b.  Harnessing fractal cuts to design robust lattice metamaterials for energy dissipation.  Addit. Manuf. 46, 102126\.
  * Zheng et al. [2022] Zheng, X., Uto, K., Hu, W.H., Chen, T.T., Naito, M., Watanabe, I., 2022.  Reprogrammable flexible mechanical metamaterials.  Appl. Mater. Today 29, 101662\.