††thanks: equal contribution††thanks: equal contribution

# Two-dimensional higher-order topology in monolayer graphdiyne

Eunwoo Lee  Department of Physics and Astronomy, Seoul National University,
Seoul 08826, Korea  Center for Correlated Electron Systems, Institute for
Basic Science (IBS), Seoul 08826, Korea  Center for Theoretical Physics (CTP),
Seoul National University, Seoul 08826, Korea  Rokyeon Kim  Department of
Physics and Astronomy, Seoul National University, Seoul 08826, Korea  Center
for Correlated Electron Systems, Institute for Basic Science (IBS), Seoul
08826, Korea  Korea Institute for Advanced Study (KIAS), Seoul 02455, Korea
Junyeong Ahn  Department of Physics and Astronomy, Seoul National University,
Seoul 08826, Korea  Center for Correlated Electron Systems, Institute for
Basic Science (IBS), Seoul 08826, Korea  Center for Theoretical Physics (CTP),
Seoul National University, Seoul 08826, Korea  Bohm-Jung Yang
<EMAIL_ADDRESS>Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Korea  Center for Correlated Electron Systems,
Institute for Basic Science (IBS), Seoul 08826, Korea  Center for Theoretical
Physics (CTP), Seoul National University, Seoul 08826, Korea

###### Abstract

Based on first-principles calculations and tight-binding model analysis, we
propose monolayer graphdiyne as a candidate material for a two-dimensional
higher-order topological insulator protected by inversion symmetry. Despite
the absence of chiral symmetry, the higher-order topology of monolayer
graphdiyne is manifested in the filling anomaly and charge accumulation at two
corners. Although its low energy band structure can be properly described by
the tight-binding Hamiltonian constructed by using only the $p_{z}$ orbital of
each atom, the corresponding bulk band topology is trivial. The nontrivial
bulk topology can be correctly captured only when the contribution from the
core levels derived from $p_{x,y}$ and $s$ orbitals are included, which is
further confirmed by the Wilson loop calculations. We also show that the
higher-order band topology of a monolayer graphdyine gives rise to the
nontrivial band topology of the corresponding three-dimensional material, ABC-
stacked graphdiyne, which hosts monopole nodal lines and hinge states.

##  I Introduction

Bulk-boundary correspondence is a fundamental property of topological phases.
In conventional topological insulators, the gapped bulk states in
$d$-dimensions support metallic states in $(d-1)$-dimensional surfaces. Such a
conventional bulk-boundary correspondence, however, is violated in a class of
topological crystalline insulators, so-called higher-order topological
insulators (HOTIs). The gapless excitations of a HOTI in $d$-dimensions are
localized, instead, in a subspace with a dimension lower than $(d-1)$, such as
corners or hinges, when the global shape of the material preserves the
crystalline symmetry relevant to the nontrivial band topology
schindler2018higher2 ; schindler2018higher ; wang2019higher ; yue2019symmetry
; wieder2018axion ; ahn2018higher ; benalcazar2017quantized ;
langbehn2017reflection ; geier2018second ; khalaf2018higher ;
benalcazar2019quantization . Recently, rhombohedral bismuth has been
identified as the first example of 3D HOTIs hosting helical hinge states
schindler2018higher2 . Also there are other candidate materials of 3D HOTIs
including SnTe with strain along the (100) direction schindler2018higher ,
transition metal dichalcogenides MoTe2 and WTe2, hosting helical hinge states
wang2019higher , and also Bi2-xSmxSe3 with chiral hinge states yue2019symmetry
.

In 2D, one can also consider second-order topological insulators (SOTIs)
hosting localized corner states benalcazar2017quantized ;
langbehn2017reflection ; khalaf2018higher ; geier2018second ;
benalcazar2019quantization ; song2017d . Contrary to 3D HOTIs, however, there
are only few theoretical proposals for candidate materials realizing 2D SOTIs
ahn2019failure ; khalaf2018higher ; benalcazar2019quantization ;
lee2019fractional ; geier2018second . For instance, phosphorene is one
candidate for 2D HOTIs proposed based on simple tight-binding model
calculations ezawa2018minimal . However, in this system, the corner states
originate from the charge polarization ezawa2018minimal , which is a 1D
topological invariant. Hence, in a strict sense, phosphorene is not a genuine
2D HOTI. Another candidate is twisted bilayer graphene, which can also exhibit
higher-order band topology at small twist angles. However, to realize the
corner states, the twist angle as well as the chemical potential should be
controlled simultaneously ahn2019failure .

In general, among the higher-order TIs, only $d$-dimensional $k$th-order TIs
with $1<k<d$ have stable band topology with Wannier obstruction while
$d$-dimensional $d$th-order TIs with zero-dimensional corner charges are
obstructed atomic insulators. Thus, in 3D, a SOTI has Wannier obstruction
while a third-order TI (TOTI) is atomic. On the other hand, in 2D, HOTIs are
always obstructed atomic insulators. In $d$-dimensional $d$th-order TIs,
additional chiral or particle-hole symmetry is required to protect corner
states which appear as in-gap states separated from bulk states
khalaf2018higher ; geier2018second ; benalcazar2019quantization ;
langbehn2017reflection . Since chiral or particle-hole symmetry usually does
not exist in insulators, realizing 2D SOTIs in materials is generally
considered difficult.

However, even if the chiral symmetry of 2D SOTIs is broken, the material still
inherits the nontrivial band topology of 2D SOTIs. The resulting state is
called a 2D TOTI cualuguaru2019higher ; okuma2019topological ;
hwang2019fragile . The key idea is that even if the in-gap states are merged
into the bulk states, the nontrivial band topology manifests in the filling
anomaly wieder2018axion ; benalcazar2019quantization ; song2017d . Namely, the
half-filling condition can never be satisfied as long as the symmetry
associated with the second-order band topology is preserved. The half-filling
condition can be satisfied only when extra electrons or holes are added, which
naturally leads to charge accumulation or depletion at corners. Therefore, the
hallmark of 2D HOTIs lacking chiral symmetry (or simply 2D TOTIs), is the
existence of filling anomaly at half-filling and the accumulation or depletion
of a half-integral charge at corners after adding extra electrons or holes
cualuguaru2019higher ; okuma2019topological ; hwang2019fragile .

In this work, we theoretically propose monolayer graphdiyne (MGD) as the first
realistic candidate material for 2D HOTIs. Namely, MGD is a 2D SOTI when
chiral symmetry exists, but in real systems, it appears as a 2D TOTI due to
the lack of chiral symmetry. Based on first-principles calculations and tight-
binding model analysis, we show that a MGD is a HOTI characterized by a 2D
topological invariant $w_{2}$, that is quantized when the system is invariant
under inversion $P$ symmetry. Although corner states are not guaranteed to be
in the gap due to the lack of chiral symmetry in MGD, the nontrivial higher-
order band topology is manifested in filling anomaly and the corner charge
accumulation. Also, it is found that although the low energy excitation can be
well-described by using the state derived from the $p_{z}$ orbitals of each
atom, the relevant band topology of such a simplified model is trivial. To
describe the nontrivial band topology, it is essential to include the
contribution from the core levels derived from $p_{x,y}$ and $s$ orbitals.
Finally, it is shown that the non-trivial $w_{2}$ of MGD gives rise to the
topological band structure of the corresponding three-dimensional material,
ABC stacked graphdiyne, hosting nodal lines with $Z_{2}$ monopole charges
nomura2018three ; ahn2018band . Due to the nonzero $w_{2}$ in a range of 2D
momentum subspaces with fixed vertical momentum $k_{z}$, the ABC stacked
graphdiyne host hinge states for the corresponding range of $k_{z}$.

##  II results

###  II.1 Topological invariant of monolayer graphdiyne

MGD is an experimentally realized planar carbon system composed of
$sp^{2}$-$sp$ carbon network of benzene rings connected by ethynyl chains
haley1997carbon ; li2010architecture . Fig. 1(a) describes the unit cell of
MGD composed of 18 carbon atoms. Since spin-orbit coupling (SOC) is
negligible, MGD can be regarded as a spinless fermion system. The lattice has
inversion $P$, time-reversal $T$, a two-fold rotation about the $x$-axis
$C_{2x}$, a six-fold rotation about the $z$-axis $C_{6z}$, and a mirror
$M_{z}:z\rightarrow-z$ symmetries. Also, due to the bipartite lattice
structure, MGD has chiral (or sublattice) symmetry when only the nearest
neighbor hopping between different sublattices is considered, similar to the
case of graphene.

Fig. 1(b) shows the band structure along high symmetry directions obtained by
first-principles calculations. Here the blue color indicates the bands derived
from $p_{z}$ orbitals while the red color denotes the bands derived from $s$,
$p_{x}$, $p_{y}$ orbitals. Since $p_{z}$ orbitals are odd whereas $s$,
$p_{x}$, $p_{y}$ orbitals are even under $M_{z}$ symmetry, the energy spectrum
from $p_{z}$ orbitals is not hybridized with that from $s$, $p_{x}$, $p_{y}$
orbitals. One can see that the low-energy band structure has approximate
chiral symmetry and can be effectively described by using only $p_{z}$
orbitals. However, to capture the higher-order band topology, the core
electronic states derived from $s$, $p_{x}$, $p_{y}$ orbitals should be
included as shown below.

