0.8pt

# Elastic ribbons in bubble columns:
when elasticity, capillarity and gravity govern equilibrium configurations

Jean Farago1 Manon Jouanlanne1 Antoine Egelé1 Aurélie Hourlier-Fargette1 1
Université de Strasbourg, CNRS, Institut Charles Sadron UPR22, F-67000
Strasbourg, France

###### Abstract

Taking advantage of the competition between elasticity and capillarity has
proven to be an efficient way to design structures by folding, bending, or
assembling elastic objects in contact with liquid interfaces. Elastocapillary
effects often occur at scales where gravity does not play an important role,
such as in microfabrication processes. However, the influence of gravity can
become significant at the desktop scale, which is relevant for numerous
situations including model experiments used to provide a fundamental physics
understanding, working at easily accessible scales. We focus here on the case
of elastic ribbons placed in two-dimensional bubble columns: by introducing an
elastic ribbon inside the central soap films of a staircase bubble structure
in a square cross-section column, the deviation from Plateau’s laws
(capillarity-dominated case dictating the shape of usual foams) can be
quantified as a function of the rigidity of the ribbon. For long ribbons,
gravity cannot be neglected. We provide a detailed theoretical analysis of the
ribbon profile, taking into account capillarity, elasticity and gravity. We
compute the total energy of the system and perform energy minimization under
constraints, using Lagrangian mechanics. The model is then validated via a
comparison with experiments with three different ribbon thicknesses.

††preprint: APS/123-QED

##  I Introduction

Systems involving both elastic and capillary forces have been attracting a
growing interest to tackle practical challenges and provide novel engineering
and design tools, compelling for the microfabrication of 2D surfaces and 3D
structures [1, 2], for the design of actuators in the emerging soft robotics
field [3], but also in the context of biomechanics or biomimetic approaches
[4, 5, 6], to cite only a few cases. The field of elastocapillarity includes
not only the deformation of soft solids by surface tension [7] but also
systems where the geometry of the elastic objects is of particular importance,
especially in the context of slender structures for which bending dominates
over stretching [8]. In appropriate conditions, capillary forces are
sufficient to force such slender structures to undergo large deformations, in
systems at multiple length scales, provided that capillary and elastic
energies are of the same order of magnitude [9, 10]. A fundamental
understanding of elastocapillary phenomena can be acquired with experiments at
the desktop scale [11, 12], for systems conceptually very close to those
occurring in microfabrication processes [13, 14]. However, when reaching
centrimetric scales, the effects of gravity on slender elastic structures can
be significant. Taking capillarity aside, a typical daily life experience of
slender structures for which gravity matters is the shape of curly hair [15].
A few examples involving elasticity, capillarity and gravity have been studied
in the literature, for instance in the context of floating thin elastic films
[16], capillary in-drop spooling [17], elastocapillary-driven snap-through
[18], or of capillary rise between elastic sheets, revisiting Jurin’s law in
the context of elastocapillarity [19].

Here we focus on the question of the influence of gravity on the equilibrium
shape of centimetric ribbons inserted inside two dimensional bubble columns.
When confined inside square cross-section tubes, bubbles rearrange themselves
into well ordered structures that depend on the confinement ratio (bubble size
divided by the width of the tube) [20, 21]. By choosing an appropriate ratio,
we can obtain the so-called staircase structure which is quasi 2D (Fig.1a). A
ribbon can then be inserted inside the central soap films, and its shape
depends on the minimisation of the total energy of the system. In our previous
work [22, 23], we placed ourselves in conditions where gravity is negligible,
by considering only the very bottom of the elastic ribbons. The total energy
of the system was thus composed of the bending energy of the ribbon and of the
interfacial energies, and we studied how the presence of the ribbon modified
the classical Plateau’s laws [24] for foams. Now, we consider long ribbons for
which a flattening due to gravity is observed, especially in the upper parts,
due to the weight of the lower parts.

##  II  Model description and notations

Figure 1: Illustration of the system studied: (a) Schematic representation of
a bubble column in a tube of square cross-section. The ratio of the bubble
size to the container dimension is chosen to obtain the so-called staircase
structure [20, 21], invariant by translation along the axis perpendicular to
the tube. (b) Experimental visualisation of the insertion of an elastic ribbon
(polydimethylsiloxane (PDMS)) in a staircase bubble structure. (c) X-ray
microtomography slice of the bubble column containing a thin PDMS ribbon
(thickness $t=$45 $\unit{\micro}$), and zoom on a half-cell presenting the
definitions of the half bubble size $\ell$ and of the amplitude $\Delta f$.
Note that the half-cell numbers are defined by starting from the bottom of the
ribbon.

As shown in Figure 1, we describe a system similar to that studied in [22]. An
elastic ribbon of length $L$ is inserted in the middle of a 2D staircase
arrangement of bubbles of longitudinal size $2\ell$. We divide the ribbon in
half-cells (HC), with numbers defined by starting from the bottom of the
ribbon (HC n°1 corresponds to the first half period containing a ribbon in the
whole central Plateau border). The ribbon is characterized by a bending
modulus per transverse length $\alpha=Et^{3}/(12(1-\nu^{2}))$ where $t$ is the
ribbon thickness, and ($E,\nu$) respectively the Young modulus and the Poisson
ratio of the ribbon, and a lineic mass per transverse length $\lambda$
($\lambda=\rho t$ with $\rho$ the ribbon mass density). Contrary to the case
studied in [22], the gravity $g$ operating in the longitudinal $Oz$ direction
is not supposed to be negligible. Its main effect is to flatten the upper
regions of the ribbon due to the weight of the lower regions pulling
downwards. The purpose of the present theory is to compute the equilibrium
shape of a long ribbon, tied at its upper end at $z=0$, in the limit of a long
ribbon $L/\ell\gg 1$ (to avoid finite size effects difficult to handle, which
moreover would obscure the understanding of the theory).

The elastocapillary energy excess $\Delta E$ per unit (transverse) length of
the system, with respect to the situation where only the bubbles are present
is given by

$\Delta E=\frac{\alpha}{2}\int_{0}^{L}{\rm d}s\ C^{2}-\lambda
g\int_{0}^{L}{\rm d}s\ z(s)\\\ +2\gamma L-\gamma\sum_{\rm HC}[\Delta
f_{i}+\sqrt{3}\ell]$ (1)

where $s$ is the curvilinear coordinate along the ribbon, starting from $s=0$
at $z=0$, $C=f^{\prime\prime}(z)/\sqrt{1+[f^{\prime}(z)]^{3/2}}$ is the local
ribbon curvature (we denote $y=f(z)$ the profile of the ribbon in the lateral
direction), and $\Delta f_{i}$ is the lateral amplitude of the ribbon
deformation in the $i$-th HC (we choose to number the HC from the bottom
($n=1$) upward, until $n=N$). The first three terms of the right hand side
(rhs) relate to the ribbon, as respectively its bending, gravitational and
interfacial energies, the latter constant term assuming a perfect wetting of
the ribbon. The second line of (1) is the excess interfacial energy the system
has to pay to substitute in each HC the original liquid staircase pattern by a
wet ribbon : if $h$ is the width of the lateral dimension of the tube (i.e.
visible in fig. 1 c), this energy difference in the $i$-th HC is indeed
$\gamma[2L+h-\Delta
f_{i}-(\frac{4\ell}{\sqrt{3}}+h-\frac{\ell}{\sqrt{3}})]=2\gamma
L-\gamma[\Delta f_{i}+\sqrt{3}\ell]$. The $\sqrt{3}\ell$ term must not be
discarded as a constant, due to the fact that the number of HC is variable in
the optimization process.

In writing (1), one assumed that the HC size $\ell$ is a constant along the
tube, moreover unaffected by the presence of the ribbon. This approximation is
justified in the Appendix at the end of the article. The direct minimization
of (1) is complicated, because the gravitational term breaks the vertical
translational invariance. However the weight of the ribbon contained in just
one HC is very modest, and it is relevant to assume that the gravitational
energy may be considered as constant at the scale of the HC. This assumption
assumes that the weight (per unit transverse length) of just one HC,
$\sim\lambda g\ell$ is relatively small with respect to the typical forces
(per unit transverse length), which can be either $\gamma$ or
$\alpha/\ell^{2}$ according to the regime considered (see later). This
assumption amounts to replace $-\lambda g\int_{{\rm HC}n}{\rm d}s\ z(s)$ by
$-\lambda g(N-n+\frac{1}{2})\ell\int_{{\rm HC}n}{\rm d}s$, where the integral
encompasses the part of the ribbon present in the $n$-th HC. In doing so, one
assumes that the center of mass of the ribbon stays always in the middle of
the HC, which is in particular provided by configurations where the ribbon
shape is centrosymmetric within each HC. This would be exactly the case for
the optimal ribbon shape in the $n$-th HC if its weight was negligible,
because in this case, the gravitational pulling forces exerted by the
neighbouring HCs on each end of the ribbon part contained in the HC would be
strictly identical.

For a long ribbon, a large upper amount is totally flattened by the weight of
its lower part, so it is relevant to remove from $\Delta E$ the constant
energy of a totally flattened ribbon, $(\Delta E)_{\rm flat}=2\gamma
L+\sum_{n=1}^{L/\ell}[-\lambda
g(\frac{L}{\ell}-n+\frac{1}{2})\ell^{2}-\gamma\sqrt{3}\ell]$, and define

$\displaystyle\Delta\mathcal{E}$ $\displaystyle=\Delta E-(\Delta E)_{\rm
flat}+\gamma\Lambda\left[\int{\rm d}s-L\right]$ (2)

where $\Lambda$ is a dimensionless Lagrangian multiplier added to encode the
constraint that the total length of the ribbon is conserved. If $N$ denotes
the total number of HC the ribbons spans, the preceding expression can be
recast into :

$\displaystyle\Delta\mathcal{E}$
$\displaystyle=\sum_{n=1}^{N}\left[\frac{\alpha}{2}\int_{{\rm HC}n}{\rm d}s\
C^{2}-\gamma\Delta f_{n}\right.$ $\displaystyle\phantom{a}\left.+[\lambda
g\ell(n-\tfrac{1}{2})+\gamma\Lambda]\left(\int_{{\rm
HC}n}\\!\\!\\!\\!\\!\\!\\!\\!{\rm
d}s-\ell\right)\right]+\gamma\ell\mathcal{R}_{N}$ (3)
$\displaystyle\mathcal{R}_{N}$ $\displaystyle=\frac{\lambda
g\ell}{2\gamma}(L/\ell-N)^{2}+(\sqrt{3}-\Lambda)(L/\ell-N)$ (4)

With this writing, one can expect that the bracketed series be convergent if
$N\rightarrow\infty$, since the large values of ${n}$ will be cancelled by the
terms $\propto(\int{\rm d}s-\ell)$, very small at the top of the ribbon (we
will see that a remaining logarithmic divergence remains, however negligible
in the limit $N\gg 1$). For sake of simplicity, the $L/\ell\rightarrow\infty$
(which entails $N\rightarrow\infty$) limit will be systematically considered
in the following. It corresponds to taking the asymptotic sum of the series in
(3), but keeping in mind that $L/\ell-N$ remains finite. If finite $N$
corrections would have to be considered, then the boundary effects associated
to the extremities of the ribbon ending in general in the middle of a HC
should be carefully estimated, a more detailed analysis which is beyond the
scope of this work.

We make the problem dimensionless by defining

$\displaystyle\widehat{\lambda}$ $\displaystyle=\frac{\lambda g\ell}{\gamma}$
(5) $\displaystyle\eta$ $\displaystyle=\frac{\alpha}{\gamma\ell^{2}}$ (6)