In general, a 2D $P$-symmetric spinless fermion system carries three $Z_{2}$
topological invariants, $w_{1x}$, $w_{1y}$, and $w_{2}$ ahn2018band . In terms
of the sewing matrix $G_{mn}(\textbf{k})=\langle
u_{m\textbf{-k}}|P|u_{n\textbf{k}}\rangle$, $w_{1a}$ ($a=x,~{}y$) is given by
its 1D winding number as

$\displaystyle w_{1a}=-\frac{i}{2}\oint_{C}\nabla_{\textbf{k}_{a}}\log\det
G(\textbf{k})\mod 2,$ (1)

which is equivalent to the quantized Berry phase $\Phi_{a}$, in a way that
$\Phi_{a}=\pi w_{1a}$ (mod $2\pi$) zak1989berry ; xiao2010berry . $w_{1a}$ can
be determined by using relation
$(-1)^{w_{1a}}=\prod_{i=1}^{2}\xi(\Gamma_{ia})$ where $\xi(\Gamma_{ia})$
($i=1,~{}2$) is the product of parity eigenvalues of occupied states at the
time-reversal invariant momentum (TRIM) $\Gamma_{ia}$, and $\Gamma_{1a}$,
$\Gamma_{2a}$ are two TRIMs along the reciprocal lattice vector $G_{a}$. Since
$w_{1a}$ is equivalent to the quantized charge polarization along the $G_{a}$
direction, it can be considered as a 1D topological invariant.

On the other hand, $w_{2}$ is a 2D topological invariant given by the 2D
winding number of the sewing matrix $G(\textbf{k})$ (modulo 2) ahn2018higher ,
and measures the higher-order band topology of $P$-symmetric spinless fermion
systems ahn2018higher ; ahn2019failure ; wieder2018axion . $w_{2}$ can be
determined by

$\displaystyle(-1)^{w_{2}}=\prod_{i=1}^{4}(-1)^{[N_{\text{occ}}^{-}(\Gamma_{i})/2]},$
(2)

where $N_{\text{occ}}^{-}(\Gamma_{i})$ is the number of occupied bands with
odd parity at the TRIM $\Gamma_{i}$ and the square bracket $[\alpha]$
indicates the greatest integer value of $\alpha$ song2018diagnosis ;
fu2007topological ; kim2015dirac ; hughes2011inversion ; turner2012quantized ;
po2017symmetry .

Using the band structure from first-principles calculations and the
corresponding parity eigenvalues at TRIMs, one can easily show that MGD is a
higher-order topological insulator with $w_{1x}=w_{1y}=0$ and $w_{2}=1$. One
interesting feature of the MGD band structure is that we get $w_{2}=0$ when
only the bands derived from $p_{z}$ orbitals are considered. Namely, although
the low-energy band structure itself can be well-described by using only
$p_{z}$ orbitals, the higher-order topology can be correctly captured only
when the core electronic states derived from $s$, $p_{x}$, $p_{y}$ orbitals
are included. This result can be confirmed by tight-binding analysis as well
as first-principles calculations. Since the bands derived from $p_{z}$
orbitals do not hybridize with those derived from other orbitals because of
$M_{z}$ symmetry, one can compute the relevant topological invariants
separately.

The orbital dependence of $w_{2}$ can also be confirmed by Wilson loop
calculations ahn2018band ; ahn2019failure ; zhao2017p ; fang2015topological .
The Wilson loop is defined as a path ordered product of the exponential of
Berry connections yu2011equivalent ; alexandradinata2014wilson ,

$\displaystyle W_{(k_{1}+2\pi,k_{2})\leftarrow(k_{1},k_{2})}$
$\displaystyle=\lim_{N\rightarrow\infty}F_{N-1}F_{N-2}\dotsb F_{1}F_{0},$ (3)

where $[F_{i}]_{mn}=\langle u_{m}(\frac{2\pi(i+1)}{N},k_{2})|u_{n}(\frac{2\pi
i}{N},k_{2})\rangle$, and its spectrum is gauge-invariant yu2011equivalent .
Computed along the $k_{1}$ direction parallel to the reciprocal lattice vector
$\bf{G_{1}}$ from ($k_{1},k_{2}$), the set of Wilson loop eigenvalues
$\\{e^{i\nu_{i}(k_{2})}\\}$ indicates the position of Wannier centers at given
$k_{2}$, and the corresponding total charge polarization is given by
$p_{1}(k_{2})=\frac{1}{2\pi}\sum_{i=1}^{N_{occ}}\nu_{i}(k_{2})$. The Wilson
loop spectrum of MGD shows that the charge polarization is zero in both
$k_{1}$ and $k_{2}$ directions, which is consistent with the parity eigenvalue
analysis (Fig. 1(c)). Also, the single linear-crossing at
$(k_{y},\nu)=(0,\pi)$ indicates $w_{2}=1$ ahn2018band . Projecting the Wilson
loop operator into the $p_{z}$ ($s$, $p_{x}$, and $p_{y}$) orbital basis, we
observe $w_{2}=0$ ($w_{2}=1$), which confirms that the higher-order band
topology originates from the core electronic levels. The orbital dependence of
$w_{2}$ can be further confirmed by computing the Wannier center of each
orbital explicitly and also by computing the nested Wilson loop that measures
the higher-order band topology wieder2018axion ; benalcazar2017quantized . In
the presence of $PT$ symmetry, the Wilson loop spectrum can be divided into
two sub-bands that are separated by a gap, each centered at $\nu=0$ or
$\nu=\pi$, respectively. When the Wilson loop operator at given $k_{2}$ is
diagonalized as

$\displaystyle
W_{(k_{1}+2\pi,k_{2})\leftarrow(k_{1},k_{2})}=\sum_{i}|\nu_{i}(k_{2})\rangle
e^{i\nu_{i}(k_{2})}\langle\nu_{i}(k_{2})|,$ (4)

then, using the set of eigenvectors corresponding to the sub-band centered at
$\nu=\pi$, the nested Wilson loop operator $\tilde{W}$ along the $\bf{G_{2}}$
vector is defined as

$\displaystyle\tilde{W}_{(k_{2}+2\pi)\leftarrow(k_{2})}=\lim_{N\rightarrow\infty}\tilde{F}_{N-1}\tilde{F}_{N-2}\dotsb\tilde{F}_{0},$
(5)

where
$[\tilde{F}_{i}]_{mn}=\langle\nu_{m}(\frac{2\pi(i+1)}{N})|\nu_{n}(\frac{2\pi
i}{N})\rangle$. Since the determinant of the nested Wilson loop operator
$\det(\tilde{W})$ is identical to $(-1)^{w_{2}}$ ahn2018band , we find that
$\det(\tilde{W})$ is nontrivial. Projecting the Wilson loop operator into the
$p_{z}$ orbital basis, we observe that the corresponding nested Wilson loop
spectrum is trivial. On the other hand, within the $s$, $p_{x}$, and $p_{y}$
orbital basis, the nested Wilson loop spectrum is nontrivial, which confirms
that the higher-order band topology originates from the core electronic
levels.

###  II.2 Filling anomaly and corner charges

The higher-order band topology of MGD with $w_{2}=1$ induces a pair of
anomalous corner states wang2019higher . To see this, we consider a finite-
size structure invariant under $P$ symmetry. When the system is chiral
symmetric, there are two robust zero-energy in-gap states related by $P$ whose
locations are fixed at $M_{y}$-invariant corners. One way to see this is to
use the topological classification proposed by Geier et al. geier2018second .
When a 2D system has commuting mirror and chiral symmetries, it belongs to
BDI${}^{M_{++}}$ class that supports nontrivial topology, and the
corresponding insulator has zero modes at mirror invariant corners. In MGD,
the $M_{y}$ symmetry commutes with chiral symmetry since $M_{y}$ symmetry does
not mix different sublattices. On the other hand, when chiral symmetry and
$M_{x}$ symmetry anti-commute, the system belongs to BDI${}^{M_{+-}}$ class,
which is topologically trivial, so there are no protected corner charges at
$M_{x}$ invariant corners.

To be specific, we can find the location of corner states by studying the low-
energy Hamiltonian of MGD near the $\Gamma$ point:

$\displaystyle H(k_{x},k_{y})$ $\displaystyle=\hbar
v(k_{x}\Gamma_{1}+k_{y}\Gamma_{2})+\Delta\Gamma_{3},$ (6)

where
$\Gamma_{1}=\sigma_{x}\tau_{x},\Gamma_{2}=\sigma_{y}\tau_{x},\Gamma_{3}=\tau_{z},$
and $\Gamma_{4}=\sigma_{z}\tau_{x}$. Corresponding symmetry representations
are $T=\sigma_{x}\tau_{z}K$, $P=\tau_{z}$, $M_{x}=\sigma_{x}\tau_{z}$,
$M_{y}=\sigma_{x}$, and chiral symmetry $S$ is described by
$S=\Gamma_{5}=\tau_{y}$. Considering a disk-shaped geometry, the low-energy
Hamiltonian in real space becomes
$H=-i\tilde{\Gamma}_{1}(\theta)\partial_{r}-i\tilde{\Gamma}_{2}(\theta)r^{-1}\partial_{\theta}+\Delta(r)\Gamma_{3}$,
where $\tilde{\Gamma}_{1}(\theta)=\cos\theta\Gamma_{1}+\sin\theta\Gamma_{2}$
and $\tilde{\Gamma}_{2}(\theta)=-\sin\theta\Gamma_{1}+\cos\theta\Gamma_{2}$.
We assume that $\Delta(r)>0$ inside the bulk and $\Delta(r)<0$ outside the
bulk. The edge Hamiltonian is obtained by applying the projection operator
$P(\theta)=\frac{1}{2}(1+i\tilde{\Gamma}_{1}(\theta)\Gamma_{3})$, leading to
$H_{\text{edge}}=-iR^{-1}\tilde{\sigma}_{z}\partial_{\theta}$ where
$\tilde{\sigma}_{z}=P\tilde{\Gamma}_{2}(\theta)P^{-1}$. On the edge, position-
dependent mass terms that open the edge gap are allowed. Mass terms that are
commuting with $P(\theta)$ and anti-commuting with
$\tilde{\Gamma}_{2}(\theta)$ and $S$ are
$H_{m}=m_{4}(R,\theta)\Gamma_{4}+im_{25}(R,\theta)\tilde{\Gamma}_{2}\Gamma_{5}$.
The corresponding edge Hamiltonian is
$H_{\text{edge}}=-iR^{-1}\tilde{\sigma}_{z}\partial_{\theta}+\tilde{m}_{1}(R,\theta)\tilde{\sigma}_{x}$,
where $\tilde{m}_{1}=m_{4}+m_{25}$,
$\tilde{\sigma}_{x}=P\Gamma_{4}P^{-1}=Pi\tilde{\Gamma}_{2}\Gamma_{5}P^{-1}$.
Here $PT$ symmetry requires $\tilde{m}_{1}(\theta+\pi)=-\tilde{m}_{1}(\theta)$
while $M_{y}$ requires $\tilde{m}_{1}(x,y)=-\tilde{m}_{1}(x,-y)$. Thus, the
mass gap is closed at $M_{y}$-invariant corners leading to localized corner
charges.

Due to the degeneracy of the two in-gap states, the half-filling condition
cannot be satisfied as long as $P$ is preserved, which is known as the filling
anomaly benalcazar2019quantization ; wieder2018axion ; song2017d . When chiral
symmetry is broken, the two degenerate in-gap states can be shifted from zero
energy and even merged to the valence (conduction) bands generating a single
hole (electron) at half-filling. Namely, the presence of an odd number of
holes (electrons) in the valence (conduction) band is the manifestation of the
higher-order band topology of generic 2D HOTIs lacking chiral symmetry
wieder2018axion ; cualuguaru2019higher ; okuma2019topological ;
hwang2019fragile . To resolve the filling anomaly while keeping $P$, one needs
to add an extra electron (or hole), which leads to the accumulation or
depletion of a half-integral charge at each $M_{y}$-invariant corner.

Explicitly, let us illustrate how the number of states can be counted in a
finite-size HOTI with $w_{2}=1$ following Ref. wieder2018axion, . When the
total number of states is $N$ and chiral symmetry exists, there are $N/2-1$
occupied and unoccupied states, respectively, since two in-gap states appear
at zero-energy. When chiral symmetry is broken while inversion is preserved,
the two in-gap states related by inversion are absorbed into either the
valence or the conduction bands simultaneously. Then the number of states
below the gap is $n_{\text{gap}}=N/2\pm 1$. Therefore
$n_{\text{gap}}=(\frac{N}{2}+$ an odd integer) when $w_{2}=1$, whereas
$n_{\text{gap}}=(\frac{N}{2}+$ an even integer) when $w_{2}=0$. When the
valence (conduction) band of a system with $w_{2}=1$ is completely occupied
(unoccupied) by adding electrons or holes additionally, a half-integral
electric charge should be accumulated (or depleted) at two corners related by
inversion.

For the finite-size MGD composed of $n\times n$ unit cells shown in Fig. 2(a),
the total number of states are $N=72n^{2}$ when $2s$, $2p_{x}$, $2p_{y}$, and
$2p_{z}$ orbitals are considered. In MGD, there is a single $sp$ electron
localized at every edge of a hexagonal unit cell. If an edge of a unit cell is
not shared with an adjacent unit cell, the relevant $sp$ orbital becomes a
half-filled zero-energy non-bonding state when chiral symmetry exists. When we
consider the geometry of a finite-size structure shown in Fig. 2(a), MGD has
$8n-2$ non-bonding states along the boundary. When chiral symmetry is broken,
$8n-2$ states can be merged into the valence band, so that
$n_{\text{gap}}=36n^{2}+4n-1=\frac{N}{2}+4n-1$. This is what is observed in
first-principles calculations and the corresponding density of states is shown
in Fig. 3. Thus, when all the states below the gap are occupied, and thus the
filling anomaly is lifted, a half-integral charge (modulo an integer charge)
is accumulated at two $M_{y}$-invariant corners. This can be easily understood
because an odd number of non-bonding states exists only at the
$M_{y}$-invariant corners. [See Fig. 2(a).]

To observe the corner states more clearly, we add a hydrogen atom to the
carbon atom at every corner, except at two $M_{y}$-invariant corners, and keep
$P$ and $M_{y}$ symmetries. Then the number of added hydrogen atoms is an
integer multiple of four, which does not alter the band topology of MGD. On
the other hand, if two hydrogen atoms are added at the $M_{y}$-invariant
corners, the total number of bands becomes $N^{\prime}=N+2$, while
$n_{\text{gap}}$ remains the same as described in Fig. 3(a). Then
$n_{\text{gap}}=\frac{N^{\prime}}{2}+4n-2$ and $w_{2}$ becomes trivial.
Therefore, the number of added hydrogen atoms should be an integer multiple of
four to keep $w_{2}$ unchanged. Maximally, $(8n-4)$ hydrogen atoms can be
added to the finite-size MGD with two non-bonding states remaining at
$M_{y}$-invariant corners. The corresponding density of states is shown in
Fig. 3(c). Here a half electric charge is accumulated at each
$M_{y}$-invariant corner when the states below the gap are fully occupied as
shown in Fig. 2(b).

###  II.3 Wannier function description

Here we show that the higher-order band topology is associated with the
chemical bonding of the $s$, $p_{x}$, and $p_{y}$ orbitals across the unit
cell boundary. In other words, Wannier functions are located away from the
atomic sites and cannot be moved to the atomic sites by a continuous
deformation that preserves inversion symmetry. Here we follow the idea
proposed by van Miert and Ortix van2018higher  to analyze the Wannier centers
in inversion-symmetric systems.

Let us take the inversion-invariant unit cell shown in Fig. 1(a) whose center
is located at the position ${\bf m}$. There are four Wyckoff positions in a
unit cell labelled as $A,B,C,$ and $D$ —located at ${\bf 0}$, $\frac{1}{2}{\bf
a_{1}}$, $\frac{1}{2}{\bf a_{2}}$, and $\frac{1}{2}{\bf a_{1}}+\frac{1}{2}{\bf
a_{2}}$ from the center of the unit cell, respectively — which are invariant
under inversion up to lattice vectors ${\bf a}_{1}$ and ${\bf a}_{2}$. The
symmetric Wannier functions centered at Wyckoff positions transform under
inversion as follows:

$\displaystyle w_{A\pm,{\bf m}}({\bf-x})$ $\displaystyle=\pm
w_{A\pm,{\bf-m}}({\bf x}),$ $\displaystyle w_{B\pm,{\bf m}}({\bf-x})$
$\displaystyle=\pm w_{B\pm,{\bf-m-a_{1}}}({\bf x}),$ $\displaystyle
w_{C\pm,{\bf m}}({\bf-x})$ $\displaystyle=\pm w_{C\pm,{\bf-m-
a_{1}-a_{2}}}({\bf x}),$ $\displaystyle w_{D\pm,{\bf m}}({\bf-x})$
$\displaystyle=\pm w_{D\pm,{\bf-m-a_{2}}}({\bf x}),$ (7)

where $w_{\alpha\pm,{\bf m}}({\bf x})$ denotes a Wannier function with the
coordinate ${\bf x}$, centered at the Wyckoff position $\alpha$ of the unit
cell with the center at ${\bf m}$. The sign $\pm$ indicates eigenvalues of the
inversion operation with respect to the Wyckoff position. We define the number
of symmetric Wannier functions centered at each Wyckoff position with
inversion eigenvalues $\pm 1$ as $N_{A\pm},N_{B\pm},N_{C\pm}$, and $N_{D\pm}$.
Let us note that two Wannier functions with inversion eigenvalues $+1$ and
$-1$ at $A$ can be expressed as a linear combination of two Wannier functions
with inversion eigenvalues $+1$ and $-1$ at any other Wyckoff position. Thus,
only the difference $\nu_{\alpha}=N_{\alpha+}-N_{\alpha-}$ at each Wyckoff
position $\alpha$ is a well-defined $\mathbb{Z}$ invariant van2018higher .

$\displaystyle\nu_{A}=N_{A+}-N_{A-},$ $\displaystyle\nu_{B}=N_{B+}-N_{B-},$
$\displaystyle\nu_{C}=N_{C+}-N_{C-},$ $\displaystyle\nu_{D}=N_{D+}-N_{D-},$
(8)

which are related to the inversion eigenvalues at TRIMs in the Brillouin zone
as van2018higher