and writing (note that $\Delta\mathcal{E}$ is an energy per transverse unit
length) :

$\displaystyle\Delta\mathcal{E}/[\gamma\ell]$
$\displaystyle=\sum_{n=1}^{N}e_{n}+\mathcal{R}_{N}$ (7) $\displaystyle e_{n}$
$\displaystyle=\frac{\eta}{2}\int_{{\rm HC}n}\frac{{\rm
d}s}{\ell}(C\ell)^{2}+\eta\kappa_{n}^{2}\left(\int_{{\rm HC}n}\frac{{\rm
d}s}{\ell}-1\right)-\frac{\Delta f_{n}}{\ell}$ (8)
$\displaystyle\kappa_{n}^{2}$
$\displaystyle=\frac{\widehat{\lambda}({n}-\frac{1}{2})+\Lambda}{\eta}$ (9)

In the preceding formula, the curvature $C$ is related to the profile $f(z)$
of the ribbon by the already mentioned classical formula
$C=f^{\prime\prime}(z)/\sqrt{1+[f^{\prime}(z)]^{3/2}}$, and ${\rm
d}s=\sqrt{df^{2}+dz^{2}}=dz\sqrt{1+[f^{\prime}(z)]^{2}}$. The optimization of
(7) consists in (i) writing the Euler-Lagrange equations for $f(z)$ of the
ribbon in each HC separately, which is possible since the optimal solution has
extrema at the boundaries of the HC, (ii) Express the constraint “ribbon
length$=\\!L$” to remove $N$ from (7) in favour of $\Lambda$, and (iii)
Extremalize the resulting expression for $\Delta\mathcal{E}$ with respect to
$\Lambda$, which is equivalent to (but easier than) extremalizing with respect
to $N$. Notice that the whole optimal solution is twice differentiable
everywhere with slight discontinuities of the second derivatives at the
boundaries of HC, which is not forbidden by the physics, as the lateral liquid
interfaces act as localized (“impulsive”) normal forces at these boundaries on
the ribbon.

In the general case, the optimization of (7) is probably intractable, but as
in [22], we can separately solve the case $\eta\gg 1$ (stiff ribbons), where
the ribbon remains everywhere almost flat, and $\eta\ll 1$ (soft ribbons),
where the ribbon perturbs only slightly the Plateau’s laws. In the latter
case, a slight difficulty must be overcome, since the gravity significantly
flattens the bubbles in the upper parts of a long enough ribbon, even if the
ribbon elasticity does not.

###  II.1 Stiff ribbon limit : $\eta\gg 1$

This case is certainly less relevant experimentally, but it is theoretically
instructive, as much simpler to proceed with. Here $f\ll\ell$ everywhere and
it is relevant to Taylor-expand $e_{n}$ up to the second order :

$\displaystyle e_{n}$
$\displaystyle\simeq\frac{\eta\ell}{2}\int_{0}^{\ell}{\rm
d}z[f^{\prime\prime}(z)]^{2}+\frac{\eta\kappa_{n}^{2}\ell^{-1}}{2}\int_{0}^{\ell}{\rm
d}z[f^{\prime}(z)]^{2}-\frac{\Delta f_{n}}{\ell}$ (10)

(notice that a translation in $z$ is made so that $z$ spans the interval
$[0,\ell]$ whatever the HC considered). One writes $f(z)=(\Delta
f_{n})g_{n}(z/\ell)$, since the particular structure of (10) allows a separate
optimization of the shape $g_{n}(z/\ell)$ and the amplitude $\Delta f_{n}$ of
the profile. Solving the 4-th order Euler-Lagrange equation with the four
boundary conditions
$(g_{n}(0)=0,g_{n}(1)=1,g_{n}^{\prime}(0)=g_{n}^{\prime}(1)=0)$ gives for the
optimal shape 111As we consider only half cells (HC), one should impose
$g_{n}(0)=1$ and $g_{n}(1)=0$ every two HC. But this point does not change any
of the conclusions in the adopted assumption of light ribbons, so we impose
$g(0)=0$ and $g(1)=1$ everywhere for sake of simplicity :

$g_{n,\rm opt}(u=\tfrac{z}{\ell})\\\
=\frac{\tanh\left(\frac{\kappa_{n}}{2}\right)[\cosh(\kappa_{n}u)-1]+\kappa_{n}u-\sinh(\kappa_{n}u)}{\kappa_{n}-2\tanh\frac{\kappa_{n}}{2}}$
(11)

To get the optimal value of $\Delta f_{n}$, we express the energy $e_{n}$
optimized with respect to $g_{n}$ as

$\displaystyle e_{n}$ $\displaystyle=\frac{\eta}{2}\left(\frac{\Delta
f_{n}}{\ell}\right)^{2}\int_{0}^{1}{\rm d}ug^{\prime}_{n,{\rm
opt}}[-g^{(3)}_{n,\rm{opt}}+\kappa_{n}^{2}g^{\prime}_{n,{\rm
opt}}]-\frac{\Delta f_{n}}{\ell}$ (12)

By optimizing $e_{n}$ with respect to $\Delta f_{n}$, we obtain

$\displaystyle\frac{\Delta f_{n,\rm opt}}{\ell}$
$\displaystyle=\eta^{-1}\frac{\kappa_{n}-2\tanh\frac{\kappa_{n}}{2}}{\kappa_{n}^{3}}$
(13) $\displaystyle e_{n,\rm opt}$ $\displaystyle=-\frac{1}{2}\frac{\Delta
f_{n,\rm opt}}{\ell}$ (14)

As announced in anticipation, it can be observed that the optimal value of the
series in (7) keeps a slight divergence, since $e_{n,\rm
opt}\sim-\frac{1}{2\widehat{\lambda}{n}}$ for large values of ${n}$. This
divergence gives a term $\propto\ln(N)$, which we make explicit by writing

$\displaystyle\sum_{n=1}^{N}e_{n,\rm
opt}\simeq\sum_{n=1}^{\infty}\left[e_{n,\rm
opt}+\frac{1}{2\widehat{\lambda}n}\right]-\frac{\ln N}{2\widehat{\lambda}}$
(15)

The constraint that the ribbon length in the optimal shape is $L$ reads
$L/\ell-N_{\rm opt}=\ell^{-1}\sum_{\rm HC}[\int{\rm ds}-\ell]$, and up to the
second order in $f_{n}$ :

$\displaystyle L/\ell-N=\frac{1}{2}\sum_{n=1}^{N_{\rm opt}}\left(\frac{\Delta
f_{n}}{\ell}\right)^{2}\int_{0}^{1}{\rm d}ug_{n,\rm
opt}^{\prime}(u)^{2}\sim\mathcal{S}$ (16) $\displaystyle\mathcal{S}$
$\displaystyle=\frac{1}{2\eta^{2}}\sum_{{n}=1}^{\infty}\frac{1}{\kappa_{n}^{5}}\left[-\tfrac{\kappa_{n}}{2}\tanh^{2}(\tfrac{\kappa_{n}}{2})+3\tfrac{\kappa_{n}}{2}-3\tanh(\tfrac{\kappa_{n}}{2})\right]$
(17)

where the summation has been pushed to $+\infty$ since the series is
convergent. Note that $\mathcal{S}=\partial_{\Lambda}\sum(e_{n,{\rm
opt}}+1/(2\widehat{\lambda}n))$, which is understood by considering Eq. (10) :
The derivative of this expression with respect to $\Lambda$ has two parts, one
implicit because $f$ depends on $\Lambda$, which is identically zero because a
first order variation of the optimal $f$ is always zero whatever its value,
and one explicit because $\kappa_{n}^{2}$ depends on $\Lambda$, which gives
exactly the second-order approximation to $L/\ell-N$, i.e. $\mathcal{S}$. Note
also that we discarded the term $\propto\ln(N)$ which induces negligible
corrections for long enough ribbons. Finally, one uses (4), (7) and (16) to
write

$\displaystyle\Delta\mathcal{E}/\gamma\ell=\sum_{n=1}^{\infty}\left[e_{n,\rm
opt}+\frac{1}{2\widehat{\lambda}n}\right]+\frac{\widehat{\lambda}}{2}\mathcal{S}^{2}+(\sqrt{3}-\Lambda)\mathcal{S}$
(18)

This expression has to be optimized with respect to $\Lambda$, a step
equivalent to finding the optimal number $N_{\rm opt}$ of HC spanned by the
ribbon at equilibrium (via Eq. (16)). The fact that
$\partial_{\Lambda}\sum_{n}[e_{n,\rm
opt}+1/(2\widehat{\lambda}n)]=\mathcal{S}$ shows that the equation fulfilled
by the optimal $\Lambda$ is either $\partial_{\Lambda}\mathcal{S}=0$ (which is
impossible because $e_{n,\rm opt}$ is a convex function of $\Lambda$) or

$\displaystyle\frac{\Lambda-\sqrt{3}}{\widehat{\lambda}}=\mathcal{S}$ (19)

To go further, one can evaluate numerically the preceding expression and solve
for $\Lambda$, but one can also make use of the fact that $\eta\gg 1$. For
small values of $\widehat{\lambda}/\eta=\lambda g\ell^{3}/\alpha$, the
parameter $\kappa_{n}$ is slowly varying with $n$, and we can use the Euler-
MacLaurin formula and replace the sum in (17) with an integral (plus the
correction terms “$\frac{f(1)}{2}-\frac{f^{\prime}(1)}{12}$” [26] to be
comprehensive at the order $\eta^{-4}$): After a cumbersome calculation, one
finds

$\displaystyle\Lambda$
$\displaystyle=\sqrt{3}+\frac{1}{24\eta}-\frac{\sqrt{3}}{240\eta^{2}}+\frac{11-\frac{17}{48}\widehat{\lambda}^{2}}{10080\eta^{3}}+o(\eta^{-3})$
(20)

It is remarkable that the first three orders of the expansion do not depend on
$\widehat{\lambda}$. This expression corresponds to the limit of large $\eta$
at constant $\widehat{\lambda}$, and can be used also in the experimentally
relevant limit of large $t$, with $\eta\propto t^{3}$ and
$\widehat{\lambda}\propto t$ : The term
$\propto\widehat{\lambda}^{2}\eta^{-3}$ must in this case be considered as
$O(t^{-7})$ and overcomes the $O(\eta^{-3})=O(t^{-9})$ one.

Another interesting case corresponds to $\widehat{\lambda}=0$ whatever $\eta$.
In this case, the lowest approximation for $\kappa_{n}$ in the case
$\widehat{\lambda}=0$ is $\kappa_{n}\sim 3^{1/4}/\sqrt{\eta}$, a result we
obtained in [22] (formula (11) in [22], $\kappa_{n}\equiv 2\kappa$. The shape
$g_{n}$ and the transverse amplitude $\Delta f_{n}$ are also the same).

Particularly interesting is the evolution of $\Delta f_{n}/\ell$ with
${n}=O(1)$ (lower part of the ribbon). In the approximation $\eta\gg 1$ we are
considering here, $\kappa_{n}$ is for the lower part of the ribbon a small
quantity and one can approximate (13) by

$\displaystyle\frac{\sqrt{3}\Delta f_{n,\rm opt}}{\ell}$
$\displaystyle=\frac{1}{4\sqrt{3}\eta}\left[1-\frac{\sqrt{3}+\widehat{\lambda}({n}-\frac{1}{2})}{10\eta}\right]$
(21)

As a result, a linear reduction of the lateral amplitude can be expected, with
a slope $\widehat{\lambda}/[40\sqrt{3}\eta^{2}]$ quite small with reasonable
values of $\widehat{\gamma}$ and $\eta\gg 1$. This limit, $\eta\gg 1$,
therefore does not appear as a relevant one to capture a signature of the
gravitational effects. However, the theory is particularly clear. In the next
paragraph, the same line of reasoning is followed in the opposite limit, when
$\eta\ll 1$, a more complicated case due essentially to the fact that it is a
expansion around a singular solution, reminiscent for instance of the zero
temperature expansion of fermionic gas.

###  II.2 Floppy ribbon limit : $\eta\ll 1$

####  II.2.1 The $\eta=0$ case

In this opposite limit, the gravitational pull and the elastocapillary
interaction oppose to each other since the capillarity dominates largely the
forces opposing the bending and forces the ribbon to a winding shape. To
tackle this case, it is useful to note that $\widehat{\lambda}\propto t$ ($t$
is the ribbon thickness), whereas $\eta\propto t^{3}$, which suggests that one
considers first the case where formally $\eta=0$ and $\widehat{\lambda}\neq
0$. We define

$\displaystyle\omega_{n}^{2}$
$\displaystyle=\widehat{\lambda}({n}-\frac{1}{2})+{\Lambda}$ (22)

instead of $\kappa_{n}^{2}=\omega_{n}^{2}/\eta$. The energy terms now read

$\displaystyle e_{n}^{\circ}$
$\displaystyle=\omega_{n}^{2}\int_{0}^{\ell}\frac{{\rm
d}z}{\ell}[\sqrt{1+[f^{\prime}(z)]^{2}}-1]-\frac{\Delta f_{n}}{\ell}$ (23)

The optimal shape is readily $g^{\circ}(u)=u$, from which we deduce the
optimal amplitude $\Delta f_{n}^{\circ}$ :

$\displaystyle\omega_{n}^{2}\frac{\Delta f_{n}^{\circ}/\ell}{\sqrt{1+(\Delta
f_{n}^{\circ}/\ell)^{2}}}-1=0\Rightarrow\frac{\Delta
f_{n}^{\circ}}{\ell}=\frac{1}{\sqrt{\omega_{n}^{4}-1}}$ (24)

and the optimal energy $e_{n}^{\circ}=\sqrt{\omega_{n}^{4}-1}-\omega_{n}^{2}$.
Once again the series in (7) has a logarithmic divergence, and we write in the
large $L/\ell$ limit:

$\displaystyle\Delta\mathcal{E}^{\circ}/[\gamma\ell]$
$\displaystyle\simeq\sum_{{n}=1}^{\infty}\left[e_{n}^{\circ}+\frac{1}{2\widehat{\lambda}{n}}\right]-\frac{1}{2\widehat{\lambda}}\ln(N)+\mathcal{R}_{N}$
(25)

One finds the same type of calculations as before, and the ribbon length
equation (16) writes here

$\displaystyle\frac{{\Lambda}-\sqrt{3}}{\widehat{\lambda}}$
$\displaystyle=\sum_{{n}=1}^{\infty}\frac{\omega_{n}^{2}-\sqrt{\omega_{n}^{4}-1}}{\sqrt{\omega_{n}^{4}-1}}$
(26)

where the series in the rhs is convergent. This equation is the self-
consistent equation setting the value of $\Lambda$ (present in the
$\omega_{n}$ as well). It is important to note (regarding the next subsection)
that as for the $\eta\gg 1$ case,
$\partial_{\Lambda}(\sum_{n}e_{n}^{\circ})=L/\ell-N$, which yields a similar
structure between Eqs. (19) and (26).

The limit $\widehat{\lambda}\rightarrow 0$ is instructive. The right hand side
series is found (using the Euler-MacLaurin formula) equivalent to
$[{\Lambda}-\sqrt{{\Lambda}^{2}-1}]/\widehat{\lambda}$, which gives
${\Lambda}\sim 2$ in this limit. This is coherent in (24) with the value
$\Delta f_{n}^{\circ}/\ell=1/\sqrt{3}$ expected everywhere in this limit
$\widehat{\lambda}=\eta=0$.

####  II.2.2 $\eta\ll 1$ is a singular expansion

To account for the cases where $\eta\ll 1$, but nonzero, the idea is to expand
from the preceding solution. A key point is that the preceding solution is
piecewise linear and therefore not twice differentiable as expected for a
ribbon with a finite bending modulus. As a result, a naive expansion is bound
to fail since the perturbative term should “cure” the zeroth-order singularity
to yield a regular solution, which is just impossible.

To tackle this difficulty, one writes formally $f_{n}(z)=\frac{\Delta
f^{\circ}_{n}}{\ell}z+\xi_{n}(z)$, where $\xi_{n}$ is a small departure from a
zeroth-order which is similar to that computed in section II.2.1, but for the
value $\Lambda$, which is not fixed at the zero-th order by (26) but will be
later on. One makes a Taylor expansion of the energy up to the second order in
(8) and (9), what gives

$e_{n}=e_{n}^{\circ}\\\
+\frac{\eta\ell}{2}\frac{\omega_{n}^{2}(\omega_{n}^{4}-1)^{1/4}}{(\omega_{n}^{4}-1)^{3/4}+1}\left[\int_{0}^{\ell}{\rm
d}z[(\xi_{n}^{\prime\prime})^{2}+k_{n}^{2}(\ell^{-1}\xi_{n}^{\prime})^{2}]\right]$
(27)

where

$\displaystyle k_{n}^{2}$
$\displaystyle=\frac{1}{\eta}\frac{(\omega_{n}^{4}-1)^{5/4}[(\omega_{n}^{4}-1)^{3/4}+1]}{\omega_{n}^{6}}$
(28)

and wherein $\omega_{n}^{2}$ the value of $\Lambda$ is yet to be determined.

At this point, it is crucial to express (27) back in terms of $f_{n}$ and
optimize the energy with respect to $f_{n}$ rather than $\xi_{n}$ (a path also
followed in [22]). The reason is that the expected solution is at least
$C^{1}$ and even twice differentiable, and that the zero-th order of our
expansion is singular with respect to this requirement. As a result, as we
cannot expect to recover with $\tfrac{\Delta f_{n}^{\circ}}{\ell}z+\xi_{n}(z)$
a twice differentiable solution, since $\xi_{n}(z)$ would be twice
differentiable as the result of an optimization process, it is compulsory to
reexpress $\xi_{n}$ in (27) in terms of $f_{n}$ and do the optimization
directly on $f_{n}$. It is interesting to note that no obvious alternate
method seems at hand to solve this issue and produce a first order expansion
in $\eta$ of the solution with the required regularity. This singularity will
be manifest in the fact that the first correction to $\Delta f^{\circ}_{n}$ in
$\Delta f$ will be $O(\sqrt{\eta})$, not $O(\eta)$.

Owing to these remarks, one write (27) making the substitution
$\xi^{\prime}_{n}=f_{n}^{\prime}-\Delta f^{\circ}_{n}/\ell$ and
$\xi_{n}^{\prime\prime}=f_{n}^{\prime\prime}$, to obtain

$e_{n}=e_{n}^{\circ}\\\
+\frac{\eta\ell}{2}\frac{\omega_{n}^{2}(\omega_{n}^{4}-1)^{1/4}}{(\omega_{n}^{4}-1)^{3/4}+1}\left[\int_{0}^{\ell}{\rm
d}z[(f_{n}^{\prime\prime})^{2}+k_{n}^{2}(\ell^{-1}f_{n}^{\prime})^{2}]\right]\\\
+\frac{1}{2\omega_{n}^{4}}\left[\sqrt{\omega_{n}^{4}-1}-2(\omega_{n}^{4}-1)\frac{\Delta
f_{n}}{\ell}\right]$ (29)

Comparison between (29) and (10) shows that $f_{n}(z)=\Delta
f_{n}h_{n}(u=z/\ell)$ where $h_{n}$ is given by the same formula as for
$g_{n}$ in (11) but with $k_{n}$ replacing $\kappa_{n}$. Likewise, the
amplitude $\Delta f_{n}$ minimizes

$e_{n,\rm
opt}-e_{n}^{\circ}=\frac{1}{2}\frac{k_{n}(\omega_{n}^{4}-1)^{3/2}}{\omega_{n}^{4}(k_{n}-2\tanh\frac{k_{n}}{2})}\left(\frac{\Delta
f_{n}}{\ell}\right)^{2}\\\
+\frac{1}{2\omega_{n}^{4}}\left[\sqrt{\omega_{n}^{4}-1}-2(\omega_{n}^{4}-1)\frac{\Delta
f_{n}}{\ell}\right]$ (30)

whence one gets

$\displaystyle\frac{\Delta f_{n,\rm opt}}{\ell}$
$\displaystyle=\left[1-\frac{\tanh\frac{k_{n}}{2}}{k_{n}/2}\right]\left(\frac{1}{\sqrt{\omega_{n}^{4}-1}}\right)$
(31)

and

$\displaystyle e_{n,{\rm
opt}}=e_{n}^{\circ}+\frac{\sqrt{\omega_{n}^{4}-1}}{\omega_{n}^{4}k_{n}}\tanh\frac{k_{n}}{2}$
(32)

which shows that the correction to the energy is consequently a convergent
series. In (31), $\Lambda$ (implicit in $\omega_{n}^{2}$ and $k_{n}$) is given
by taking into account the constraint equation expressing that the length of
the ribbon is $L$ : Up to the first order in
$\xi_{n}^{\prime}=f_{n}^{\prime}-\Delta f_{n}^{\circ}/\ell$, this is

$\mathcal{S}^{\dagger}(\Lambda)=\frac{L}{\ell}-N\\\
=\sum_{n}\left[\frac{\omega_{n}^{2}-\sqrt{\omega_{n}^{4}-1}}{\sqrt{\omega_{n}^{4}-1}}-\frac{1+3\omega_{n}^{4}}{4\omega_{n}^{6}\sqrt{\omega_{n}^{4}-1}}\frac{\tanh\frac{k_{n}}{2}}{k_{n}/2}\right]\\\
+\sum_{n}\frac{\sqrt{\omega_{n}^{4}-1}}{4\omega_{n}^{6}}\left[1-\tanh^{2}\frac{k_{n}}{2}\right]$
(33)

The value of $\Lambda$ is finally found as the one minimizing

$\displaystyle\frac{\Delta\mathcal{E}(\Lambda)}{\gamma\ell}$
$\displaystyle=\sum_{n=1}^{\infty}\left[e_{n,\rm
opt}+\frac{1}{2\widehat{\lambda}n}\right]+\frac{\widehat{\lambda}}{2}{\mathcal{S}^{\dagger}}^{2}+(\sqrt{3}-\Lambda)\mathcal{S}^{{\dagger}}$
(34)

which is simply done using for instance Matlab. It is worth noting that here
the structure of the equation is more complicated than before ($\eta\gg 1$ or
$\eta=0$ cases). This is due to the fact that in (29), the explicit (i.e. not
in $f_{n}$) dependence of the expression with respect to $\Lambda$ is not of
the form $\Lambda(L/\ell-N)$ (as before), but involves also the curvature
energy term. Consequently, we have no longer
$\partial_{\Lambda}\sum_{n}(e_{n,\rm opt}+1/(2\widehat{\lambda}n))$ equating
$\mathcal{S}^{{\dagger}}$, and (34) does not boil down to
$(\Lambda-\sqrt{3})/\widehat{\lambda}=\mathcal{S}^{{\dagger}}$ anymore.