$\displaystyle\nu_{A}$
$\displaystyle=-\Gamma_{1}-\frac{1}{2}\Gamma_{-1}+\frac{1}{2}M_{-1}+\frac{1}{2}X_{-1}+\frac{1}{2}Y_{-1}=-1,$
$\displaystyle\nu_{B}$
$\displaystyle=\frac{1}{2}\Gamma_{-1}-\frac{1}{2}M_{-1}+\frac{1}{2}X_{-1}-\frac{1}{2}Y_{-1}=-1,$
$\displaystyle\nu_{C}$
$\displaystyle=\frac{1}{2}\Gamma_{-1}-\frac{1}{2}M_{-1}-\frac{1}{2}X_{-1}+\frac{1}{2}Y_{-1}=-1,$
$\displaystyle\nu_{D}$
$\displaystyle=\frac{1}{2}\Gamma_{-1}+\frac{1}{2}M_{-1}+\frac{1}{2}X_{-1}-\frac{1}{2}Y_{-1}=-1.$
(9)

In monolayer graphdiyne (MGD), we find that
$\nu_{A}=\nu_{B}=\nu_{C}=\nu_{D}=-1$ from the inversion eigenvalues at TRIM.
In the viewpoint of chemistry, $|\nu_{A}|=1$ comes from the electrons on the
benzene ring around the unit cell center, $|\nu_{B}|=|\nu_{C}|=|\nu_{D}|=1$
originates from the $sp$ bonding across the unit cell boundary.

To relate these indices with the two-dimensional topological invariant $w_{2}$
defined before, we define electric dipole moments $p_{i=1,2}$ and an electric
quadrupole moment $q_{12}$ for the unit cell as follows,

$\displaystyle p_{x}$ $\displaystyle\equiv\sum_{\begin{subarray}{c}i\in\rm
occ,\\\ \rm unit\;cell\end{subarray}}X_{i}=\frac{1}{2}(N_{B}+N_{C})=0\mod 1,$
$\displaystyle p_{y}$ $\displaystyle\equiv\sum_{\begin{subarray}{c}i\in\rm
occ,\\\ \rm unit\;cell\end{subarray}}Y_{i}=\frac{1}{2}(N_{C}+N_{D})=0\mod 1,$
$\displaystyle q_{xy}$ $\displaystyle\equiv\sum_{\begin{subarray}{c}i\in\rm
occ,\\\ \rm
unit\;cell\end{subarray}}X_{i}Y_{i}=\frac{1}{4}N_{C}=\frac{1}{4}\mod\frac{1}{2},$
(10)

where $X_{i}$ and $Y_{i}$ are the ${\bf a}_{1}$ and ${\bf a}_{2}$ components
of the Wannier center for the $i$th occupied Wannier function. $p_{x}$ and
$p_{y}$ are nothing but the polarizations along ${\bf a}_{1}$ and ${\bf
a}_{2}$, respectively. $w_{2}$ is related to $p_{x,y}$, $q_{xy}$ via
$w_{2}=4(p_{x}p_{y}-q_{xy})$ mod 2 ahn2018band ; ahn2019failure , so

$\displaystyle w_{2}$ $\displaystyle=(N_{B}+N_{C})(N_{C}+N_{D})-N_{C}$
$\displaystyle=N_{B}N_{C}+N_{C}N_{D}+N_{D}N_{B}\mod 2$
$\displaystyle=\nu_{B}\nu_{C}+\nu_{C}\nu_{D}+\nu_{D}\nu_{B}\mod 2,$ (11)

because $\nu=N_{+}-N_{-}=N_{+}+N_{-}-2N_{-}=N-2N_{-}=N$ mod 2. This shows that
$w_{2}=1\mod 2$ when $|\nu_{B}|=|\nu_{C}|=|\nu_{D}|=1$ due to the $sp$ bonding
across the unit cell boundary.

##  III Discussion

Since MGD has $w_{2}=1$, when layers of MGD are stacked vertically, the
resulting 3D insulator with $PT$ symmetry becomes a 3D weak topological
insulator ahn2018band , dubbed a 3D weak Stiefel-Whitney insulator, when
inter-layer coupling is weak. When inter-layer coupling becomes large enough,
an accidental band crossing can happen at a TRIM at which a pair of nodal
lines with $Z_{2}$ monopole charge are created. The existence of two nodal
lines at $k_{z}=\pm k_{c}$ $(k_{c}>0)$ was predicted in ABC-stacked graphdiyne
nomura2018three , and their $Z_{2}$ monopole charge was also confirmed
ahn2018band  based on first-principles calculations. In ABC-stacked
graphdiyne, the 2D subspace with a fixed $k_{z}$ carries $w_{2}=1$ ($w_{2}=0$)
when $|k_{z}|<k_{c}$ ($k_{c}<|k_{z}|<\pi$) since the band inversion is
happened at $\bm{k}=(0,0,\pi)$. Since a 2D HOTI with $w_{2}=1$ possesses
corner charges, similar corner charges are also expected in ABC-stacked
graphdiyne in the subspace with a fixed $k_{z}$ where $w_{2}=1$ wang2019higher
, which leads to hinge modes of the 3D structure shown in Fig. 4(a).

To demonstrate the hinge modes of ABC-stacked graphdiyne, we study a tight-
binding model by using only $p_{z}$ orbitals to reduce numerical costs. The
energy spectrum of a finite-size system with $PT$ and $C_{2x}$ symmetries is
shown in Fig. 4(b), which clearly shows the presence of hinge modes in the
momentum region where $w_{2}=1$ as expected. When only $p_{z}$ orbitals are
included, although the low-energy band structure from the tight-binding model
is consistent with that from first-principles calculations, the topological
property is different. Namely, the tight-binding model predicts $w_{2}=0$
($w_{2}=1$) when $|k_{z}|<k_{c}$ ($k_{c}<|k_{z}|<\pi$), which is opposite to
the result from the first-principles calculations. Such a discrepancy can be
remedied when the core electronic levels are included in the tight-binding
calculation, which is confirmed by separate calculations.

To sum up, we have shown that the nontrivial band topology of a 2D HOTI
lacking chiral symmetry can be revealed by using the filling anomaly, which
can be explicitly demonstrated by counting the number of states of a finite-
size system. This idea can be straightforwardly generalized to any $d$-th
order topological insulator in $d$-dimensions hosting zero-dimensional corner
states. We believe that our theoretical study provides a general way to extend
the scope of HOTI materials to a wider class.

##  IV Methods

###  IV.1 DFT calculations

We have carried out the first-principles density-functional theory
calculations using the Vienna ab initio simulation package kresse1996efficient
. The projector augmented-wave method blochl1994projector  and the exchange-
correlation functional of the generalized gradient approximation in the
Perdew, Burke, and Ernzerhof perdew1996generalized  scheme were used. The
self-consistent total energy was evaluated with an $12\times 12\times 1$
$k$-point mesh, and the cutoff energy for the plane-wave basis set was 500 eV.
For the finite-size structure, we used only $\Gamma$ point. A single unit cell
contains 18 carbon atoms, where a benzene ring is connected to six neighboring
ones by linear chains. In the optimized structure, we have $L=9.46$ Å,
$b_{1}=1.43$ Å, $b_{2}=1.40$ Å, $b_{3}=1.23$ Å, and $b_{4}=1.34$ Å.

###  IV.2 Tight-binding calculation

Low-energy band structure of ABC-stacked graphdiyne can be properly described
by $p_{z}$ orbital basis. The intralayer coupling is described by the nearest
neighbor hopping parameters expressed as the transfer integral
$V_{pp\pi}=V_{pp\pi}^{0}e^{(-R-a_{0})/r_{0}}$, where
$V_{pp\pi}^{0}\approx-2.7eV$ is the hopping parameter for graphene’s nearest
neighbor interaction ($a_{0}=0.142nm$ and $r_{0}=0.09nm$) nomura2018three .
The interlayer coupling is described by two hopping parameters of two shortest
vertical bonds which are expressed as the transfer integral
$V_{pp\sigma}=V_{pp\sigma}^{0}e^{(-R-d_{0})/r_{0}}$, where
$V_{pp\sigma}\approx 0.48eV$ is the hopping parameter at the interlayer
distance of graphite ($d_{0}=0.335nm$) nomura2018three .

##  V Notes

During the preparation of our manuscript, we have found a related work
sheng2019two . In [sheng2019two, ], the authors also observed HOTI phase in
MGD. However, some of their theoretical arguments are inconsistent with our
results. For instance, in [sheng2019two, ], it is stated that crystalline
symmetries are not essential for the protection of corner states, which is in
sharp contrast to our paper where the role of $P$ symmetry is emphasized.
Also, in [sheng2019two, ], it is mentioned that double band inversion of
$p_{z}$ orbitals at the $\Gamma$ point is the origin of the HOTI phase. This
is different from our result in that higher-order band topology originates
from the core levels derived from $s,p_{x},p_{y}$ orbitals.

##  VI Data Availability

The datasets generated during and/or analyzed during the current study are
available from the corresponding author on reasonable request.

##  VII Acknowledgements

We acknowledge the helpful discussions with Dongwook Kim and Youngkuk Kim.

##  VIII Competing Interests

The Authors declare no competing interests.

##  IX Author Contributions

E.L. performed tight-binding calculations, R.K. performed first-principles
calculations, and E.L. and J.A. did theoretical analysis. B.-J.Y. supervised
the project. All authors analyzed the data and wrote the manuscript.
∗These authors contributed equally.

##  X Funding

E. L., R. K. and J.A. were supported by IBS-R009-D1. B.-J.Y. was supported by
the Institute for Basic Science in Korea (Grant No. IBS-R009-D1) and Basic
Science Research Program through the National Research Foundation of Korea
(NRF) (Grant No.0426-20180011), and the POSCO Science Fellowship of POSCO TJ
Park Foundation (No.0426-20180002). This work was supported in part by the
U.S. Army Research Office under Grant Number W911NF-18-1-0137.