##  III Comparison with experiments

###  III.1 Materials and Methods

A bubble column is produced by bubbling air via a pressure controller through
a nozzle in a detergent solution (composed of water with 4.5 vol% Fairy
liquid, 1.5 vol% Glycerol (Sigma-Aldrich) and 10 g.L-1 J-Lube lubricant
(Jorgensen Labs)). Once the square-cross section tube of width 15 mm is filled
with bubbles, a polydimethylsiloxane thin ribbon (PDMS Sylgard 184, Dow
Corning, 10:1 base to curing agent ratio, cured at 60°C) of width 14.5 $\pm$
0.1 mm, hydrophilised via plasma cleaning treatment on both sides, is inserted
in the central Plateau border. All details of the preparation procedure can be
found in our previous article [22], and are visually highlighted in the video
provided in the Supplementary Information. The surface tension of the solution
is $\gamma=26\pm 1$ mN/m, the thickness of the ribbons are varied using
different spin-coating speeds ($t$=45, 55 and 127 $\unit{\micro}$), and the
elastic material parameters are taken as $E=1.7\pm 0.2$ MPa for the Young’s
modulus, 0.45 for the Poisson’s ratio, and $\rho=1027$ kg.m-3 for the mass
density.

A specific care has been given in the current work to avoid any 3D effect in
the case of long ribbons. To do so, the top of the ribbon has been pierced
with six holes of 2 mm diameter, equally spaced on a horizontal line, to
release lateral contraints at the upper clamping point. We also noticed that
when the bottom of the ribbon arrives right at a node of the bubble column,
there can be a pinning force exerted on the system, so we avoided this case
and considered systems where the end of the ribbon is not at the extremity of
a half-cell.

The structure is allowed to drain for 30 minutes, to avoid motion artefacts,
and is then scanned with a X-ray microtomograph EasyTom 150/160 from RX
Solutions. A helical scan of 48 minutes is performed at resolution 12
$\unit{\micro}$ to capture the deformation of the ribbon. The 3D image is then
converted into 100 slices equally spaced, taken on planes perpendicular to the
ribbon. On each slice, the amplitude of the deformation $\Delta f$ and the
half-cell size $\ell$ is measured, and values are averaged over the 100
slices.

###  III.2 Results and discussion

Figure 2: Decrease of $\Delta f_{n}/\ell$ as a function of $n$ for three
different ribbon thicknesses : $t=45$ $\unit{\micro}$ (red), $t=55$
$\unit{\micro}$ (blue), and $t=127$ $\unit{\micro}$ (yellow). The symbols are
results of measurements from experiments (two different samples for $t=45$
$\unit{\micro}$), the solid line are the result of the theory of section II.2
($\eta\ll 1$) without adjustable parameters, the dashed lines the theory where
the dimensionless mass density $\widehat{\lambda}$ has been adjusted (with a
factor 1.4 and 1.5 for the red and blue curves respectively). The yellow dash-
dotted curve is result of the theory of section II.1 ($\eta\gg 1$). Notice
that for sake of readability, the red curves $y$-axis is located to the right,
whereas the other’s is to the left. For the three experiments, the average
values of $\ell$ are respectively (in mm) 5.7, 6.3 and 6.0 and the associated
$\eta$ are $$1.9\text{\times}{10}^{-2}$,$2.8\text{\times}{10}^{-2}$$ and
$0.38$.

We compared the prediction of the above theory with experiments, and the
results are plotted in Figure 2, where three different ribbons with different
thicknesses were tested. For the thicknesses $t=45$ $\unit{\micro}$ and $t=55$
$\unit{\micro}$, the parameters $\eta$ are respectively
$1.9\text{\times}{10}^{-2}$ and $2.8\text{\times}{10}^{-2}$ and it is relevant
to use the theory of paragraph II.2. Interestingly, the theory predicts
correctly $\Delta f_{1}/\ell$, but underestimates slightly the decrease of
$\Delta f_{n}/\ell$, which is correctly taken into account only if the
dimensionless mass density $\widehat{\lambda}$ is allowed to be renormalized
by a factor $1.4-1.5$. This discrepancy is probably due to the fact that the
ribbon could be coated by 2 thin liquid layers of thickness of order
$t_{w}\sim 10$ $\unit{\micro}$ each. The viscoelastic nature of the detergent
liquid used allows long-lasting structures but also slows down drainage, which
makes this hypothesis likely. In the case of a thicker ribbon
$t=\qty{127}{\micro}$, we used the theory of section II.1 ($\eta\gg 1$), which
captures well the experimental data, as shown in yellow in Figure 2.

##  IV Conclusion

In this work, we have tackled the question of the equilibrium shape of elastic
ribbons in 2D bubble columns, in the case where gravity, capillarity and
elasticity compete to set the equilibrium architecture. We have performed the
analysis in multiple limit cases, and highlighted the most relevant one
(floppy ribbon limit, for which gravity, capillarity and elasticity all
significantly matter in the shape optimisation). This specific problem is
fundamentally interesting as it involves an optimisation around a singular
shape, provided by Plateau’s laws when bubbles only are present. Finally, we
have successfully compared our analysis to experiments with multiple ribbon
thicknesses, providing evidence of the importance of considering gravity in
such systems and assessing the validity of our modeling.

##  V Acknowledgements

We are grateful to Wiebke Drenckhan, Thierry Charitat and François Schosseler
for fruitful discussions, to Damien Favier for his help in some preliminary
tomography experiments, to ICS PECMAT team and especially Fouzia Boulmedais
for the access to a spin-coater tool and to a plasma cleaner, and to
Christophe Lambour for the fabrication of the square-section tube. We also
thank Capucine Loth and Thibaut Schutz for their help in the recording of the
Supplementary Material video. This work of the Interdisciplinary Institute
HiFunMat, as part of the ITI 2021-2028 program of the University of
Strasbourg, CNRS and Inserm, was supported by IdEx Unistra (ANR-10-IDEX-0002)
and SFRI (STRAT’US project, ANR-20-SFRI-0012) under the framework of the
French Investments for the Future Program. We also acknowledge funding from
the IdEx Unistra framework (A. Hourlier-Fargette), and from the ANR (FOAMINT
project, ANR-23-CE06-0014-01).

This research was funded, in whole or in part, by the ANR. A CC-BY public
copyright license has been applied by the authors to the present document and
will be applied to all subsequent versions up to the Author Accepted
Manuscript arising from this submission, in accordance with the grant’s open
access conditions.

##  VI Appendix: Justification of the approximation $\ell$=constant

The volume variations induced by the ribbon are tiny and negligible, and
$\ell={\rm const.}$ is a good approximation, provided the experimental setup
warrants the generation of identical bubbles. To see this, one balances the
deformation work done on a bubble by the typical energy of curvature of a
zigzagging ribbon : $\delta V/V\simeq\alpha/P_{0}V$, where $V$ is a bubble
volume (the geometry of the staircase bubbles is such that all relevant
typical length scales, either direct or transverse, have the same order of
magnitude $\ell$). With typical values for $(\alpha,P_{0},V)$, one gets
$\delta\ell/\ell\sim\delta V/V\sim${10}^{-5}$$, which validates the constant
$\ell$ assumption. Another even more important source of volume variation is
associated to the fact that the ribbon thickness is not entirely negligible
with respect to the tube lateral size : for a ribbon of size
$\sim\qty{100}{\micro}$, the induced (negative) volume variations are a priori
$\delta V/V\sim${10}^{-2}$$. If the pressure stays constant, $\delta V$ must
actually vanish, which would occur via a (positive) variation
$\delta\ell/\ell\sim${10}^{-2}$$ between a HC with a ribbon and a HC without.
If the pressure actually increases slightly, $\delta\ell/\ell$ is even
smaller. All in all, we remark that $\delta\ell/\ell$ is always quite small
and we will henceforth neglect its variations.

## References

  * De-Volder and Hart [2013] M. De-Volder and A. J. Hart, Engineering hierarchical nanostructures by elastocapillary self-assembly, Angewandte Chemie International Edition 52, 2412 (2013).
  * Kwok  _et al._ [2020] K. S. Kwok, Q. Huang, M. Mastrangeli, and D. H. Gracias, Self-folding using capillary forces, Advanced Materials Interfaces 7, 1901677 (2020).
  * Hines  _et al._ [2016] L. Hines, K. Petersen, G. Z. Lum, and M. Sitti, Soft actuators for small-scale robotics, Advanced Materials , 1603483 (2016).
  * Heil  _et al._ [2008] M. Heil, A. L. Hazel, and J. A. Smith, The mechanics of airway closure, Respiratory Physiology & Neurobiology 163, 214 (2008).
  * Gernay  _et al._ [2016] S. Gernay, W. Federle, P. Lambert, and T. Gilet, Elasto-capillarity in insect fibrillar adhesion, Journal of The Royal Society Interface 13, 20160371 (2016).
  * Elettro  _et al._ [2016] H. Elettro, S. Neukirch, F. Vollrath, and A. Antkowiak, In-drop capillary spooling of spider capture thread inspires hybrid fibers with mixed solid–liquid mechanical properties, Proceedings of the National Academy of Sciences 113, 6143 (2016).
  * Style  _et al._ [2017] R. W. Style, A. Jagota, C.-Y. Hui, and E. R. Dufresne, Elastocapillarity: Surface tension and the mechanics of soft solids, Annual Review of Condensed Matter Physics 8, 99 (2017).
  * Holmes [2019] D. P. Holmes, Elasticity and stability of shape-shifting structures, Current Opinion in Colloid & Interface Science 40, 118 (2019).
  * Roman and Bico [2010] B. Roman and J. Bico, Elasto-capillarity: deforming an elastic structure with a liquid droplet, Journal of Physics: Condensed Matter 22, 493101 (2010).
  * Bico  _et al._ [2018] J. Bico, É. Reyssat, and B. Roman, Elastocapillarity: when surface tension deforms elastic solids, Annual Review of Fluid Mechanics 50, 629 (2018).
  * Bico  _et al._ [2004] J. Bico, B. Roman, L. Moulin, and A. Boudaoud, Adhesion: Elastocapillary coalescence in wet hair, Nature 432, 690 (2004).
  * Py  _et al._ [2007] C. Py, P. Reverdy, L. Doppler, J. Bico, B. Roman, and C. N. Baroud, Capillary origami: Spontaneous wrapping of a droplet with an elastic sheet, Phys. Rev. Lett. 98, 156103 (2007).
  * Lau  _et al._ [2003] K. K. S. Lau, J. Bico, K. B. K. Teo, M. Chhowalla, G. A. J. Amaratunga, W. I. Milne, G. H. McKinley, and K. K. Gleason, Superhydrophobic carbon nanotube forests, Nano Letters 3, 1701 (2003).
  * Legrain  _et al._ [2014] A. Legrain, T. G. Janson, J. W. Berenschot, L. Abelmann, and N. R. Tas, Controllable elastocapillary folding of three-dimensional micro-objects by through-wafer filling, Journal of Applied Physics 115, (2014).
  * Miller  _et al._ [2014] J. T. Miller, A. Lazarus, B. Audoly, and P. M. Reis, Shapes of a suspended curly hair, Phys. Rev. Lett. 112, 068103 (2014).
  * Huang  _et al._ [2010] J. Huang, B. Davidovitch, C. D. Santangelo, T. P. Russell, and N. Menon, Smooth cascade of wrinkles at the edge of a floating elastic film, Phys. Rev. Lett. 105, 038302 (2010).
  * Elettro  _et al._ [2018] H. Elettro, A. Antkowiak, and S. Neukirch, The heavy windlass: Buckling and coiling of an elastic rod inside a liquid drop in the presence of gravity, Mechanics Research Communications 93, 58 (2018).
  * Fargette  _et al._ [2014] A. Fargette, S. Neukirch, and A. Antkowiak, Elastocapillary snapping: Capillarity induces snap-through instabilities in small elastic beams, Physical review letters 112, 137802 (2014).
  * Kim and Mahadevan [2006] H.-Y. Kim and L. Mahadevan, Capillary rise between elastic sheets, Journal of Fluid Mechanics 548, 141 (2006).
  * Hutzler  _et al._ [2009] S. Hutzler, J. Barry, P. Grasland-Mongrain, D. Smyth, and D. Weaire, Ordered packings of bubbles in columns of square cross-section, Colloids and Surfaces A: Physicochemical and Engineering Aspects 344, 37 (2009).
  * Reinelt  _et al._ [2001] D. Reinelt, P. Boltenhagen, and N. Rivier, Deformed foam structure and transitions in a tube, The European Physical Journal E 4, 299 (2001).
  * Jouanlanne  _et al._ [2022] M. Jouanlanne, A. Egelé, D. Favier, W. Drenckhan, J. Farago, and A. Hourlier-Fargette, Elastocapillary deformation of thin elastic ribbons in 2d foam columns, Soft Matter 18, 2325 (2022).
  * Shishkov  _et al._ [2022] O. Shishkov, A. Blonder, S. Gowen, T. J. Jones, M. Jouanlanne, Q. Zhang, E. Del Gado, and I. Bischofberger, First annual aps dsoft gallery of soft matter, Physical Review E 106, 050001 (2022).
  * Plateau [1873] J. A. F. Plateau, Statique expérimentale et théorique des liquides soumis aux seules forces moléculaires, Statique expérimentale et théorique des liquides soumis aux seules forces moléculaires, Gauthier-Villars (1873).
  * Note [1] As we consider only half cells (HC), one should impose $g_{n}(0)=1$ and $g_{n}(1)=0$ every two HC. But this point does not change any of the conclusions in the adopted assumption of light ribbons, so we impose $g(0)=0$ and $g(1)=1$ everywhere for sake of simplicity.
  * [26] DLMF, https://dlmf.nist.gov/2.10#i, .


4\sqrt{2}\omega^{2}}\hat{a}^{\dagger}_{1}+{1\over\sqrt{2}}\hat{a}_{2}+{a\over
4\sqrt{2}\omega^{2}}\hat{a}^{\dagger}_{2}\ ,$ $\displaystyle\hat{c}$
$\displaystyle\approx$ $\displaystyle{1\over\sqrt{2}}\hat{a}_{1}-{a\over
4\sqrt{2}\omega^{2}}\hat{a}^{\dagger}_{1}-{1\over\sqrt{2}}\hat{a}_{2}-{a\over
4\sqrt{2}\omega^{2}}\hat{a}^{\dagger}_{2}\ ,$ $\displaystyle\hat{b}^{\dagger}$
$\displaystyle\approx$
$\displaystyle{1\over\sqrt{2}}\hat{a}^{\dagger}_{1}-{a\over
4\sqrt{2}\omega^{2}}\hat{a}_{1}+{1\over\sqrt{2}}\hat{a}^{\dagger}_{2}+{a\over
4\sqrt{2}\omega^{2}}\hat{a}_{2}\ ,$ $\displaystyle\hat{c}^{\dagger}$
$\displaystyle\approx$
$\displaystyle{1\over\sqrt{2}}\hat{a}^{\dagger}_{1}-{a\over
4\sqrt{2}\omega^{2}}\hat{a}_{1}-{1\over\sqrt{2}}\hat{a}^{\dagger}_{2}-{a\over
4\sqrt{2}\omega^{2}}\hat{a}_{2}\ .$ (B.24)

Now we note that the combination

$\hat{b}+{a\over 4\omega^{2}}\hat{c}^{\dagger}\ ,$ (B.25)

has only annihilation operators $\hat{a}_{1}$ and $\hat{a}_{2}$. Thus

$\left(\hat{b}+{a\over 4\omega^{2}}\hat{c}^{\dagger}\right)|0\rangle_{a}=0\ ,$
(B.26)

which has the solution

$|0\rangle_{a}=Ce^{-{a\over
4\omega^{2}}\hat{b}^{\dagger}\hat{c}^{\dagger}}|0\rangle_{b}\otimes|0\rangle_{c}\
.$ (B.27)

We see that if we have two oscillators with the same frequency, weakly coupled
together and then we remove the coupling suddenly, the ground state of the
initial system becomes an entangled state of the two uncoupled oscillators.

###  B.5 The entangled nature of the final state

We can now see the entangled nature of the state (B.27). We can expand the
exponential in (B.27) to find

$|0\rangle_{a}=C\left[\leavevmode\nobreak\
|0\rangle_{b}\otimes|0\rangle_{c}-\left({a\over
4\omega^{2}}\right)|1\rangle_{b}\otimes|1\rangle_{c}+\left({a\over
4\omega^{2}}\right)^{2}|2\rangle_{b}\otimes|2\rangle_{c}+\dots\right]\ .$
(B.28)

Thus the above model with two oscillators gives a bit model for Hawking
radiation. When the wavelength of the mode is small compared to the horizon
radius (say $\lambda\sim r_{h}/10$), the parts of the mode that are just
inside and just outside the horizon are strongly coupled and such coupled
modes are described by the oscillators (B.16). When the wavelength of the mode
is large ($\lambda\gtrsim r_{h}$) then the inside and outside parts are weakly
coupled and are described by the mode expansion (B.20). The relation (B.21)
relates these two mode expansions and encodes the essence of the Hawking pair
creation process.