## References

  * (1) Schindler, F. _et al._ Higher-order topology in bismuth.  _Nature Physics_ 14, 918 (2018).
  * (2) Schindler, F. _et al._ Higher-order topological insulators.  _Science Advances_ 4, eaat0346 (2018).
  * (3) Wang, Z., Wieder, B. J., Li, J., Yan, B. & Bernevig, B. A.  Higher-order topology, monopole nodal lines, and the origin of large fermi arcs in transition metal dichalcogenides XTe2 (X= Mo, W).  _Physical Review Letters_ 123, 186401 (2019).
  * (4) Yue, C. _et al._ Symmetry-enforced chiral hinge states and surface quantum anomalous hall effect in the magnetic axion insulator Bi2-xSmxSe3.  _Nature Physics_ 1 (2019).
  * (5) Wieder, B. J. & Bernevig, B. A.  The axion insulator as a pump of fragile topology.  _Preprint at http://arXiv.org/abs/1810.02373_ (2018).
  * (6) Ahn, J. & Yang, B.-J.  Symmetry representation approach to topological invariants in $c_{2z}t$-symmetric systems.  _Physical Review B_ 99, 235125 (2019).
  * (7) Benalcazar, W. A., Bernevig, B. A. & Hughes, T. L.  Quantized electric multipole insulators.  _Science_ 357, 61–66 (2017).
  * (8) Langbehn, J., Peng, Y., Trifunovic, L., von Oppen, F. & Brouwer, P. W.  Reflection-symmetric second-order topological insulators and superconductors.  _Physical Review Letters_ 119, 246401 (2017).
  * (9) Geier, M., Trifunovic, L., Hoskam, M. & Brouwer, P. W.  Second-order topological insulators and superconductors with an order-two crystalline symmetry.  _Physical Review B_ 97, 205135 (2018).
  * (10) Khalaf, E.  Higher-order topological insulators and superconductors protected by inversion symmetry.  _Physical Review B_ 97, 205136 (2018).
  * (11) Benalcazar, W. A., Li, T. & Hughes, T. L.  Quantization of fractional corner charge in $c_{n}$-symmetric higher-order topological crystalline insulators.  _Physical Review B_ 99, 245151 (2019).
  * (12) Song, Z., Fang, Z. & Fang, C.  (d- 2)-dimensional edge states of rotation symmetry protected topological states.  _Physical review letters_ 119, 246402 (2017).
  * (13) Ahn, J., Park, S. & Yang, B.-J.  Failure of nielsen-ninomiya theorem and fragile topology in two-dimensional systems with space-time inversion symmetry: application to twisted bilayer graphene at magic angle.  _Physical Review X_ 9, 021013 (2019).
  * (14) Lee, E., Furusaki, A. & Yang, B.-J.  Fractional charge bound to a vortex in two dimensional topological crystalline insulators.  _Preprint at http://arXiv.org/abs/1903.02737_ (2019).
  * (15) Ezawa, M.  Minimal models for wannier-type higher-order topological insulators and phosphorene.  _Physical Review B_ 98, 045125 (2018).
  * (16) Călugăru, D., Juričić, V. & Roy, B.  Higher-order topological phases: A general principle of construction.  _Physical Review B_ 99, 041301 (2019).
  * (17) Okuma, N., Sato, M. & Shiozaki, K.  Topological classification under nonmagnetic and magnetic point group symmetry: Application of real-space atiyah-hirzebruch spectral sequence to higher-order topology.  _Physical Review B_ 99, 085127 (2019).
  * (18) Hwang, Y., Ahn, J. & Yang, B.-J.  Fragile topology protected by inversion symmetry: Diagnosis, bulk-boundary correspondence, and wilson loop.  _Physical Review B_ 100, 205126 (2019).
  * (19) Nomura, T., Habe, T., Sakamoto, R. & Koshino, M.  Three-dimensional graphdiyne as a topological nodal-line semimetal.  _Physical Review Materials_ 2, 054204 (2018).
  * (20) Ahn, J., Kim, D., Kim, Y. & Yang, B.-J.  Band topology and linking structure of nodal line semimetals with Z2 monopole charges.  _Physical Review Letters_ 121, 106403 (2018).
  * (21) Haley, M. M., Brand, S. C. & Pak, J. J.  Carbon networks based on dehydrobenzoannulenes: synthesis of graphdiyne substructures.  _Angewandte Chemie International Edition in English_ 36, 836–838 (1997).
  * (22) Li, G. _et al._ Architecture of graphdiyne nanoscale films.  _Chemical Communications_ 46, 3256–3258 (2010).
  * (23) Zak, J.  Berry’s phase for energy bands in solids.  _Physical Review Letters_ 62, 2747 (1989).
  * (24) Xiao, D., Chang, M.-C. & Niu, Q.  Berry phase effects on electronic properties.  _Reviews of Modern Physics_ 82, 1959 (2010).
  * (25) Song, Z., Zhang, T. & Fang, C.  Diagnosis for nonmagnetic topological semimetals in the absence of spin-orbital coupling.  _Physical Review X_ 8, 031069 (2018).
  * (26) Fu, L. & Kane, C. L.  Topological insulators with inversion symmetry.  _Physical Review B_ 76, 045302 (2007).
  * (27) Kim, Y., Wieder, B. J., Kane, C. & Rappe, A. M.  Dirac line nodes in inversion-symmetric crystals.  _Physical Review Letters_ 115, 036806 (2015).
  * (28) Hughes, T. L., Prodan, E. & Bernevig, B. A.  Inversion-symmetric topological insulators.  _Physical Review B_ 83, 245132 (2011).
  * (29) Turner, A. M., Zhang, Y., Mong, R. S. & Vishwanath, A.  Quantized response and topology of magnetic insulators with inversion symmetry.  _Physical Review B_ 85, 165120 (2012).
  * (30) Po, H. C., Vishwanath, A. & Watanabe, H.  Symmetry-based indicators of band topology in the 230 space groups.  _Nature Communications_ 8, 50 (2017).
  * (31) Zhao, Y. & Lu, Y.  $pt$-symmetric real dirac fermions and semimetals.  _Physical Review Letters_ 118, 056401 (2017).
  * (32) Fang, C., Chen, Y., Kee, H.-Y. & Fu, L.  Topological nodal line semimetals with and without spin-orbital coupling.  _Physical Review B_ 92, 081201 (2015).
  * (33) Yu, R., Qi, X. L., Bernevig, A., Fang, Z. & Dai, X.  Equivalent expression of Z2 topological invariant for band insulators using the non-abelian berry connection.  _Physical Review B_ 84, 075119 (2011).
  * (34) Alexandradinata, A., Dai, X. & Bernevig, B. A.  Wilson-loop characterization of inversion-symmetric topological insulators.  _Physical Review B_ 89, 155114 (2014).
  * (35) Van Miert, G. & Ortix, C.  Higher-order topological insulators protected by inversion and rotoinversion symmetries.  _Physical Review B_ 98, 081110 (2018).
  * (36) Kresse, G. & Furthmüller, J.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  _Physical Review B_ 54, 11169 (1996).
  * (37) Blöchl, P. E.  Projector augmented-wave method.  _Physical Review B_ 50, 17953 (1994).
  * (38) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Physical Review Letters_ 77, 3865 (1996).
  * (39) Sheng, X.-L. _et al._ Two-dimensional second-order topological insulator in graphdiyne.  _Preprint at http://arXiv.org/abs/1904.09985_ (2019).

Figure 1:  Lattice structure, band structure, and the Wilson loop spectrum.
(a) A schematic figure describing the unit cell and the relevant atomic
orbitals of a monolayer graphdiyne (MGD) composed of 18 carbon atoms. Two
sublattices are marked with yellow and blue colors, respectively. Since MGD
has a bipartite lattice structure, chiral symmetry exists when only the
nearest neighbor hoppings are considered. (b) The band structure of MGD
obtained by first-principles calculations, which shows approximate chiral
symmetry near the Fermi level. Since $p_{z}$ orbitals are odd while $s$,
$p_{x}$, $p_{y}$ orbitals are even under the mirror $M_{z}:z\rightarrow-z$
operation, the energy spectrum from $p_{z}$ orbitals (blue) is not hybridized
with that from $s$, $p_{x}$, $p_{y}$ orbitals (red). (c) The Wilson loop
spectrum computed including $s$, $p_{x}$, $p_{y}$, $p_{z}$ orbitals. The
spectrum exhibits a crossing point at $k_{y}=0$, $\theta=\pi$, which indicates
$w_{2}=1$.  Figure 2:  Geometry of a finite-size structure and corner charges.
(a) A finite-size structure of MGD preserving $P$ and $M_{y}$ symmetries
designed to observe the higher-order band topology protected by $P$ symmetry.
The red horizontal line indicates the $M_{y}$-invariant line. (b) Electron
density distribution for a finite-size MGD composed of $9\times 9$ unit cells
where hydrogen atoms are attached at every corner except at the two
$M_{y}$-invariant corners. To resolve the filling anomaly, we add one electron
to MGD additionally and fill the valence band completely. The accumulated or
depleted charges with respect to the half-filled configuration are plotted.
Here a half-integral charge accumulated at each $M_{y}$-invariant corner
appears due to the nontrivial 2D topological invariant $w_{2}=1$.  Figure 3:
Filling anomaly and density of states. (a) A schematic figure displaying the
energy spectrum of a finite-size MGD composed of $n\times n$ unit cells, where
the total number of states is $N=72n^{2}$. Here a thick gray strip near the
Fermi level $E_{F}$ indicates the non-bonding states arising from carbon atoms
along the boundary. The upper panel corresponds to the case at half-filling
where an odd number of holes is below the gap, which demonstrates the higher-
order band topology with $w_{2}=1$. The lower panel shows the case when two
hydrogen atoms are added to the finite-size MGD. Here the total number of
states $N^{\prime}$ increases by two ($N^{\prime}=N+2$). The hybridization
between two hydrogen states and two non-bonding states generates two bonding
and two anti-bonding states, so that the number of states below the gap
changes from $\frac{N}{2}+4n-1$ to $\frac{N^{\prime}}{2}+4n-2$. Then, $w_{2}$
also changes from 1 to 0. Thus, to maintain the value of $w_{2}$, the number
of hydrogen atoms attached for passivation should be an integer multiple of
four. (b) Density of states (DOS) of a finite-size MGD without hydrogen
passivation. Here the carbon atom at each corner has a non-bonding state.
There are $(4n-1)$ holes below the gap at half-filling. The green color
indicate the gapped region. (c) DOS of a finite-size MGD with hydrogen
passivation where $(8n-4)$ hydrogen atoms are attached along the boundary.
Here only two carbon atoms at $M_{y}$-invariant corners have non-bonding
states. There is a single hole below the gap at half-filling. Since both
systems in (b) and (c) have an odd number of holes, both are HOTIs with
$w_{2}=1$.  Figure 4:  Hinge modes of ABC-stacked graphdiyne with monopole
nodal lines. (a) A schematic figure describing the 3D geometry of ABC-stacked
graphdiyne which is finite in the $xy$ plane but periodic in the vertical
direction. The structure must preserve $C_{2x}$ symmetry to carry hinge
states. (b) Hinge modes, located in the region $k_{c}<|k_{z}|<\pi$, obtained
by tight-binding Hamiltonian using only $p_{z}$ orbitals. In real materials
including all core electronic levels, the hinge modes should appear in the
region $0<|k_{z}|<k_{c}$ where $w_{2}=1$.


# Roughness-induced friction in liquid foams

Manon Marchand  Université Paris-Saclay, CNRS, Laboratoire de Physique des
Solides, 91405, Orsay, France Frédéric Restagno  Université Paris-Saclay,
CNRS, Laboratoire de Physique des Solides, 91405, Orsay, France Emmanuelle Rio
Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405,
Orsay, France François Boulogne  Université Paris-Saclay, CNRS, Laboratoire de
Physique des Solides, 91405, Orsay, France

###### Abstract

Complex liquids flow is known to be drastically affected by the roughness
condition at the interfaces. We combined stresses measurements and
observations of the flow during the motion of different rough surfaces in dry
liquid foams. We visually show that three distinct friction regimes exists:
slippage, stick-slip motion, and anchored soap films. Our stress measurements
are validated for slippage and anchored regimes based on existing models, and
we propose a leverage rule to describe the stresses during the stick-slip
regime. We find that the occurrence of the stick-slip or anchored regimes is
controlled by the roughness factor, defined as the ratio between the size of
the surface asperities and the radius of curvature of the Plateau borders.

The flow of complex fluids – polymers [1, 2], emulsions [3], granular
materials [4, 5], bubbly liquids [6], foams [7, 8] – near a solid surface
displays different behaviors depending on the surface properties. Often, these
systems exhibit a slippage on smooth surfaces but deform and flow on rough
surfaces. Therefore, asperities are commonly added on the measurement
apparatus to suppress the interfacial dissipations [9, 10, 11]. Here, we focus
our study on the flow of foams for which the elementary elements, the bubbles,
can be observed directly.

The friction of an elongated bubble in a smooth round capillary has been
studied in the pioneering work of Bretherton at low capillary numbers
$\mathrm{Ca}=\mu V/\gamma$, where $\mu$ and $\gamma$ are the dynamic viscosity
and air-liquid surface tension of the solution, and $V$ is the meniscus
velocity in the capillary [12]. It is shown experimentally and theoretically
that the stress follows a $\mathrm{Ca}^{2/3}$ power law. This prediction is
also valid for a single bubble in a flat cell [13], and can be extended for
foams [14]. The surface rheology of the foaming solution influences the power
law, but in this study we work with a solution for which this issue is
negligible [14, 15]. Also, experimental measurements provide a correction
accounting for the liquid fraction of the foam $\varphi_{\ell}$ [16].

In rheology, the wall slippage of the probed material must be avoided to
ensure that the dissipation is localized in the bulk of the sample. In the
case of foams, a commonly adopted rule-of-thumb to prevent slippage consists
in adding asperities larger than the bubble size on the walls [17], and to
assert either visually [18] or by showing that the rheological measurements do
not depend on the confinement [15, 10], that there is no remaining slippage.
This trial and error process is not always possible and the characterization
of the effect of different roughness sizes has been studied for single wet
bubbles on walls [6], but not for dry liquid foams.

In this Letter, we explore systematically the effect of different sizes of
wall asperities $a$ on the flow of dry liquid foams and we identify three
friction regimes at zero, high or intermediate roughness. The contribution of
the roughness size $a$, of the liquid fraction $\varphi_{\ell}$, and of the
inserting velocity in the foam $V$ to the stresses are rationalized in each
regime. The roughness size is shown to change the stresses over one order of
magnitude. We show that the roughness factor $a/r_{\mathrm{pb}}$ is a single
criterion for the crossover between wall and bulk dissipation regimes.
Interestingly, this parameter is independent of the imposed velocity.

Our experiment consists in inserting horizontally a controlled rough surface
in a dry monodisperse foam generated in a transparent container (Fig. 1a and
b). We produce the foam by blowing air through needles in a foaming solution.
The obtained bubble radius is $R=0.70\pm 0.05$ mm, and the liquid fraction
$\varphi_{\ell}$ of the foam can be varied between $0.03$ and $8$ % by
changing the working height above the liquid-foam interface [19] (Fig. S2).
For such dry foams, the radius of curvature of the Plateau borders is given by
$r_{\mathrm{pb}}=R\sqrt{\varphi_{\ell}/0.33}$ [20]. It varies in the range
$0.02-0.4$ mm in our experiment. Model rough surfaces are obtained by gluing
glass beads of mean radii $a$ on microscope glass slides (Fig. S1), Fig. 1a).
We insert the surfaces in foams at constant velocities $V$ from $0.5$ to $20$
mm/s, corresponding to capillary numbers $\mathrm{Ca}$ in the range $1.7\times
10^{-5}$ to $6.7\times 10^{-4}$. Simultaneously, we measure the tangential
force $F$ exerted by the foam on the surface (Fig. S3).

Figure 1: Foam flow dynamics near a solid surface decorated with glass beads.
(a) Side view of the surface of roughness size $a$ moving at constant speed
$V$ in a foam of bubble radius $R$ and radius of curvature of Plateau borders
$r_{\mathrm{pb}}$. The force $F$ exerted by the foam on the surface is
recorded. (b) Photograph corresponding to the schematic (a). The scale bar is
$5$ mm. (c) Raw force measurements as a function of time $t$ and position of
the surface $p$ for one smooth and two patterned surfaces translating at $V=1$
mm/s from left to right. (d) Visualization of the Plateau borders in contact
with the patterned surfaces for three different $a/r_{\mathrm{pb}}$. The
bottom of each picture is the surface moving at $1$ mm/s. We mark the position
$p_{\mathrm{pb}}$ of a Plateau border with a circle. A plastic event is
circled in white. The scale bar represents $2$ mm. (e) Position of the Plateau
borders $p_{\mathrm{pb}}$ for the three friction regimes, reported from (d).

By tracking the position of the Plateau borders $p_{\mathrm{pb}}$ in contact
with the patterned surface during the motion, we identify three distinct
friction regimes (Fig. 1d), which are discriminated by the dimensionless
roughness factor $a/r_{\mathrm{pb}}$. The emergence of this parameter can be
intuited. When a bead is smaller than the size of the Plateau border
$r_{\mathrm{pb}}$, it enters the Plateau border and slides inside it whereas
bigger beads experience a pinning force at the liquid air interface (Fig.2e).
For negligible roughness sizes compared to the characteristic size of a
Plateau border, $a\ll r_{\mathrm{pb}}$, the Plateau borders slide on the solid
surface (Supplementary Movie S1). For $a\sim r_{\mathrm{pb}}$, Plateau borders
are anchored on the surfaces but occasionally jump back (Supplementary Movie
S2), while for $a\gtrsim r_{\mathrm{pb}}$ Plateau borders remain anchored and
the stress is released by plastic events between bubbles in the bulk of the
foam (Fig. 1d, e, Supplementary Movie S3). We will respectively refer to these
regimes as sliding, stick-slip, and anchored regimes in the following
discussion. For $a/r_{\mathrm{pb}}$ between $1$ and $3$, the visualizations
reveal that the stress is released by mixed slip or plastic events. This
provides a first estimate of the transition between the stick-slip regime and
the anchored one.