## References

  * [1] S. W. Hawking, “Particle Creation by Black Holes,” Commun. Math. Phys. 43 (1975) 199–220. [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * [2] S. W. Hawking, “Breakdown of Predictability in Gravitational Collapse,” Phys. Rev. D 14 (1976) 2460–2473.
  * [3] O. Lunin and S. D. Mathur, “AdS / CFT duality and the black hole information paradox,” Nucl. Phys. B 623 (2002) 342–394, arXiv:hep-th/0109154.
  * [4] S. D. Mathur, A. Saxena, and Y. K. Srivastava, “Constructing ‘hair’ for the three charge hole,” Nucl. Phys. B 680 (2004) 415–449, arXiv:hep-th/0311092.
  * [5] O. Lunin, J. M. Maldacena, and L. Maoz, “Gravity solutions for the D1-D5 system with angular momentum,” arXiv:hep-th/0212210.
  * [6] I. Kanitscheider, K. Skenderis, and M. Taylor, “Fuzzballs with internal excitations,” JHEP 06 (2007) 056, arXiv:0704.0690 [hep-th].
  * [7] S. D. Mathur, “The Fuzzball proposal for black holes: An Elementary review,” Fortsch. Phys. 53 (2005) 793–827, arXiv:hep-th/0502050.
  * [8] I. Bena and N. P. Warner, “Black holes, black rings and their microstates,” Lect. Notes Phys. 755 (2008) 1–92, arXiv:hep-th/0701216.
  * [9] B. D. Chowdhury and A. Virmani, “Modave Lectures on Fuzzballs and Emission from the D1-D5 System,” in 5th Modave Summer School in Mathematical Physics.  1, 2010.  arXiv:1001.1444 [hep-th].
  * [10] I. Bena, S. Giusto, R. Russo, M. Shigemori, and N. P. Warner, “Habemus Superstratum! A constructive proof of the existence of superstrata,” JHEP 05 (2015) 110, arXiv:1503.01463 [hep-th].
  * [11] I. Bena, E. Martinec, D. Turton, and N. P. Warner, “Momentum Fractionation on Superstrata,” JHEP 05 (2016) 064, arXiv:1601.05805 [hep-th].
  * [12] I. Bena, S. Giusto, E. J. Martinec, R. Russo, M. Shigemori, D. Turton, and N. P. Warner, “Smooth horizonless geometries deep inside the black-hole regime,” Phys. Rev. Lett. 117 no. 20, (2016) 201601, arXiv:1607.03908 [hep-th].
  * [13] J. Maldacena and L. Susskind, “Cool horizons for entangled black holes,” Fortsch. Phys. 61 (2013) 781–811, arXiv:1306.0533 [hep-th].
  * [14] S. D. Mathur, “The Information paradox: A Pedagogical introduction,” Class. Quant. Grav. 26 (2009) 224001, arXiv:0909.1038 [hep-th].
  * [15] G. W. Gibbons and N. P. Warner, “Global structure of five-dimensional fuzzballs,” Class. Quant. Grav. 31 (2014) 025016, arXiv:1305.0957 [hep-th].
  * [16] W. Z. Chua and N. Afshordi, “Electromagnetic Albedo of Quantum Black Holes,” JHEP 07 (2021) 006, arXiv:2103.05790 [hep-th].
  * [17] [See for example the talk and discussion by Maldacena at the workshop ‘Black-hole Microstructure’, June 7-11, 2021 (Paris) http://www.youtube.com/watch?v=-vHrjKGgtMk].
  * [18] C. G. Callan, Jr., S. B. Giddings, J. A. Harvey, and A. Strominger, “Evanescent black holes,” Phys. Rev. D 45 no. 4, (1992) R1005, arXiv:hep-th/9111056.
  * [19] J. G. Russo, L. Susskind, and L. Thorlacius, “Black hole evaporation in (1+1)-dimensions,” Phys. Lett. B 292 (1992) 13–18, arXiv:hep-th/9201074.
  * [20] E. Keski-Vakkuri and S. D. Mathur, “Evaporating black holes and entropy,” Phys. Rev. D 50 (1994) 917–929, arXiv:hep-th/9312194.
  * [21] T. M. Fiola, J. Preskill, A. Strominger, and S. P. Trivedi, “Black hole thermodynamics and information loss in two-dimensions,” Phys. Rev. D 50 (1994) 3987–4014, arXiv:hep-th/9403137.
  * [22] S. D. Mathur, “Resolving the black hole causality paradox,” Gen. Rel. Grav. 51 no. 2, (2019) 24, arXiv:1703.03042 [hep-th].
  * [23] S. D. Mathur, “Spacetime has a “thickness”,” Int. J. Mod. Phys. D 26 no. 12, (2017) 1742002, arXiv:1705.06407 [hep-th].
  * [24] G. Penington, S. H. Shenker, D. Stanford, and Z. Yang, “Replica wormholes and the black hole interior,” arXiv:1911.11977 [hep-th].
  * [25] A. Strominger and C. Vafa, “Microscopic origin of the Bekenstein-Hawking entropy,” Phys. Lett. B 379 (1996) 99–104, arXiv:hep-th/9601029.
  * [26] M. E. Agishtein, L. Jacobs, A. A. Migdal, and J. L. Richardson, “Geometric Characterization of States in Two-dimensional Quantum Gravity,” Mod. Phys. Lett. A 5 (1990) 965.
  * [27] S. Jain and S. D. Mathur, “World sheet geometry and baby universes in 2-D quantum gravity,” Phys. Lett. B 286 (1992) 239–246, arXiv:hep-th/9204017.
  * [28] V. G. Knizhnik, A. M. Polyakov, and A. B. Zamolodchikov, “Fractal Structure of 2D Quantum Gravity,” Mod. Phys. Lett. A 3 (1988) 819.
  * [29] F. David, “Conformal Field Theories Coupled to 2D Gravity in the Conformal Gauge,” Mod. Phys. Lett. A 3 (1988) 1651.
  * [30] J. Distler and H. Kawai, “Conformal Field Theory and 2D Quantum Gravity,” Nucl. Phys. B 321 (1989) 509–527.
  * [31] A. Karlsson, “Replica wormhole and island incompatibility with monogamy of entanglement,” arXiv:2007.10523 [hep-th].
  * [32] A. Karlsson, “Concerns about the replica wormhole derivation of the island conjecture,” arXiv:2101.05879 [hep-th].
  * [33] A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian, and A. Tajdini, “The entropy of Hawking radiation,” Rev. Mod. Phys. 93 no. 3, (2021) 035002, arXiv:2006.06872 [hep-th].
  * [34] D. N. Page, “Average entropy of a subsystem,” Phys. Rev. Lett. 71 (1993) 1291–1294, arXiv:gr-qc/9305007.
  * [35] H. Liu and S. Vardhan, “Entanglement entropies of equilibrated pure states in quantum many-body systems and gravity,” P. R. X. Quantum. 2 (2021) 010344, arXiv:2008.01089 [hep-th].
  * [36] S. D. Mathur, “Fuzzballs and black hole thermodynamics,” arXiv:1401.4097 [hep-th].
  * [37] S. W. Hawking and W. Israel, General Relativity: An Einstein Centenary Survey.  Univ. Pr., Cambridge, UK, 1979.
  * [38] E. Witten, “Anti-de Sitter space, thermal phase transition, and confinement in gauge theories,” Adv. Theor. Math. Phys. 2 (1998) 505–532, arXiv:hep-th/9803131.
  * [39] [See for example comments by Maldacena in the discussion session chaired by Engelhardt and Myers at the conference ‘Strings 2021’, June 21- July 2, 2021 (São Paulo) https://www.youtube.com/watch?v=hJjrbX_7WJY].
  * [40] D. Marolf and H. Maxfield, “Transcending the ensemble: baby universes, spacetime wormholes, and the order and disorder of black hole information,” JHEP 08 (2020) 044, arXiv:2002.08950 [hep-th].
  * [41] D. Marolf and H. Maxfield, “Observations of Hawking radiation: the Page curve and baby universes,” JHEP 04 (2021) 272, arXiv:2010.06602 [hep-th].
  * [42] K. Papadodimas and S. Raju, “An Infalling Observer in AdS/CFT,” JHEP 10 (2013) 212, arXiv:1211.6767 [hep-th].
  * [43] S. Raju, “Failure of the split property in gravity and the information paradox,” arXiv:2110.05470 [hep-th].
  * [44] P. Gao, D. L. Jafferis, and A. C. Wall, “Traversable Wormholes via a Double Trace Deformation,” JHEP 12 (2017) 151, arXiv:1608.05687 [hep-th].
  * [45] J. Maldacena, D. Stanford, and Z. Yang, “Diving into traversable wormholes,” Fortsch. Phys. 65 no. 5, (2017) 1700034, arXiv:1704.05333 [hep-th].
  * [46] S. D. Mathur, “What prevents gravitational collapse in string theory?,” Int. J. Mod. Phys. D 25 no. 12, (2016) 1644018, arXiv:1609.05222 [hep-th].
  * [47] O. Lunin and S. D. Mathur, “Rotating deformations of AdS(3) x S**3, the orbifold CFT and strings in the pp wave limit,” Nucl. Phys. B 642 (2002) 91–113, arXiv:hep-th/0206107.
  * [48] J. Gomis, L. Motl, and A. Strominger, “PP wave / CFT(2) duality,” JHEP 11 (2002) 016, arXiv:hep-th/0206166.
  * [49] E. Gava and K. S. Narain, “Proving the PP wave / CFT(2) duality,” JHEP 12 (2002) 023, arXiv:hep-th/0208081.
  * [50] S. Hampton, S. D. Mathur, and I. G. Zadeh, “Lifting of D1-D5-P states,” JHEP 01 (2019) 075, arXiv:1804.10097 [hep-th].
  * [51] S. D. Mathur and D. Turton, “The flaw in the firewall argument,” Nucl. Phys. B 884 (2014) 566–611, arXiv:1306.5488 [hep-th].
  * [52] A. Almheiri, D. Marolf, J. Polchinski, and J. Sully, “Black Holes: Complementarity or Firewalls?,” JHEP 02 (2013) 062, arXiv:1207.3123 [hep-th].
  * [53] S. D. Mathur, “A model with no firewall,” arXiv:1506.04342 [hep-th].
  * [54] S. D. Mathur, “The VECRO hypothesis,” arXiv:2001.11057 [hep-th].


# A Distributed Sparse Channel Estimation Technique for mmWave Massive MIMO
Systems

Maria Trigka, Christos Mavrokefalidis, Kostas Berberidis
Department of Computer Engineering and Informatics, University of Patras,
Greece
<EMAIL_ADDRESS>
###### Abstract

In this paper, we study the problem of sparse channel estimation via a
collaborative and fully distributed approach. The estimation problem is
formulated in the angular domain by exploiting the spatially common sparsity
structure of the involved channels in a multi-user scenario. The sparse
channel estimation problem is solved via an efficient distributed approach in
which the participating users collaboratively estimate their channel sparsity
support sets, before locally estimate the channel values, under the assumption
that global and common support subsets are present. The performance of the
proposed algorithm, named WDiOMP, is compared to DiOMP, local OMP and a
centralized solution based on SOMP, in terms of the support set recovery error
under various experimental scenarios. The efficacy of WDiOMP is demonstrated
even in the case in which the underlining sparsity structure is unknown.

##  I Introduction

In 5G (and beyond) wireless communication systems, the millimeter wave
(mmWave) spectrum is considered for increasing, among others, the transmission
capacity [1],[2],[3]. Signals in the mmWave bands, however, experience high
sensitivity to path loss, limiting the communication range. To mitigate this
loss, massive Multi-Input, Multi-Output (MIMO) and beamforming have been
suggested in order to guarantee reliable communication. To benefit from the
massive MIMO array gain and perform high directional beamforming, channel
information is required between all antenna pairs, increasing the training
overhead. On the other hand, due to limited scattering [4], a small number of
significant paths actually contribute to the transmission. Thus, in the
angular domain, channel estimation can be achieved with reduced training
overhead by identifying the channel gains, the Angles of Arrival (AoA) and
Angles of Departure (AoD), related only to those paths via Compressed Sensing
(CS) [5].

Exploiting CS for channel estimation has been an active research area in
recent years [6]. Most of the proposed works employ CS algorithms either at
individual nodes [5] or at a central point (namely, a Base Station - BS), in
which multiple measurements of groups of nodes are fed back [4], [7]. In the
latter case, the CS algorithms are able to exploit common sparsity patterns
that manifest in the node measurements due to the transmission environment
and, thus, achieve improved estimation performance or reduced training
overhead. In this paper, a fully distributed channel estimation algorithm is
proposed that (a) exploits common sparsity patterns among the collaborating
nodes for improving performance and (b) requires no central point for
gathering measurements and performing the associated processing. The
deployment of the proposed algorithm is able to reduce the channel uses for
training by the BS and improve the channel estimation performance of the
nodes, compared to the case they operate as individual entities (especially,
at the low SNR regime).

In more detail, in this paper, it is assumed that the collaborating nodes
receive, in the downlink, a number of paths from identical angles while, in
the remaining ones, the angles might be different [4]. This means that the
involved channels have a global and, possibly, common sparsity patterns (or
support sets), observed by all and groups of nodes, respectively. Capitalizing
on the distributed sparse coding literature, where the previous “sparsity
model” was introduced [8] and solved [9], [10], the proposed technique extends
the Distributed Orthogonal Matching Pursuit (DiOMP) algorithm [10] by
employing a weighted majority voting mechanism with adaptive combination
weights to identify the desired support sets. The voting mechanism, contrary
to the simple majority voting rule used in [10], assumes no previous knowledge
about the structure of the global and common sparsity patterns (which in the
distributed literature answers also to the term “node-specific” [11]), apart
from the total number of involved paths. The performance of the proposed
algorithm is extensively assessed via simulations in three different scenarios
and it turns out that the proposed scheme exhibits considerable improvement
over the existing ones.

The rest of this paper is organized as follows. In Section II, the data model
is described. In Section III, the proposed distributed MIMO downlink channel
modeling and estimation, with the aid of the CS theory, is elaborated in
detail. Finally, in Section IV, several experiments are presented to verify
the efficiency of the proposed approach.

##  II System Model

In this section, a brief description of the adopted data model is provided.
Consider a massive MIMO system with $K$ single antenna users and a BS equipped
with $N$ antennas which broadcasts a sequence of $T$ training pilot vectors,
denoted by $\boldsymbol{X}\in\mathbb{C}^{T\times N}$, towards each user for
estimating the downlink channel. Then, the downlink received signal
$\boldsymbol{y}_{k}\in\mathbb{C}^{T\times 1}$ at the $k$-th user is given by

$\boldsymbol{y}_{k}=\boldsymbol{X}\boldsymbol{h}_{k}+\boldsymbol{n}_{k},$ (1)

where $\boldsymbol{n}_{k}\in\mathbb{C}^{\textrm{T}\times 1}$ stands for
additive noise which is modeled as random vector with elements being
independent and identically distributed as complex Gaussian random variables
with zero mean and variance equal to $\sigma^{2}$.

Assuming a block, flat fading channel
$\boldsymbol{h}_{k}\in\mathbb{C}^{N\times 1}$ from the BS to the $k$-th user,
the Geometry-Based Stochastic Channel Model (GSCM) [12] is adopted which,
after dropping the index $k$ for clarity, can be written as

$\boldsymbol{h}=\sum_{c=1}^{N_{c}}\sum_{s=1}^{N_{c,s}}\gamma_{c,s}\boldsymbol{a}(\theta_{c,s}),$
(2)

where $N_{c}$ and $N_{c,s}$ stand for the number of scattering clusters and
the number of sub-paths per scattering cluster, respectively, $\gamma_{c,s}$
is the complex gain of the $s$-th sub-path in the $c$-th scattering cluster,
$\theta_{c,s}$ is the corresponding AoD and

$\boldsymbol{a}(\theta_{c,s})=\frac{1}{\sqrt{N}}\left[1,\,e^{j\frac{2d\pi}{\lambda_{c}}\sin{\theta_{c,s}}},\dots,e^{j\frac{2(N-1)d\pi}{\lambda_{c}}\sin{\theta_{c,s}}}\right]^{\textrm{T}}$
(3)

is a steering vector, assuming a Uniform Linear Array (ULA). From now on, we
denote the true AoDs as $\theta_{l},l=1,2,\ldots,L$, with $L=N_{c}N_{c,s}$.

Exploiting the limited number of scattering clusters in the BS side and few
active transmit directions for each user, the channel can be approximated by
$\boldsymbol{h}\approx\boldsymbol{\Psi}\boldsymbol{w}$, where $\boldsymbol{w}$
is a sparse vector and $\boldsymbol{\Psi}\in\mathbb{C}^{N\times\hat{L}}$ is a
sparsifying dictionary [12]. A typical $\boldsymbol{\Psi}$ is the normalized
square DFT matrix, leading to the known “virtual channel model”, which maps
the spatial channel response to the angular domain. Alternatively, an
overcomplete dictionary $\boldsymbol{\Psi}$ can be employed via the
overcomplete DFT matrix, where $\hat{L}\gg N$ is the number of atoms in
$\boldsymbol{\Psi}$ corresponding to the AoDs in a finer grid. In this case,
$\boldsymbol{w}\in\mathbb{C}^{\hat{L}\times 1}$ is able to approximate more
closely the true angles at $\theta^{l}$. Thus, $\boldsymbol{y}_{k}$ can be
written as

$\displaystyle\boldsymbol{y}_{k}\approx\boldsymbol{X}\boldsymbol{\Psi}\boldsymbol{w}_{k}+\boldsymbol{n}_{k},$
(4)

where $\boldsymbol{X}\boldsymbol{\Psi}$ is the measurement matrix.To achieve
the CS recovery, the measurement matrix $\boldsymbol{X}\boldsymbol{\Psi}$
should satisfy specific conditions like the ones determined by the Restricted
Isometry Property (RIP). For example, this can be achieved if the training
matrix $\boldsymbol{X}$ has random elements which are identically and
independently distributed [12].

##  III Cooperative and distributed sparse channel estimation

In this section, the problem of mmWave MIMO sparse channel estimation is
discussed. First, the adopted, joint sparsity model for the involved channels
will be presented. Then, the proposed cooperative and distributed algorithm
will be described.

###  III-A Sparse Channel Modeling

Let us consider a set $\mathcal{V}=\\{1,2,\ldots,K\\}$ of $K$ geographically
distributed users with each one cooperating with all its connected neighbors.
Due to the scattering nature of the environment (see Fig. 1), the users belong
to groups depending on which scatterers are involved to their channels.
Following [4], the users of a group face the same scattering structure because
they are affected by the same scatterers, which means that the involved sparse
vectors $\boldsymbol{w}$’s have the same sparsity support sets. However, the
values of the coefficients at the corresponding positions may be different,
due to the surrounding environment of each user. Between two groups, the
scattering structure might be similar, namely, part of the sparsity support
sets of the involved $\boldsymbol{w}$’s can be the same.

The performance of the downlink sparse channel estimation can be improved if
we exploit the aforementioned relation on the sparsity support sets among the
involved channels of the users [4], [7]. In more detail, the transmission
environment at the BS side, as illustrated in Fig.1, can be organized into
global and common (denoted as $g$ and $c_{j}$, respectively) scattering
clusters, leading to respective AoDs that are global to all users, while
others are only common to particular groups. For each user $k$, there is a
subset $\boldsymbol{I}_{k}$ of indices $j=\\{1,2,\ldots,J\\}$ indicating the
corresponding common interest support sets and groups in the network. This
situation translates to the respective atoms of the employed dictionary
$\boldsymbol{\Psi}$ which can be distinguished as global (to all users) and
common (to groups of users), as well.

Adopting the aforementioned sparsity model, the sparse representation vector
$\boldsymbol{w}_{k}$ of user $k$ can be formulated as

$\boldsymbol{w}_{k}=\boldsymbol{w}_{k}^{g}+\sum_{j\in
I_{k}}\boldsymbol{w}_{k}^{c_{j}},$ (5)

where $\boldsymbol{w}_{k}^{g}$ stands for the sparse representation vector
whose support set is globally shared among all users in the network and
$\boldsymbol{w}_{k}^{c_{j}}$ stands for the sparse representation vector whose
support set is commonly shared among users in a particular group in the
network. Hence, a Global and Multiple Common support sets model arises, which
captures the sparse properties of the involved channels among the BS and each
user in the network.

Moreover, it is assumed that $\boldsymbol{w}_{k}^{g}$ of each user $k$
consists of $L_{g}$ non-zero parameters and $\boldsymbol{w}_{k}^{c_{j}}$’s
consist of $L_{c_{j}}$ non-zero parameters. Hence, the sparsity profile of a
user $k$ consists of global and multiple common interest support subsets. In
the following, some definitions will be provided. In the following, some
definitions will be provided. There, $\Omega=\\{1,2,\ldots,\hat{L}\\}$ denotes
the set of atom indices.

###### Definition III.1

(Global Interest Support Set): Let the globally shared sparse representation
vector $\boldsymbol{w}^{g}_{k}$ with only $L_{g}$ nonzero entries. The global
interest support set is defined as

$\mathcal{S}^{g}=supp(\boldsymbol{w}^{g}_{m})=supp(\boldsymbol{w}^{g}_{n}),\forall
m,n\in\mathcal{V},$ (6)

where $\mathcal{S}^{g}\subset\Omega$ and $|\mathcal{S}^{g}|=L_{g}$ ($l_{0}$
sparsity).

###### Definition III.2

(Common Interest Support Set): Let the sparse representation vector
$\boldsymbol{w}^{c_{j}}_{k}$ with only $L_{c_{j}}$ nonzero entries. The common
interest support set is defined as

$\mathcal{S}^{c_{j}}=supp(\boldsymbol{w}^{c_{j}}_{m})=supp(\boldsymbol{w}^{c_{j}}_{n}),\forall
m,n\in\mathcal{C}_{j},$ (7)

where $\mathcal{S}^{c_{j}}\subset\Omega$, $|\mathcal{S}^{c_{j}}|=L_{c_{j}}$
($l_{0}$ sparsity) and $\mathcal{C}_{j}$ a set of users indices that are
impacted by the common scattering cluster $c_{j}$, for $j=1,2,\ldots,J$, and
are concerned about $\mathcal{S}^{c_{j}}$.

Therefore, the support set of a user $k$ can be formulated as
$\mathcal{S}_{k}=\mathcal{S}_{k}^{g}\cup\\{\mathcal{S}_{k}^{c_{j}}\\}$, for
all $j\in I_{k}$, meaning that the joint sparsity $L\leq L_{g}+\sum_{j\in
I_{k}}L_{c_{j}}$. Here, although $\mathcal{S}^{g}$ is the same for all users
in the network or $\mathcal{S}^{c_{j}}$ is similar among users in the same
sub-group $\mathcal{C}_{j}$, the corresponding non-zero values of
$\boldsymbol{w}_{k}^{g}$ and $\boldsymbol{w}_{k}^{c_{j}}$ are still individual
and possibly independent among the users.

Figure 1: Downlink transmission scenario: $\mathcal{C}_{1}=\\{1,2,3,4,5,6\\}$,
$\mathcal{C}_{2}=\\{5,6,7,8,9,10\\}$.

In summary, the channel $\boldsymbol{h}_{k}$ of user $k$ is written as

$\boldsymbol{h}_{k}=\boldsymbol{\Psi}(\boldsymbol{w}_{k}^{g}+\sum_{j\in
I_{k}}\boldsymbol{w}_{k}^{c_{j}}).$ (8)

Considering (8), our aim is to exploit the joint sparsity of the support sets
of neighboring users to estimate the per user sparse channel in a cooperative
and distributed way. Such an approach will be elaborated in the following.

###  III-B Proposed cooperative and distributed algorithm

The problem of estimating $\boldsymbol{h}_{k}$ in a collaborative and
distributed manner, using the predefined dictionary $\boldsymbol{\Psi}$, is
cast to the estimation of $\boldsymbol{w}_{k}$. The algorithm estimates
collaboratively the support sets and locally, via Least-Squares (LS), the non-
zero values at each user.

In more detail, the algorithm is based on DiOMP [10] which is improved in a
key aspect by considering a new, weighted majority voting mechanism with
adaptive weights, instead of the simple majority voting that DiOMP considers.
In particular, via the new mechanism, the proposed algorithm does not require
any knowledge about the global and common sparsity structure of the involved
channels, apart from the total number of non-zero AoDs. Also, it is noted
that, to the best of the authors’ knowledge, distributed sparse coding
approaches (such as DiOMP) are for the first time employed for a channel
estimation task. The proposed algorithm, called WDiOMP, consists of an
initialization stage, an iterative stage, where the users collaborate to
estimate the sparsity support sets, and, finally, a local channel estimation
stage.

### Algorithm WDiOMP

First stage:

  1. 1.

For $T$ time slots, the BS transmits (broadcasts) a number of pilots,
constituting matrix $\boldsymbol{X}$, to all users in the network.

  2. 2.

Each user utilizes the collected measurements $\boldsymbol{y}_{k}$ and
independently gets an initial estimation of its complete support set
$\hat{\mathcal{S}}_{k}$, performing the standard OMP, initialized with an
empty support set.

Second stage (iterates over the following steps for $i=1$ to $L$, namely, up
to the total sparsity level):

  1. 1.

User $k$, $\forall k$, transmits its $\hat{\mathcal{S}}_{k}$ to its
neighboring users (using an appropriate D2D communication protocol [13]),
denoted by the set $\mathcal{L}_{k}^{out}$, and receives the estimates
$\hat{\mathcal{S}}_{l}$’s from its neighbors, $\forall
l\in\mathcal{L}_{k}^{in}$.

  2. 2.

Each user utilizes the received $\hat{\mathcal{S}}_{l}$’s and applies a
weighted majority voting mechanism to select the $i$ best indices.

  3. 3.

The remaining non-zero indices, up to $L$, are acquired locally by each user
by employing OMP.

Third stage:
Each user utilizes the estimated set $\hat{\mathcal{S}}_{k}$ on (4) by keeping
only the relevant columns of $\boldsymbol{X\Psi}$ and finds the non-zero
values by employing LS.

The motivation for the weighted voting mechanism in the second step of the
second stage of the Algorithm WDiOMP stems from the fact that a user may not
contribute in the same manner to the support set recovery for several reasons
like (a) different communication conditions (noise level), (b) different
common support sets among groups and (c) lack of knowledge concerning the
global and common parts of the support sets and the grouping of users. To cope
with these cases, a weighted majority voting mechanism is employed which
combines the votes of the cooperating users with a weight depicting the
suitability of a user participating in a group. In more detail, at iteration
$i$, user $k$ utilizes the vector $\boldsymbol{z}_{k,i}$ of size $\hat{L}$,
initially, with zero elements. For each received
$\\{\hat{\mathcal{S}}_{l,i}\\}$, $l\in\mathcal{L}_{k}^{in}$, user $k$ adds the
weight $a_{lk}$ to the positions of $\boldsymbol{z}_{k,i}$, indexed by
$\hat{\mathcal{S}}_{l,i}$. User $k$, incorporating the information from all
users $l\in\mathcal{L}_{k}^{in}$, it subsequently selects the indexes of the
$i$ largest elements of the final $\boldsymbol{z}_{k,i}$. The weight $a_{lk}$
is updated with a mechanism similar to the one in [14] as follows.

$\displaystyle b_{lk,i}$
$\displaystyle=(1-v)b_{lk,i-1}+v|\hat{\mathcal{S}}_{l,i}\setminus\hat{\mathcal{S}}_{k,i-1}|.$
(9) $\displaystyle a_{lk,i}$
$\displaystyle=\frac{\frac{1}{b_{lk,i}}}{\sum_{m\in\mathcal{L}_{k}^{in}}\frac{1}{b_{mk,i}}},l\in\mathcal{L}_{k}^{in},k\in
K.$ (10)

In (9), the second term captures the number of different elements in the
involved sets and $v\in(0,1)$ is a small positive factor. It is noted at this
point that, for a sufficient number of measurements and high enough SNR
values, when only global support set exists, it has been observed
experimentally that the weights tend to
$a_{lk}=\frac{1}{|\mathcal{L}_{k}^{in}|}$ for all users, while when both a
global and common support sets exist,
$a_{lk}=\frac{1}{|\mathcal{L}_{k}^{in}\cap\mathcal{C}_{j}|}$ for users in the
same group and for users outside the group may tend to zero.

Before concluding this section, the transmission efficiency in bits of the
cooperative schemes, namely DiOMP and WDiOMP, is compared to the case when the
measurements are gathered in a central point (for processing by a centralized
algorithm like SOMP [15]). Using $q$ bits for the real and imaginary parts of
the measurements in $\boldsymbol{y}$ for feeding them back to the BS, $2\cdot
K\cdot T\cdot q$ \+ $L\cdot\lceil{log_{2}(\hat{L})}\rceil$ bits are required,
where the second term captures the transmission overhead of the recovered
support sets indices back to the users. For the (W)DiOMP-based schemes,
$K\cdot|\mathcal{L}_{k}^{out}|\cdot\lceil{log_{2}(\hat{L})}\rceil\cdot L^{2}$
bits are required during cooperation. For example, if $q=36$bits, $T=20$,
$K=10$, $\hat{L}=200$, $L=5$, $|\mathcal{L}_{k}^{out}|=6$ $\forall
k\in\mathcal{V}$, SOMP demands $14.4$ Kbits while (W)DiOMP-based schemes
require $12$Kbits.

##  IV Simulations

In this section, the performance of WDiOMP will be presented under three
experiments. $K=10$ single antenna users are considered along with a BS with
$N=128$ antennas. An overcomplete, predefined dictionary $\boldsymbol{\Psi}$
(using DFT) is employed with $\hat{L}=200$ atoms, which define a corresponding
grid for $[-\pi,\pi]$, while $\theta_{l}\in[-\pi,\pi]$, and $T\in[10,40]$. The
elements of the training matrix $\boldsymbol{X}$ are i.i.d. random variables
that follow $\mathcal{CN}(0,\frac{1}{T})$ (as in [16]). The $L$ none-zero
elements in $\boldsymbol{w}_{k}$ are also i.i.d. random variables following
$\mathcal{N}(0,1)$ (as in [9]). Each user has knowledge of only the total
sparsity level $L$ (i.e., not the individual parameters $L_{g}$, $L_{c}$).
Finally, it is assumed that all users are connected via a network and can
share the positions of the non-zero elements in their $\boldsymbol{w}_{k}$’s,
namely, their support sets. Thus,
$\mathcal{L}_{k}^{in}=\mathcal{L}_{k}^{out}=\mathcal{V}$ holds.

The Average Support set Cardinality Error (ASCE) is considered for the
evaluation. It takes values in the range $[0,1]$ and is defined as follows

$ASCE=1-\frac{1}{K}\frac{1}{M}\sum_{m=1}^{M}\sum_{k=1}^{K}\bigg{\\{}\frac{|\hat{\mathcal{S}}_{k}^{m}\cap\mathcal{S}_{k}|}{|\mathcal{S}_{k}|}\bigg{\\}}.$
(11)

For the results, $M=1000$ Monte Carlo simulations were performed.

In the following, the performance of WDiOMP (using $v=0.1$ for the calculation
of the weights) will be compared with DiOMP (which employs a simple majority
voting), the centralized algorithm SOMP [15] and the per user (locally)
employed OMP algorithm. Three scenarios will be presented and the $ASCE$ of
the various approaches will be assessed versus the number of measurements $T$.

In the first experiment, depicted in Fig. 2, $ASCE$ is evaluated assuming that
all users share a global support set and considering two $SNR$ values, namely,
$10$ and $20$ dB. It is observed that both the DiOMP and WDiOMP achieve a
performance close to the one achieved by SOMP, which, being the point of
reference for this experiment, was expected to demonstrate the best, though
similar, performance. Additionally, the previous schemes perform considerably
better than the per-user (locally) employed OMP algorithm, resulting, for
instance, in about $40\%$ reduction in required training measurements for
achieving near zero ASCE, for $SNR=20$ dB. Moreover, for $SNR=10$ dB, the
reduction is even more pronounced as the OMP approach does not reach a near
zero ASCE for the $T$ values considered. This experiment demonstrates the
benefit of the users when they participate in a collaborative channel
estimation procedure that exploits the adopted sparsity model, as opposed to
the case in which each user operates in isolation. Finally, it should be
mentioned that although DiOMP, WDiOMP and SOMP perform similarly, the first
two do not utilize the BS during channel estimation.

In the second experiment, shown in Fig. 3, an “unbalanced” scenario is
examined in which, although the users share only a global support set, they
operate in different SNR levels, namely, $SNR=20$ dB and $SNR=0$ dB for 7 and
3 users, respectively. For the OMP case, in which the users operate at
isolation, apart from the overall mean performance (depicted by the “OMP”
curve), the performances when focusing on the groups that face $SNR=20$ dB and
$SNR=0$ dB, are also presented. It is observed that the collaboration of the
users (either centrally via SOMP or in a distributed fashion via DiOMP and
WiOMP) improves considerably the estimation performance over the case in which
the users operate individually (via OMP) and, in particular, for users facing
both the better and (especially) the worse conditions. Finally, it is observed
that WDiOMP is able to start from a lower $ASCE$ which is attributed at its
capability of weighting the contribution of each user, a behaviour that will
be even more pronounced in the following experiment.

Figure 2: $L_{g}=5$, $L_{c}=0$, $SNR=\\{10,20\\}$dB Figure 3: Different $SNR$
conditions per user. $L_{g}=5$, $L_{c}=0$.

In the third experiment (see Fig. 4), another “unbalanced” scenario is
examined. In more detail, the users are divided into two groups that share a
global support set (assuming $L_{g}=5$), while the users, of each group, share
also a common support set of size $L_{c}=3$. Additionally, all users face
$SNR=20$ dB. Let us recall that only the total sparsity level ($L=8$) is known
(namely, $L_{g}$ and $L_{c}$ are considered unknown). It is observed that
WDiOMP perfoms similarly to the previous experiments and considerably better
than SOMP and DiOMP. For the latter, the floor in their performance is a
result of the unknown $L_{g}$ and $L_{c}$. WiOMP does not have this issue
because it is able to infer (in a wide sense) the sparsity structure via the
employed weighting voting mechanism. Also, it is noted that SOMP, here, is
shown in order to assess the performance of DiOMP and should not be considered
as the centralized version of the examined scenario. Furthermore, WDiOMP is
better than OMP and requires about $40\%$ less measurements for achieving a
near-zero $ASCE$ (in the case of OMP no floor appears as only $L$ is
relevant). Furthermore, it is mentioned that, although not shown here, the
weights calculated by a particular user that employs WDiOMP, tend to similar
(and larger) values for users in the same group while the weights are smaller
for the remaining users, as hinted in Sec III-B. As a final remark, it is
noted that no $MSE$ curves for estimating $\boldsymbol{h}$ in (8) are
presented due to space limitations. However, the conclusions are similar to
the ones for $ASCE$. From the previous results, WDiOMP demonstrates a robust
performance for all considered experimental scenarios via the proposed
weighting voting mechanism.

Figure 4: Two groups of users face global and common support sets. $L_{g}=5$,
$L_{c}=3$), $SNR=20$dB.

##  V Conclusion

A fully distributed algorithm, WDiOMP, for the downlink channel estimation has
been presented that exploits a general sparsity support model for the involved
channels. The improved performance of WDiOMP has been assessed in various
scenarios, even when the structure of the sparsity support is unknown. Our aim
is to further investigate the entailed collaboration benefits in terms of
communication and energy consumption.

##  VI Acknowledgment

This research is co-financed by Greece and the European Union via
’Strengthening Human Resources Research Potential via Doctorate Research’
(MIS-5000432), implemented by the State Scholarships Foundation (IKY). This
work is also supported by the ERDF and the Republic of Cyprus through the
Project INFRASTRUCTURES/1216/0017 IRIDA.

## References

  * [1] A. Alkhateeb, J. Mo, N. Gonzalez-Prelcic, and R. W. Heath, “MIMO precoding and combining solutions for millimeter-wave systems,” _IEEE Communications Magazine_ , vol. 52, no. 12, pp. 122–131, 2014.
  * [2] J. A. Zhang, X. Huang, V. Dyadyuk, and Y. J. Guo, “Massive hybrid antenna array for millimeter-wave cellular communications,” _IEEE Wireless Communications_ , vol. 22, no. 1, pp. 79–87, 2015.
  * [3] R. W. Heath, N. Gonzalez-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave MIMO systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 10, no. 3, pp. 436–453, 2016.
  * [4] X. Rao and V. K. Lau, “Distributed compressive csit estimation and feedback for FDD multi-user massive MIMO systems,” _IEEE Trans. on Signal Processing_ , vol. 62, no. 12, pp. 3261–3271, 2014.
  * [5] A. Alkhateeb, O. El Ayach, G. Leus, and R. W. Heath, “Channel estimation and hybrid precoding for millimeter wave cellular systems,” _Journal of Selected Topics in Signal Processing_ , vol. 8, 2014.
  * [6] W. U. Bajwa, J. Haupt, A. M. Sayeed, and R. Nowak, “Compressed channel sensing: A new approach to estimating sparse multipath channels,” _Proceedings of the IEEE_ , vol. 98, no. 6, pp. 1058–1076, 2010.
  * [7] Z. Gao, L. Dai, Z. Wang, and S. Chen, “Spatially common sparsity based adaptive channel estimation and feedback for FDD massive MIMO,” _IEEE Trans. on Signal Processing_ , vol. 63, no. 23, pp. 6169–6183, 2015\.
  * [8] D. Baron, M. B. Wakin, M. F. Duarte, S. Sarvotham, and R. G. Baraniuk, “Distributed compressed sensing,” _Rice University, Depart. Electrical and Computer Engineering Technical Report TREE-0612_ , 2006.
  * [9] T. Wimalajeewa and P. K. Varshney, “OMP based joint sparsity pattern recovery under communication constraints,” _IEEE Trans. on Signal Processing_ , vol. 62, no. 19, pp. 5059–5072, 2014.
  * [10] D. Sundman, S. Chatterjee, and M. Skoglund, “Distributed greedy pursuit algorithms,” _Signal Processing_ , vol. 105, pp. 298–315, 2014.
  * [11] J. Plata-Chaves, N. Bogdanovic, and K. Berberidis, “Distributed diffusion-based LMS for node-specific adaptive parameter estimation,” _IEEE Trans. on Signal Processing_ , vol. 63, no. 13, pp. 3448–3460, 2015\.
  * [12] Y. Ding and B. D. Rao, “Dictionary learning-based sparse channel representation and estimation for FDD massive MIMO systems,” _IEEE Trans. on Wireless Communications_ , vol. 17, no. 8, pp. 5437–5451, 2018.
  * [13] F. Jameel, Z. Hamid, F. Jabeen, S. Zeadally, and M. A. Javed, “A survey of device-to-device communications: Research issues and challenges,” _IEEE Communications Surveys & Tutorials_, vol. 20, 2018.
  * [14] A. H. Sayed, S.-Y. Tu, J. Chen, X. Zhao, and Z. J. Towfic, “Diffusion strategies for adaptation and learning over networks: an examination of distributed strategies and network behavior,” _IEEE Signal Processing Magazine_ , vol. 30, no. 3, pp. 155–171, 2013.
  * [15] G. Li, T. Wimalajeewa, and P. K. Varshney, “Decentralized and collaborative subspace pursuit: A communication-efficient algorithm for joint sparsity pattern recovery with sensor networks,” _IEEE Trans. on Signal Processing_ , vol. 64, no. 3, pp. 556–566, 2015.
  * [16] R. Zhang, H. Zhao, and J. Zhang, “Distributed compressed sensing aided sparse channel estimation in FDD massive MIMO system,” _IEEE Access_ , vol. 6, pp. 18 383–18 397, 2018.