Typical force measurements as a function of time $t$, or equivalently of the
position of the surface $p$, are presented for each regime in Fig. 1c. The
linear evolution of the force indicates that the probed phenomena are
independent of the penetration depth of the surface, and that a constant
stress value $\tau_{\mathrm{p}}$ can be extracted from $F(t)$ and $p(t)$ for
each experiment (Fig. S5).

Figure 2: Stress measurements and identification of the three friction
regimes. (a) Stress $\tau_{\mathrm{p}}$ obtained from the force measurements
as a function of the roughness factor $a/r_{\mathrm{pb}}$. The solid black
line represents the equation (3) for $V\rightarrow 0$, and the dashed black
line is $\tau_{\mathrm{ys}}$ with $K=2.9$. The gray area shows the transition
between the stick-slip and the anchored regime as deduced from observation of
the videos, and the orange area is the transition obtained by optimizing $k$
in Eq. (3). (b) For $a/r_{\mathrm{pb}}=0$, stress measurements as a function
of the equation (1). (c) Stress values for the sliding (green circles) and the
stick-slip (blue triangles) conditions, as a function of equation (3). (d) In
the anchored regime, measured stresses as a function of the imposed velocity
$V$. The solid line is the mean of all stresses $\tau_{\mathrm{ys}}$ and the
gray area represents the standard deviation. The deviation for the higher
velocities indicates a more significant contribution of the viscous term in
the Hershell-Bulkley model. (e) Schematics of a Plateau border moving on a
rough surface. The beads spend a time proportional to $a/V$ going through the
liquid air interface, and a time proportional to $r_{\mathrm{pb}}/V$ in the
bulk of the Plateau border.

We reproduce the experiment for various roughness sizes $a$, radius of
curvature of Plateau borders $r_{\mathrm{pb}}$, and insertion velocities $V$.
In Fig. 2a, we plot the stress $\tau_{\mathrm{p}}$ as a function of the
roughness factor $a/r_{\mathrm{pb}}$. The different point colors represent the
three regimes identified from the images. In this representation, the data for
a single velocity collapses on a trend curve (Fig. S4). Complementary
representations of the stick-slip and the anchored regimes presenting the
influence of the asperity size $a$ and the velocity $V$ are provided in Fig.
S7. Interestingly, tuning $a/r_{\mathrm{pb}}$ leads to a stress variation of
about an order of magnitude, which highlights the crucial effect of the
boundary conditions on stresses. In the following, we examine each regime to
rationalize the dynamics.

The sliding regime, for which $a\ll r_{\mathrm{pb}}$, has already been
investigated [12, 16, 14], and these studies validate experimentally that the
stress at the solid wall can be written

$\tau_{\mathrm{slip}}=3.8\frac{\gamma}{R}\varphi_{\ell}^{-0.25}\mathrm{Ca}^{2/3}.$
(1)

In Fig. 2b, the stress measurements are plotted as a function of the
prediction given by equation (1), which shows a good agreement and validates
this experimental approach.

Now, we consider the anchored regime for the large roughness factors. The
common description of the foam rheology suggests that for shear stresses lower
than a yield stress $\tau_{\mathrm{ys}}$, the foam behaves as an elastic
solid, whereas for stresses exceeding the yield stress, the foam flows.
Phenomenologically, this behavior is described by a Hershell-Bulkley law
$\tau_{\mathrm{HB}}=\tau_{\mathrm{ys}}+\eta(\dot{\gamma})\dot{\gamma}$, where
$\dot{\gamma}\sim V/\delta$ is the shear rate that scales as the velocity over
a characteristic shear length $\delta$ and $\eta(\dot{\gamma})$ is the
effective foam viscosity that depends on $\dot{\gamma}$ [21, 11]. The ratio of
the yield stress and the viscous term defines the Bingham number
$\textrm{Bi}=\frac{\tau_{\mathrm{ys}}}{\eta(\dot{\gamma})\dot{\gamma}}$ [22].
The yield stress value is described by a semi-empirical law
$\tau_{\mathrm{ys}}=K\frac{\gamma}{R}(\varphi_{\mathrm{c}}-\varphi_{\ell})^{2}$,
where $\varphi_{\mathrm{c}}=0.26$ is the fraction of gaps remaining in a
close-packing of hard spheres [19], and $K$ is a proportionality factor which
reported values vary between $0.5$ and $6$ [23, 24]. Since we are working with
dry foams, we have $\varphi_{\ell}\ll\varphi_{\mathrm{c}}$, which leads to

$\tau_{\mathrm{ys}}=K\frac{\gamma}{R}\varphi_{\mathrm{c}}^{2}.$ (2)

In the limit of low velocities, corresponding to $\mathrm{Bi}>1$, the
contribution of the yield stress is dominant (see SI). The Fig. 2d shows that
the stress applied by the foam on the rough surfaces is nearly independent on
the velocity and its value is predicted by the equation (2) with a prefactor
$K=2.9\pm 0.3$, which is in agreement with previously recorded values [23,
24]. We notice that the measured stress increases slightly for the larger
velocities, which is reminiscent of the viscous term in the Hershell-Bulkley
model.

Figure 3: Edible foams and texture sensation. (a) Albumen foam in a spoon. The
scale bar is $1$ cm. A drop of food dye is added for visualization. (b) Rough
surface with beads of mean radius $a=225$ $\mu$m. The scale bar is $5$ mm. (c)
Close up picture that shows the papillae on a human tongue. The scale bar is
$5$ mm. (d) Albumen foam sheared between a table spoon and a glass surface
covered with beads of mean radius $a=225$ $\mu$m. The left side is a wet foam
with $\varphi_{\ell}\simeq 0.7$ %, corresponding to $\tilde{a}\simeq 0.5$. The
right side is a dryer foam with a liquid fraction $\varphi_{\ell}\simeq 0.02$
%, meaning $\tilde{a}\simeq 3$. In those two cases the last bubble on the
surface is circled to enlighten the motion of foam on the model tongue. There
is a $1$ s interval between two consecutive images. The scale bar is $1$ cm.

Up to this point, we analyzed the two extreme limits of negligible and large
dimensionless roughness factor $a/r_{\mathrm{pb}}$. For roughness sizes
comparable to the radius of the Plateau borders, we observe a stick-slip
regime that has been observed but not characterized [25], to the best of our
knowledge. Our experiments indicate that the transition between the stick-slip
and the anchored regimes is continuous (Fig. 2a). This suggests that the
stick-slip regime is a combination of sliding, where the energy is released by
viscous dissipation in the Plateau borders, and anchoring, where the energy is
stored by deforming elastically the foam.

In the stick-slip regime, (Fig. 2e) an asperity moving in a foam goes through
the liquid air interface of a Plateau border over a duration
$t_{\mathrm{anchored}}\propto a/V$, experiencing a pinning force. Then, the
asperity moves inside the Plateau border over a typical duration
$t_{\mathrm{slip}}\propto r_{\mathrm{pb}}/V$. Therefore, we introduce the
ratio of these durations, $\tilde{a}=a/kr_{\mathrm{pb}}$, where the prefactor
$k$ represents the value of the transition. This two-step description happens
for each asperity on the whole probing surface. The stick-slip events are not
occurring simultaneously (Supplementary Movie S2), thus we average the
individual behaviors of the Plateau borders. Hence, we construct a leverage
rule stating that the total stress is the sum of the contributions of a
pinning stress for a relative duration $\tilde{a}$ and of a viscous stress for
the complementary relative duration $(1-\tilde{a})$. This writes

$\tau_{\mathrm{stick-
slip}}=\tau_{\mathrm{ys}}\tilde{a}+\tau_{\mathrm{slip}}(1-\tilde{a}).$ (3)

To determine more precisely the value of the transition $k$, we optimized
equation (3) against this parameter and we obtained a transition at $k=1.9\pm
0.2$ (Fig. S6a). This value of order unity is in agreement with our estimation
of the transition from direct visualizations of the Plateau border dynamics
(see the orange and gray areas in Fig. 2a). As shown in Fig. 2c, the proposed
equation (3) is in excellent agreement with the stress measurements for the
different liquid fractions and surface velocities explored in this study. In
the vicinity of small velocities, we expect $\tau_{\mathrm{stick-
slip}}\underset{V\to 0}{\rightarrow}\tau_{\mathrm{ys}}\tilde{a}$, and this is
observed in Fig. 2a and Fig. S6b.

Thus, we identified a practical criterion $\tilde{a}$ to discriminate between
surface and bulk dissipation for the friction of a foam on a solid surface.
The anchored regime is of particular interest for rheological studies focusing
on the flow of foams where a no slip condition is necessary to transfer shear
to the bulk of the material. It was already known that adding some asperities
by gluing sand paper or engraving grooves on a measuring apparatus reduces the
wall slippage [17], but the size of these asperities had not been tuned
systematically. Our study validates this approach, since the stress in the
anchored regime are independent of the size of the asperities beyond a
critical roughness factor (Fig. 2a). When the grit size of the sand paper is
insufficient to totally eliminate the wall slippage, the mixed surface and
bulk dissipation leads to a difficult interpretation of the data [8].
Differences in the measured stresses are also reported when the grit size is
increased [26]. This behavior is likely the stick-slip regime we report in
Fig. 2c. Our contribution allows one to choose a suitable roughness to be
added on the measuring systems regarding the foam properties. The results then
assess that the obtained measurements in the anchored regime will be
uncorrelated with the size of the roughness.

The stick-slip regime allows, for instance, to understand the mastication of
an aerated food (Fig. 3a). Adding air bubbles in a food product affects
mouthfeel perception, improves digestibility, and aids mastication [27]. In
the case of eating disorders, such as dysphagia, more viscous foods are proven
to be safer [28]. With viscosities typically $10^{3}$ times bigger than the
viscosity of water [29], foams could help patients. We make a visualization
experiment to illustrate the tunability of perception during the mastication
of a foamed food. Human tongue asperities have a fixed size, therefore the
remaining adjustable parameters in this case are the bubble size and the
liquid fraction of the foam. A glass slide covered with asperities of mean
radius $a=225$ $\mu$m, a size comparable to human papillae [30], can be
considered as a model human tongue (Fig. 3a, b and c). We produce two edible
albumen foams with a bubble size $R=0.7$ mm and two different liquid
fractions, using the same technique than soap foams, but with a solution
containing $1$ g of egg white powder (Louis François) and $200$ mL of
deionized water. The two foams are sheared between a table spoon moving at
$10$ mm/s and the static model tongue. In the two cases, the last bubble is
circled to illustrate the slippage of the foam (Fig. 3d). For the wetter foam
for which $\tilde{a}\simeq 0.5<1$, the conditions predict a stick-slip regime
and the foam moves on the surface (Supplementary Movie S4), whereas for the
dryer foam, for which $\tilde{a}\simeq 3>1$, the foam sticks to the surface
(Supplementary Movie S5). This illustrates that the results of our controlled
experiment are valid in a more general context, with different shearing
geometries and surfactants. When the normalized roughness $\tilde{a}$ varies
from $0$ to $1$ in the stick-slip regime, we see in Fig. 2a that the stress
varies between $0.8$ and $9.8$ Pa, a range of stresses detectable by a human
tongue [31, 32]. Thus, tuning the macroscopic properties of a product, namely
the bubble size and the liquid-gas ratio, must lead to different sensations in
mouth.

In conclusion, we reveal that the foam flow near surfaces exhibit different
behaviors: sliding, anchored, and stick-slip. By recording the stresses
exerted during these regimes, we show that there exists a roughness,
normalized by the curvature radius of the Plateau borders of the foam,
$\tilde{a}=a/kr_{\mathrm{pb}}$ beyond which the dissipation is transferred
from the surface to the bulk. This transition does not depend on the imposed
velocity in the explored range. In the stick slip regime, we propose a
leverage rule describing the influence of speed, radius of curvature of
Plateau borders, and liquid fraction on the stresses exerted on the walls. We
measure that these stresses vary on one order of magnitude when the normalized
roughness is increased. Therefore, we contribute to an improved description of
wall roughness effect on flowing foams, which is of particular interest to
design surfaces given the foam properties in industrial applications. In
contrast, when surface asperities are unchangeable, as in the case of tongue
papillae, our description of friction forces allows to finely tune the
properties of edible foams to get creamy or gooey mouth feelings. This is
useful for inventing new foamed food products, either for the ease of feeding
for dysphagic patients, or for the gourmets’ pleasure.

## Acknowledgements

M.M. acknowledges EDPIF for funding support. We thank J. Sanchez for
programming the computer interface of the experiment, and G. Guillier for
support manufacturing the force sensor. We thank B. Dollet and C. Raufaste for
discussions.

## References

  * [1] R. D. Priestley, C. J. Ellison, L. J. Broadbelt, and J. M. Torkelson.  Structural relaxation of polymer glasses at surfaces, interfaces, and in between.  Science, 309(5733):456–459, 2005.
  * [2] P. E. Boukany, O. Hemminger, S.-Q. Wang, and L. J. Lee.  Molecular imaging of slip in entangled dna solution.  Phys. Rev. Lett., 105:027802, 7 2010.
  * [3] J. Goyon, A. Colin, G. Ovarlez, A. Ajdari, and L. Bocquet.  Spatial cooperativity in soft glassy flows.  Nature, 454(7200):84—87, 7 2008.
  * [4] S. Siavoshi, A. V. Orpe, and A. Kudrolli.  Friction of a slider on a granular layer: Nonmonotonic thickness dependence and effect of boundary conditions.  Phys. Rev. E, 73:010301, 1 2006.
  * [5] O. Pouliquen.  On the shape of granular fronts down rough inclined planes.  Physics of Fluids, 11(7):1956–1958, 1999.
  * [6] D. Germain and M. Le Merrer.  Bubbles slipping along a crenelated wall.  EPL, 115(6):64005, 2016.
  * [7] A. Kabla and G. Debrégeas.  Local stress relaxation and shear banding in a dry foam under shear.  Phys. Rev. Lett., 90:258303, 6 2003.
  * [8] K. Golemanov, S. Tcholakova, N. D. Denkov, K. P. Ananthapadmanabhan, and A. Lips.  Breakup of bubbles and drops in steadily sheared foams and concentrated emulsions.  Phys. Rev. E, 78:051405, 11 2008.
  * [9] H. A. Barnes.  A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure.  Journal of Non-Newtonian Fluid Mechanics, 56(3):221 – 251, 1995\.
  * [10] M. Cloitre and R. T. Bonnecaze.  A review on wall slip in high solid dispersions.  Rheologica Acta, 56(3):283–305, 3 2017.
  * [11] S. Cohen-Addad, R. Höhler, and O. Pitois.  Flow in foams and flowing foams.  Annual Review of Fluid Mechanics, 45(1):241–267, 2013.
  * [12] F. P. Bretherton.  The motion of long bubbles in tubes.  Journal of Fluid Mechanics, 10:166–188, 3 1961.
  * [13] C.-W. Park and G. M. Homsy.  Two-phase displacement in hele shaw cells: theory.  Journal of Fluid Mechanics, 139:291–308, 1984.
  * [14] I. Cantat.  Liquid meniscus friction on a wet plate: Bubbles, lamellae, and foams.  Physics of Fluids, 25(3):031303, 2013.
  * [15] N. D. Denkov, V. Subramanian, D. Gurovich, and A. Lips.  Wall slip and viscous dissipation in sheared foams: Effect of surface mobility.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 263(1):129 – 145, 2005.
  * [16] C. Raufaste, A. Foulon, and B. Dollet.  Dissipation in quasi-two-dimensional flowing foams.  Physics of Fluids, 21(5):053102, 2009.
  * [17] S. A. Khan, C. A. Schnepper, and R. C. Armstrong.  Foam rheology: III. measurement of shear flow properties.  Journal of Rheology, 32(1):69–92, 1988.
  * [18] G. Katgert, B. P. Tighe, M. E. Möbius, and M. van Hecke.  Couette flow of two-dimensional foams.  EPL, 90(5):54002, 6 2010.
  * [19] A. Maestro, W. Drenckhan, E. Rio, and R. Hohler.  Liquid dispersions under gravity: volume fraction profile and osmotic pressure.  Soft Matter, 9:2531–2540, 2013.
  * [20] S. A. Koehler, S. Hilgenfeldt, and H. A. Stone.  A generalized view of foam drainage: Experiment and theory.  Langmuir, 16(15):6327–6341, 2000.
  * [21] N. D. Denkov, S. Tcholakova, K. Golemanov, K. P. Ananthpadmanabhan, and A. Lips.  The role of surfactant type and bubble surface mobility in foam rheology.  Soft Matter, 5:3389–3408, 2009.
  * [22] E. C. Bingham.  An investigation of the law of plastic flow.  Bulletin of the Bureau of Standards, 13:309–353, 6 1916.
  * [23] M. Lexis and N. Willenbacher.  Yield stress and elasticity of aqueous foams from protein and surfactant solutions – the role of continuous phase viscosity and interfacial properties.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 459:177 – 185, 2014.
  * [24] F. Rouyer, S. Cohen-Addad, and R. Höhler.  Is the yield stress of aqueous foam a well-defined quantity?  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 263(1):111 – 116, 2005.
  * [25] M. Buchgraber, L. M. Castanier, and A. R. Kovscek.  Microvisual investigation of foam flow in ideal fractures: Role of fracture aperture and surface roughness.  SPE Annual Technical Conference and Exhibition, 1 2012.
  * [26] C. A. Jimenez-Junca, J. C. Gumy, A. Sher, and K. Niranjan.  Rheology of milk foams produced by steam injection.  Journal of Food Science, 76(9):E569–E575, 2011.
  * [27] G. M. Campbell and E. Mougeot.  Creation and characterisation of aerated food products.  Trends in Food Science & Technology, 10(9):283 – 296, 1999.
  * [28] K. Nishinari, M. Turcanu, M. Nakauma, and Y. Fang.  Role of fluid cohesiveness in safe swallowing.  NPJ Science of Food, 3(5), 2019.
  * [29] S. P. L. Marze, A. Saint-Jalmes, and D. Langevin.  Protein and surfactant foams: linear rheology and dilatancy effect.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 263(1):121 – 128, 2005.
  * [30] G. K. Essick, A. Chopra, S. Guest, and F. McGlone.  Lingual tactile acuity, taste perception, and the density and diameter of fungiform papillae in female subjects.  Physiology & Behavior, 80(2):289 – 302, 2003.
  * [31] E. A. Foegeding, C. J. Vinyard, G. Essick, S. Guest, and C. Campbell.  Transforming structural breakdown into sensory perception of texture.  Journal of Texture Studies, 46(3):152–170, 2015.
  * [32] B. Linne and C. T. Simons.  Quantification of oral roughness perception and comparison with mechanism of astringency perception.  Chemical Senses, 42(7):525–535, 06 2017.