# Probing active nematics with in-situ microfabricated elastic inclusions

Ignasi Vélez-Cerón  Department of Materials Science and Physical Chemistry,
Universitat de Barcelona, 08028 Barcelona, Spain  Institute of Nanoscience and
Nanotechnology, IN2UB, Universitat de Barcelona, 08028 Barcelona, Spain  Pau
Guillamat  Institute for Bioengineering of Catalonia (IBEC), The Barcelona
Institute for Science and Technology (BIST), Barcelona, 08028, Spain  Francesc
Sagués  Department of Materials Science and Physical Chemistry, Universitat de
Barcelona, 08028 Barcelona, Spain  Institute of Nanoscience and
Nanotechnology, IN2UB, Universitat de Barcelona, 08028 Barcelona, Spain  Jordi
Ignés-Mullol  Department of Materials Science and Physical Chemistry,
Universitat de Barcelona, 08028 Barcelona, Spain  Institute of Nanoscience and
Nanotechnology, IN2UB, Universitat de Barcelona, 08028 Barcelona, Spain
<EMAIL_ADDRESS>
###### Abstract

In this work, we report a direct measurement of the forces exerted by a
tubulin/kinesin active nematic gel as well as its complete rheological
characterization, including the quantification of its shear viscosity, $\eta$,
and its activity parameter, $\alpha$. For this, we develop a novel method that
allows us to rapidly photo-polymerize compliant elastic inclusions in the
continuously remodelling active system. Moreover, we quantitatively settle
long-standing theoretical predictions, such as a postulated relationship
encoding the intrinsic time scale of the active nematic in terms of $\eta$ and
$\alpha$. In parallel, we infer a value for the nematic elasticity constant,
$K$, by combining our measurements with the theorized scaling of the active
length scale. On top of the microrheology capatilities, we demonstrate novel
strategies for defect encapsulation, quantification of defect mechanics, and
defect interactions, enabled by the versatility of the new microfabrication
strategy that allows to combine elastic motifs of different shape and
stiffness that are fabricated in-situ and on-time.

Active matter embraces a wide variety of systems composed of dense collections
of units that continuously convert energy into forces, thus being inherently
out of equilibrium [1]. Interactions between these force-generating units lead
to collective phenomena such as flocking, seemingly-turbulent flows, and other
emergent mechanodynamical patterns [2, 3]. Besides the many natural examples
that can be found at different length scales [4, 5, 6], similar phenomena can
be obtained in-vitro by bottom up approaches based on the densification of
active colloidal particles including, for instance, monolayers of mammalian
cells [7, 8, 9, 10], self-propelled particles [11, 12, 13, 14], living
bacteria [15, 16, 17, 18] or reconstitutions of cytoskeletal filaments and
motor proteins [19, 20].

The latter systems are among the most fascinating examples in the field of
synthetic active matter [21]. At the expenses of adenosine triphosphate (ATP)
consumption, bound molecular motor clusters internally shear protein filament
bundles, leading to continuous chaotic flows. In a dense, quasi-two-
dimensional phase, known as active nematic (AN) [22, 23], these active gels
feature continuously changing textures with liquid-crystal-like orientational
(nematic) order. The latter is lost in local regions that configure
topological defects.

Experiments have been able to characterize [24, 25, 26], even control [27, 28,
23] and manipulate the dynamics of AN realisations [29, 30, 31], and
theoretical and numerical works have successfully rationalized experimental
observations [32, 33, 34, 35, 36]. In spite of remarkable recent studies with
monolayers of mammalian cells featuring nematic order, [7] little is known
about the active mechanics of AN materials. More fundamentally, an accurate
quantitative determination of the most important constitutive parameters of
the microtubule-based system is lacking. Here, we address this challenge by
embedding microfabricated elastic inclusions with pre-designed size and
geometry within the AN layer. By studying the deformations of such objects we
have been able to probe the system’s rheology and mechanics.

Figure 1: In situ photo-polymerization of hydrogel structures in an active
nematic. (A) (Top) Sketch of the experimental cell and polymerization of four
rigid hydrogel rings within the AN (See SI Video S1 and also Fig. S1 for
further details of the setup). (Bottom) Fluorescence micrographs before
irradiation and during the emergence of the hydrogel structures, and the
subsequent adaptation of the AN to the sudden confinement. Elapsed times are
overlaid on each frame. The inset in the leftmost micrograph includes a sketch
with the microstructure of the AN filaments, where antiparallel microtubules
slide due to the action of kinesin dimers. The presence of a depletant agent
favours the condensation of the dispersed biomaterials. (B) Shorter
irradiation times lead to the formation of a deformable columnar ring
(fluorescence micrographs; see also SI Video S2). Elapsed times from the first
panel, where the ring has already been created, are overlaid on each frame.
The rightmost panel illustrates the deformation of the inclusion over time.
Scale bars, $50\;\mu$m.  Figure 2: Probing AN forces with elastic micro-
pillars. (A) Fluorescence micrograph of an AN layer at the instant when a
lattice of soft columns has been photo-polymerized (see also SI Movie S4).
Irradiation time is 1.5 s. Scale bar is $50\;\mu$m. (B) Color-coded location
of the columns’ free surface during the first 5 s after polymerization. (C)
The flexible columns are deflected by $\vec{\delta}$ due to the action of
$\vec{F}$, the resultant of the viscous drag and active forces acting at free
end of each column, which are submerged for a distance $\Delta z$ in the AN
layer. (D) Measured activity parameter and viscosity of the AN layer vs. the
average speed of the AN flows. In the inset, the square of the speed is
plotted vs. the activity parameter. (E) Comparison between the inverse of the
active time scale, $\alpha/\eta$, computed from the parameters reported in
panel (D), and the average vorticity, $\omega$, obtained from the active flow
velocity map. In (D) and (E), error bars are the standard deviation of the
mean for $v$ and $\omega$, and confidence intervals for the fitted parameters
$\alpha$ and $\eta$.

To this end, we tune our active preparations by adding to the standard
formulation of a kinesin/tubulin active gel the precursors of a hydrogel that
can be photo-polymerized with UV light. This new mixture is introduced in a
flow cell with a gap of 50 microns between two glass plates (see SI, Sec. I,
for full details). One plate incorporates a superhydrophilic polyacrylamide
(PAM) brush, which prevents protein adsorption. The other plate is
functionalized with fluorinated alkane chains that adsorb a thin layer of a
fluorinated oil. The active nematic forms at the water/oil interface, which is
decorated with a polyethyleneglycol-based surfactant that prevents
biomolecules from being in direct contact with the oil. As stated above, apart
from the standard components of the active preparation, the aqueous subphase
incorporates polyethyleneglycol-acrylate monomers and lithium
phenyl-2,4,6-trimethylbenzoylphosphinate (LAP)[37], a photosensitive molecule
(see Table S1). Irradiation with UV light induces the dissociation of LAP into
radical species that initiate a chain polymerization reaction leading to the
formation of a PEG hydrogel in the regions exposed to the UV light. In our
experiments, we use a custom illumination setup based on a digital micromirror
device (DMD) projector integrated into the light path of an inverted
fluorescence microscope (Fig. S1). Once the AN is formed, patterns of UV light
are focused on the AN layer by means of the microscope objective, which sets
the maximum size and lateral resolution of the built objects, as well as the
range of power densities of the irradiation.

With this setup, we can generate in-situ within the active gel columnar
hydrogel microstructures of arbitrary shape, set by the illumination pattern,
and stiffness, according to the irradiation time and footprint of the objects.
Columns are grafted on the PAM-coated plate and extend across the aqueous
phase, up to the AN layer in contact with the oil phase. In Fig. 1A we show
schematically the process to create rigid hydrogel columns of arbitrary shape
anywhere and at any time on an evolving AN layer. As an example, in Fig. 1A
four adjacent rings are created, as seen in the fluorescence micrographs below
the sketches. The AN layer adapts to the new boundary conditions by
transforming the local defect pattern to accommodate the imposed topology.
Notice that the fluorescence probe contained in the hydrogel volume is
photobleached by the action of UV light, thus appearing dark in the
fluorescence micrographs upon polymerization. In Fig. 1B, we demonstrate the
system’s capability of generating also soft columnar objects. In this case, a
flexible ring that is continuously deformed by forces, generated by the
interplay of the ANs inside and outside the inclusion. Both rigid and flexible
walls trigger the formation of active boundary layers [38] that determine the
geometry of the enclosed AN. The high versatility of this protocol allows to
print, either in sequence or simultaneously, both rigid and soft, deformable
structures. This can be demonstrated by imprinting a combination of flexible
pillars enclosed by a rigid ring; see SI. Video S3. This allows both to
monitor the adaptation of the AN to a sudden change in boundary conditions,
and to probe the local stresses exerted by the active material, as we discuss
below.

In this study, we take advantage of the elasticity of the soft hydrogel to
imprint an array of compliant pillars with circular cross-section as
microrheological probes for the AN layer (Fig. 2A,B). Each pillar will act as
a cantilever, anchored on the hydrophilic plate and bent by the action of the
AN layer. The active material exerts a combination of hydrodynamic, active,
and elastic stresses that results in the so-called active turbulent regime
[32, 26]. Under mechanical equilibrium, column bending will be the result of
the net forces arising in the active fluid. Previous studies with this
material suggested that hydrodynamic and active stresses have similar orders
of magnitude, while elastic stresses, being of higher order in the gradients
of the orientational field, are subdominant except very near topological
defects [39, 40]. We will take this into account and assume in our system that
column deflection is balanced by the effect of active and hydrodynamic (shear)
stresses alone.

Figure 3: Active nematic forces. (A) Force density map, obtained from
$f_{i}=2\eta\partial_{j}\varepsilon_{ij}-\alpha\partial_{j}Q_{ij}$, in the
vicinity of a $+1/2$ defect. The force density accumulates at the tip of the
self-propelled defect. (B) In the vicinity of a $-1/2$ defect, the force
density exhibits nearly three-fold symmetry. (C) Force density map during the
unbinding of a defect pair. The force density accumulates at the tip of the
bending filaments, leading to the formation of a $+1/2$ and $-1/2$ defect
pair. pair. Elapsed times from the first frame are 2 s, 8 s, and 12 s (see
also SI Video S5). (D) Force density map during the annihilation of a defect
pair. The high force density at the tip of the moving $+1/2$ defect vanishes
after coalescence. Elapsed times from the first frame are 9 s, 14 s, and 19 s
(see also SI Video S6). Scale bars, 50 $\mu$m.

For the characterization of active forces, we balance the elastic pillar
deflection force (see SI Sec. II.1) with the stresses within the active
material, $\sigma^{AN}_{ij}$, computed on the thin corona of the hydrogel
cylinder of thickness $\Delta z$ that is submerged in the AN layer (Fig. 2C).
In the following, we also assume that the force is applied at the edge of the
free surface of the pillar. Moreover, we ignore the viscous drag due to the
oil layer in contact with the upper cylinder’ surface and that due to the
aqueous volume that embraces most of the column, both negligible (see SI Sec.
II.2 for details of this analysis).

The force exerted by the AN on a pillar is computed as
$F_{i}=\oint\sum_{j=1,2}\sigma^{AN,2d}_{ij}n_{j}dl$, where $i,j$ span the two
spatial dimensions, $n_{i}$ is the unit vector normal to the lateral surface
of the column, the line integral spans the circular cross section of the
column, and $\sigma^{AN,2d}_{ij}$ is the two-dimensional AN stress tensor (see
SI Sec. II.1). In terms of the dominant viscous and active stresses,
$\sigma^{AN,2d}_{ij}=2\eta\varepsilon_{ij}-\alpha Q_{ij}$, where
$\varepsilon_{ij}$ and $Q_{ij}$ denote, respectively, the 2d strain rate
tensor and nematic tensor order parameter. The parameters $\eta$ and $\alpha$,
which we measure here, stand respectively for the 2d shear viscosity and
activity coefficient of the AN (assumed positive for our extensile system
[22]).

On the other hand, the force required for the bending deflection of the pillar
is $F_{i}^{def}=\frac{3\delta_{i}EI}{L^{3}}$, where $\delta_{i}$ denotes the
in-plane component of the deflection, while $E$, $I$, and $L$ correspond
respectively to the Young’s modulus (which we determine using fluids of known
viscosity, see SI, Sec.II.3), the second moment of inertia of the cylindrical
column, and the height of the pillars, set by the spacing between glass plates
(see Fig. 2C).

After their photopolymerization within the AN, we will assume that each pillar
achieves mechanical equilibrium either when its deflection reaches a steady-
state value or just before it changes direction (see SI Sec. II.1). Moreover,
local shear and active stress will be computed from the instantaneous
orientation and velocity fields of the AN, as defined above, just before
pillar polymerization. In our analysis, we have included columns of different
diameters in the range $10.4\;\mu$m-$18\;\mu$m, to achieve a deflection at
mechanical equilibrium in the range $2\;-14\;\mu$m for all ATP concentrations.
Our measurements show that the total force exerted on a single column is in
the range 1 - 20 pN, depending on the column diameter and local fields of the
AN. In all cases, the drag force due to active flows is typically 3-4 times
larger than the corresponding force due to active stress.

For each pillar, we write an equation for each of the two spatial components
in the force balance condition, $F_{i}=F^{def}_{i}$, which depends on $\alpha$
and $\eta$ (see SI Sec. II.1). We find the optimum values for $\alpha$ and
$\eta$ that more closely satisfy the force balance equations obtained for a
set including dozens of pillars for different experimental realizations
performed under the same experimental conditions (see SI Sec. II.1). The
results of our measurements are shown in Fig. 2B, where $\alpha$ and $\eta$
are plotted as a function of the average speed of the active flow. The latter
is used here as a reliable proxy for the system’s activity. A calibration
curve of speed vs. ATP concentration is presented in Figure S3A, which allows
to plot $\alpha$ and $\eta$ as a function of ATP concentration (Fig. S3B). Our
results show that $\alpha$ is around $2\;10^{-7}$Pa m, only weakly dependent
on activity, while $\eta$ is in the range $(4-10)\;10^{-6}$ Pa s m, revealing
a shear thinning nature of the active material, as previously reported [41].

We have compared our direct measurement of $\eta$ with an early strategy
reported by Martínez-Prat et al. [26] based on a model for the kinetic energy
spectra, and have found excellent agreement (see SI Sec. II.4 and Figure S4).
Moreover, earlier theoretical studies predicted a scaling relationship between
the activity parameter and the average speed of the form $v^{2}\sim\alpha$
[32], which we have validated with our experimental data (see inset in Fig.
2D). We have also verified that our measurements agree with the theoretical
expectation that $\alpha/\eta$ is related to the inverse of the intrinsic time
scale of the AN, which we can extract from the average vorticity, $\omega$.
Indeed, Fig. 2D suggests that these two parameters are proportional,
$\omega\sim\alpha/\eta$ [32, 39].

Having obtained reliable values for $\alpha$ and $\eta$ at different
activities, we can estimate a third material parameter, $K$, which encodes the
bending rigidity of the active filaments. For this purpose, we resort to the
scaling relation for the active length scale, $l_{\alpha}^{2}\sim K/\alpha$
[22]. In our experiments, we estimate $l_{\alpha}\sim n^{-1/2}$, where $n$ is
the average defect density, from which $K\sim\alpha/n$ can be estimated for
each concentration of ATP (Fig. S5). Our results suggest an average value
$K\simeq 6\;10^{-16}$N m, with a marginal increase with the concentration of
ATP, indicating that microtubule bundling may be enhanced by ATP.

Figure 4: Mixed rigid and compliant objects for probing forces in confined AN
arrangements. (A) Two wide rings impose rigid confinement on the AN layer,
leading to the reorganization of the active flow into a pair of rotating
$+1/2$ defects. Once the flow is steady, we imprint additional compliant
columns, either at the center (top panel) or off-centered (botom panel) with
respect to the rigid ring (see also SI Video S7). (B) Compliant columns
deformation during the first 14 seconds after printing. (C) A rigid ring and a
set of flexible columns are simultaneously imprinted around a $-1/2$ defect
using 1.5 s of irradiation (see also SI Video S8). (D) Column deformation
during the first 8 s after printing. The scale bar is 50 $\mu$m.

Having quantified the material parameters that determine the main local stress
of the AN, the presented method enables us to quantitatively map, for any AN
configuration, the resultant force density field, obtained as $f_{i}\simeq
2\eta\partial_{j}\varepsilon_{ij}-\alpha\partial_{j}Q_{ij}$. Note that here we
have also neglected the contribution from the AN elasticity, as discussed
above. In Fig. 3 we apply this result to map the total force density field
around a $+1/2$ (Fig. 3A) and a $-1/2$ (Fig. 3B) defect. The force patterns
for each configuration, respectively parabolic and triangular, are consistent
with previously reported theoretical results [22, 42]. Experimental force
patterns show that $+1/2$ defects self-propel from tail to head dragging the
surrounding material. Conversely, we see that negative defects balance the
local forces acting in their surroundings, leading to no net propulsion and
exhibiting advected motion only.

It is interesting to compare the relative strengths of the two dominant force
contributions in the AN layer in the absence of obstacles. To this purpose, we
have resolved the force density field into its active and viscous components
near a $+1/2$ defect (SI Fig. S6). The balance in this case is dominated by
active forces, contrary to what was observed when computing the force balance
for column deflection. In that case, viscous drag arising from no-slip
boundary conditions on the hydrogel cylindrical surface resulted in strong
forces in the direction of the average flow. Without obstacles, viscous force
density is relatively weak and it opposes the self-propelling flows generated
by the local active stress.

The propelling effect from accumulated internal stress within aligned
microtubule bundles is best illustrated during the unbinding of defect pairs
of opposite topological charge, one of the most distinctive events in the
dynamics of the quasi-two-dimensional ANs. One of this events is illustrated
in Fig. 3C, where the total force density increases at the tip of an unstable
bending region. This promotes its growth, until a $+1/2$ defect unbinds, being
self-propelled by a strong force density at its tip, while a $-1/2$ defect is
left behind, surrounded by a stagnating force density field. The complementary
scenario, namely that of defect annihilation, is shown in Fig. 3D. The $+1/2$
defect is propelled by a strong force field acting on its tip towards a nearly
stationary $-1/2$ defect. Recombination leads to a locally defect-free region
where the force field vanishes.

The versatility of the employed method allows us to polymerize a combination
of rigid and compliant inclusions. This opens the door to measuring the force
fields generated in confined AN arrangements. Some examples are presented in
Fig. 4. Fig. 4A and B includes a system of two rotating positive defects
created after active flows adapt to sudden ring confinement. Stresses are
probed with one compliant column centered in one case, and off-center in the
other. Column bending in each situation clearly reflects the stress pattern
created by the AN, with small deformations near the center (top) and large
deformations in the off-centered case (bottom). Fig. 4C and D explore, in a
similar fashion, a $-1/2$ defect captured within a rigid ring. Here, the AN is
prone to destabilization under confinement. An array of five compliant pillars
seeded in a centered-square arrangement exhibit an inappreciable deformation
of the central probe and small, albeit measurable, deformations of the outer
columns that indicate a torque pattern.

The reported strategy of hybridizing the microtubule/kinesin system with a
photopolymerizable monomer not only deepens significatively our understanding
of the constitutive characteristics of this active gel but opens side by side
endless possibilities to control the active nature of this particular and
other similar materials. The incorporation of active materials with
externally-tunable characteristics into microchips represents a significant
step forward in the design of future micro-machines built on the principles of
active matter.

Data availability. All data is available form the corresponding author upon
reasonable request.

###### Acknowledgements.

The authors are indebted to the Brandeis University MRSEC Biosynthesis
facility for providing the tubulin. We thank M. Pons, A. LeRoux, and G. Iruela
(Universitat de Barcelona) for their assistance in the expression of motor
proteins. I.V. acknowledges funding from Generalitat de Catalunya though a
FI-2020 PhD. Fellowship. P.G. acknowledges support from Generalitat de
Catalunya through the Beatriu de Pin´os program (grant number 2020 BP 00248).
I.V., J.I.-M., and F.S. acknowledge funding from
MICINN/AEI/10.13039/501100011033 (Grant No. PID2019-108842GB-C22). Brandeis
University MRSEC Biosynthesis facility is supported by NSF MRSEC 2011846. The
authors acknowledge helpul discussions with C. S. and N. I-M.

## References

  * Sagués [2023] F. Sagués,  _Colloidal Active Matter: Concepts, Experimental Realizations, and Models_ (Taylor & Francis, London, 1st. ed., 2023).
  * Ramaswamy [2010] S. Ramaswamy, The Mechanics and Statistics of Active Matter, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Reviews of Modern Physics 85, 1143 (2013), 1207.2929 .
  * Herbert-Read  _et al._ [2011] J. E. Herbert-Read, A. Perna, R. P. Mann, T. M. Schaerf, D. J. T. Sumpter, and A. J. W. Ward, Inferring the rules of interaction of shoaling fish, Proceedings of the National Academy of Sciences 108, 18726 (2011).
  * Bialek  _et al._ [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Statistical mechanics for natural flocks of birds, Proceedings of the National Academy of Sciences 109, 4786 (2012).
  * Buhl  _et al._ [2006] J. Buhl, D. J. T. Sumpter, I. D. Couzin, J. J. Hale, E. Despland, E. R. Miller, and S. J. Simpson, From disorder to order in marching locusts, Science 312, 1402 (2006).
  * Saw  _et al._ [2017] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Topological defects in epithelia govern cell death and extrusion, Nature 544, 212 (2017).
  * Kawaguchi  _et al._ [2017] K. Kawaguchi, R. Kageyama, and M. Sano, Topological defects control collective dynamics in neural progenitor cell cultures, Nature 545, 327 (2017).
  * Blanch-Mercader  _et al._ [2018] C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi, and P. Silberzan, Turbulent dynamics of epithelial cell cultures, Physical Review Letters 120, 208101 (2018), arXiv:1711.01568 .
  * Duclos  _et al._ [2018] G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J. F. Joanny, J. Prost, and P. Silberzan, Spontaneous shear flow in confined cellular nematics, Nature Physics 14, 728 (2018).
  * Theurkauff  _et al._ [2012] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, Dynamic clustering in active colloidal suspensions with chemical signaling, Physical Review Letters 108, 268303 (2012).
  * Palacci  _et al._ [2013] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and P. M. Chaikin, Living crystals of light-activated colloidal surfers, Science 339, 936 (2013).
  * Buttinoni  _et al._ [2013] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles, Physical Review Letters 110, 238301 (2013).
  * Wang  _et al._ [2015] W. Wang, W. Duan, S. Ahmed, A. Sen, and T. E. Mallouk, From One to Many: Dynamic Assembly and Collective Behavior of Self-Propelled Colloidal Motors, Accounts of Chemical Research 48, 1938 (2015).
  * Zhang  _et al._ [2010] H. P. Zhang, A. Be’er, E.-L. Florin, and H. L. Swinney, Collective motion and density fluctuations in bacterial colonies, Proceedings of the National Academy of Sciences 107, 13626 (2010).
  * Wensink  _et al._ [2012] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Meso-scale turbulence in living fluids, Proceedings of the National Academy of Sciences 109, 14308 (2012).
  * Li  _et al._ [2019] H. Li, X. Shi, M. Huang, X. Chen, M. Xiao, C. Liu, H. Chaté, and H. P. Zhang, Data-driven quantitative modeling of bacterial active nematics, Proceedings of the National Academy of Sciences 116, 777 (2019).
  * Turiv  _et al._ [2020] T. Turiv, R. Koizumi, K. Thijssen, M. M. Genkin, H. Yu, C. Peng, Q. H. Wei, J. M. Yeomans, I. S. Aranson, A. Doostmohammadi, and O. D. Lavrentovich, Polar jets of swimming bacteria condensed by a patterned liquid crystal, Nature Physics 16, 481 (2020).
  * Sanchez  _et al._ [2012] T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Spontaneous motion in hierarchically assembled active matter, Nature 491, 431 (2012).
  * Martínez-Prat  _et al._ [2019] B. Martínez-Prat, J. Ignés-Mullol, J. Casademunt, and F. Sagués, Selection mechanism at the onset of active turbulence, Nature Physics 15, 362 (2019).
  * Needleman and Dogic [2017] D. Needleman and Z. Dogic, Active matter at the interface between materials science and cell biology, Nat. Rev. Mat. 2, 17048 (2017).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Active nematics, Nat. Commun. 9, 3246 (2018).
  * Zhang  _et al._ [2021] R. Zhang, S. A. Redford, P. V. Ruijgrok, N. Kumar, A. Mozaffari, S. Zemsky, A. R. Dinner, V. Vitelli, Z. Bryant, M. L. Gardel, and J. J. de Pablo, Spatiotemporal control of liquid crystal structure and dynamics through activity patterning, Nature Materials 20, 875 (2021).
  * Guillamat  _et al._ [2016a] P. Guillamat, J. Ignés-Mullol, S. Shankar, M. C. Marchetti, and F. Sagués, Probing the shear viscosity of an active nematic film, Phys. Rev. E 94, 060602 (2016a).
  * Tan  _et al._ [2019] A. J. Tan, E. Roberts, S. A. Smith, U. A. Olvera, J. Arteaga, S. Fortini, K. A. Mitchell, and L. S. Hirst, Topological chaos in active nematics, Nature Physics (2019).
  * Martínez-Prat  _et al._ [2021] B. Martínez-Prat, R. Alert, F. Meng, J. Ignés-Mullol, J.-F. Joanny, J. Casademunt, R. Golestanian, and F. Sagués, Scaling regimes of active turbulence with external dissipation, Physical Review X 11, 031065(1 (2021).
  * Guillamat  _et al._ [2016b] P. Guillamat, J. Ignés-Mullol, and F. Sagués, Control of active liquid crystals with a magnetic field, Proceedings of the National Academy of Sciences 113, 5498 (2016b).
  * Guillamat  _et al._ [2017] P. Guillamat, J. Ignés-Mullol, and F. Sagués, Taming active turbulence with patterned soft interfaces, Nat. Commun. 8 (2017).
  * Keber  _et al._ [2014] F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch, Topology and dynamics of active nematic vesicles, Science 345, 1135 (2014).
  * Wu  _et al._ [2017] K. T. Wu, J. B. Hishamunda, D. T. N. Chen, S. J. DeCamp, Y. W. Chang, A. Fernández-Nieves, S. Fraden, and Z. Dogic, Transition from turbulent to coherent flows in confined three-dimensional active fluids, Science 355, 10.1126/science.aal1979 (2017).
  * Hardoüin  _et al._ [2019] J. Hardoüin, R. Hughes, A. Doostmohammadi, J. Laurent, T. Lopez-Leon, J. M. Yeomans, J. Ignés-Mullol, and F. Sagués, Reconfigurable Flows and Defect Landscape of Confined Active Nematics, Commun. Phys. 2, 121 (2019).
  * Giomi [2015] L. Giomi, Geometry and topology of turbulence in active nematics, Phys. Rev. X 5, 031003 (2015).
  * Zhang  _et al._ [2016] R. Zhang, Y. Zhou, M. Rahimi, and J. J. De Pablo, Dynamic structure of active nematic shells, Nature Communications 7, 13483 (2016).
  * Alert  _et al._ [2020] R. Alert, J.-F. Joanny, and J. Casademunt, Universal scaling of active nematic turbulence, Nature Physics 16, 682 (2020), arXiv:1906.04757 .
  * Pearce  _et al._ [2019] D. J. G. Pearce, P. W. Ellis, A. Fernandez-Nieves, and L. Giomi, Geometrical Control of Active Turbulence in Curved Topographies, Physical Review Letters 122, 168002 (2019).
  * Thijssen  _et al._ [2021] K. Thijssen, D. A. Khaladj, S. A. Aghvami, M. A. Gharbi, S. Fraden, J. M. Yeomans, L. S. Hirst, and T. N. Shendruk, Submersed micropatterned structures control active nematic flow, topology, and concentration, Proceedings of the National Academy of Sciences of the United States of America 118, 1 (2021), arXiv:2102.10184 .
  * Fairbanks  _et al._ [2009] B. D. Fairbanks, M. P. Schwartz, C. N. Bowman, and K. S. Anseth, Photoinitiated polymerization of peg-diacrylate with lithium phenyl-2,4,6-trimethylbenzoylphosphinate: polymerization rate and cytocompatibility, Biomaterials 30, 6702 (2009).
  * Hardoüin  _et al._ [2022] J. Hardoüin, C. Doré, J. Laurent, T. Lopez-Leon, J. Ignés-Mullol, and F. Sagués, Active boundary layers in confined active nematics, Nature Communications 13, 6675 (2022).
  * Joshi  _et al._ [2022] C. Joshi, S. Ray, L. M. Lemma, M. Varghese, G. Sharp, Z. Dogic, A. Baskaran, and M. F. Hagan, Data-driven discovery of active nematic hydrodynamics, Phys Rev Lett 129, 258001 (2022).
  * Golden  _et al._ [2023] M. Golden, R. O. Grigoriev, J. Nambisan, and A. Fernandez-Nieves, Physically informed data-driven modeling of active nematics, Sci Adv 9, eabq6120 (2023).
  * Rivas  _et al._ [2020] D. P. Rivas, T. N. Shendruk, R. R. Henry, D. H. Reich, and R. L. Leheny, Driven topological transitions in active nematic films, Soft Matter 16, 9331 (2020).
  * Serra  _et al._ [2023] M. Serra, L. Lemma, L. Giomi, Z. Dogic, and L. Mahadevan, Defect-mediated dynamics of coherent structures in active nematics, Nature Physics (2023).
  * Ellis  _et al._ [2020] P. W. Ellis, J. Nambisan, and A. Fernandez-Nieves, Coherence-enhanced diffusion filtering applied to partially-ordered fluids, Molecular Physics , 1 (2020).
  * Hibbeler [2010] R. C. Hibbeler,  _Mechanics of Materials_ (Prentice Hall, New Jersey, 8th ed., 2010).
  * Landau and M. [1987] L. D. Landau and L. E. M.,  _Fluid Mechanics_ (Pergamon Press, Oxford, 2nd ed., 1987).

SUPPLEMENTARY MATERIAL

Probing active nematics with in-situ microfabricated elastic inclusions
Vélez-Cerón et al.

##  I Methods

###  I.1 Polymerization of microtubules

Microtubules ($\sim 1\;\mu$m) were polymerized from heterodimeric
($\alpha$,$\beta$)-tubulin from bovine brain (Brandeis University Biological
Materials Facility). Tubulin was incubated at 37 ∘C for 30 min in aqueous M2B
buffer (80 mM Pipes (piperazine-N,N’-bis(2-ethanesulfonic acid)), 1 mM EGTA
(ethylene glycol-bis($\beta$-aminoethyl ether)-N,N,N’,N’-tetraacetic acid),
2mM MgCl2) (Sigma; P1851, E3889 and M4880, respectively) prepared with Milli-Q
water and supplemented with the antioxidant agent dithiothrethiol (DTT) (Sigma
43815) and non-hydrolysable guanosine triphosphate (GTP) analog GMPCPP
(guanosine-5’-[($\alpha$,$\beta$)-methyleno]triphosphate) (Jena Biosciences,
NU-405) up to a concentration of 1 mM and 0.6 mM, respectively. To allow the
observation under fluorescence microscopy, 3% of the total tubulin
concentration (8 mg ml-1) was labelled with Alexa 647. Afterward, the mixture
was annealed at room temperature for 4 h and kept at -80 ∘C until use.

###  I.2 Kinesin expression

The kinesin used in the experiments was the heavy-chain kinesin-1 from
Drosophila melanogaster truncated at residue 401, fused to biotin carboxyl
carrier protein (BCCP), and labelled with six histidine tags (K401-BCCP-6His).
This was expressed in Escherichia coli using the plasmid WC2 from the Gelles
Laboratory (Brandeis University) and purified with a nickel column. After
dialysis against 500 mM imidazole aqueous buffer, kinesin concentration was
estimated by means of absorption spectroscopy. Finally, the kinesin was stored
in a 40% (wt/vol) aqueous sucrose solution at -80 ∘C until use.

###  I.3 Assembly of the microtubule-based active mixture

Biotinylated kinesin motor proteins were incubated on ice for 30 min with
tetrameric streptavidin (Invitrogen; 43-4301) at the specific ratio of 2:1 to
obtain kinesin-streptavidin motor clusters. The standard preparation of the
active mixture consisted of a M2B aqueous solution that contains ATP
(adenosine 5’-triphosphate) (Sigma; A2383), the motor clusters, microtubules,
and the depleting agent PEG (poly(ethylene glycol)) (20 kDa, Sigma;95172). To
maintain a constant concentration of ATP during the experiments, an enzymatic
ATP-regeneration system was added, consisting on phosphoenolpyruvate (PEP)
(Sigma; P7127) that fueled the enzyme pyruvate kinase/lactate dehydrogenase
(PKLDH) (Sigma; P0294) to convert ADP (adenosine 5’-diphosphate) back to ATP.
The mixture also contained several antioxidants to avoid photobleaching during
the fluorescence imaging: DTT, Glucose Oxidase (Sigma; GT141), catalase
(Sigma; C40), Trolox (6-hydroxy-2,5,7,8-tetramethylchroman-2-carboxylic acid)
(Sigma; 238813), and glucose (Sigma; G7021). The constituents of the
photopolymerizing hydrogel were incorporated in solid form to the active
mixture. The former consisting on the photoinitiator lithium
phenyl-2,4,6-trimethylbenzoylphosphinate (LAP, TCI; L0290) and the monomer
4-ArmPEG-Acrylate 5 kDa (4PEG5k, Biochempeg; A44009-5k). The final
concentrations of each reagent used in the preparation of the active material
is shown in Table S1.

Compound | Buffer | Final Conc. | Units
---|---|---|---
PEG (20 kDa) | M2B | 1.54 | % w/v
PEP | M2B | 25.68 | mM
MgCl2 | M2B | 3.12 | mM
ATP | M2B | 60 - 296 |  $\mu$M
DTT | M2B | 5.21 | mM
Streptavidin | M2B | 0.01 | mg/mL
Trolox | Phosphate | 1.93 | mM
Catalase | Phosphate | 0.04 | mg/mL
Glucose | Phosphate | 3.20 | mg/mL
Glucose Oxidase | Phosphate | 0.21 | mg/mL
PK | Original | 25.01 | u/mL
LDH | Original | 24.91 | u/mL
Kinesin | Original | 0.08 | mg/mL
Microtubules | Original | 1.85 | mg/mL
LAP | Mixture | 0.226 | % w/v
4PEG5k | Mixture | 4.52 | % w/v
Table S1: Composition of all solutions (including buffer used for their
preparation), and concentration of the different species in the final mixture.
Acronyms used in this table are: PEG (Poly-ethylene glycol); PEP (Phosphoenol
pyruvate); ATP (Adenosin triphosphate); PK (Pyruvate Kinase); LDH (Lactic
Dehydrogenase); DTT (1,4-dithiothreitol); 4PEG5k (4-ArmPEG-Acrylate 5 kDa);
LAP (lithium phenyl-2,4,6-trimethylbenzoylphosphinate). The two latter
compounds where directly incorporated to the active mixture. M2B buffer: 80 mM
PIPES (piperazine-N,N’-bis(2-ethanesulfonic acid)) pH 6.8, 2 mM MgCl2 1 mM
EGTA (egtazic acid). Phosphate buffer: 20 mM Phosphate buffer (6.68 mM KH2PO4,
12.32 mM K2HPO4) pH 7.2; Original: species is obtained already dissolved in
its custom buffer.

###  I.4 Active nematic cell

The experiments were carried out using flow cells with a channel of 3-5 mm
width and 50 $\mu$m height. The cell was composed of a bioinert and
superhydrophilic polyacrylamide-coated glass and a hydrophobic Aquapel-coated
glass, separated using a 50 $\mu$m thick double-sided tape. The cell was first
filled by capillarity with fluorinated oil (HFE7500; Fluorochem 051243) which
contains 2% of a fluorosurfactant copolymer (RanBiotechnologies, 008
Fluorosurfactant). The active material was subsequently introduced in the cell
by capillarity. To avoid evaporation, the cell was sealed using petroleum
jelly. This preparation resulted in the approximate thickness for each phase
(as determined using fluorescence confocal microscopy): 10 $\mu$m for the oil
layer, 5-10 $\mu$m for the AN layer, and 30-35 $\mu$m for the aqueous
subphase.

###  I.5 Observation of the active nematic

The observation of the active nematic layer was carried out by means of
fluorescence microscopy. We used a Nikon Eclipse Ti2-e equipped with a white
LED source (Thorlabs MWWHLP1) and a Cy5 filter cube (Edmund Optics). Images
were captured using an ExiBlue CCD camera (QImaging) operated with the open-
source software ImageJ $\mu$-Manager.

###  I.6 Hydrogel Polymerization

The fluorescence microscope was modified in order to incorporate a digital
micromirror device (DMD) UV light projector (Texas Instruments LightCrafter
4500, with a 2W 385nm LED; EKB Technologies, Ltd), see Fig. S1. Projected
patterns are incorporated into the light path of the inverted microscope by
means of a collimating lens (f = +150 mm) and a 505 nm dichroic mirror
(Thorlabs DMLP505R), and are focused on the sample by means of the microscope
objectives, reaching a lateral resolution up to a few microns. The DMD
projector is connected as an external monitor to a computer, thus enabling
real time control of the projected patterns using MS-PowerPoint slides. For
the experiments shown in this work, the hydrogel was polymerized using an x20
objective, resulting in a light power density of 3.1 W cm-2. The resulting
columns have a height of $40\;\mu$m, of which 5-10 $\mu$m are embedded in the
AN layer.

Figure S1: Optical setup as described in the text.

The rigidity of the hydrogel structures can be adjusted by a combination of
the irradiation time and the power of the led light source that feeds the DMD
device. Since the AN features self-sustained flows that drag the water
subphase, polymerization must be performed quickly to prevent the formation of
distorted shapes, as the incipient hydrogel is advected by the active flows.
In practice, we require that the maximum flow advection during polymerization
is under 2 $\mu$m s-1, which limits ATP concentration to be below 450 $\mu$M.
With the used concentration of monomer and initiator, we find that the
polymerization time must be in the range 750 - 1250 ms for soft and 1500 -
3000 ms for rigid structures. As an alternative to adjusting the
polymerization time, the lateral compliance of the structures can be set by
controlling their in-plane cross section.

###  I.7 Sample characterization

The director field of the active nematic layer was obtained using a coherence-
enhanced diffusion filtering (CEDF) MatLab code [43] applied on fluorescence
images. Tracer-free velocimetry analysis of the active nematic was performed
with the public MatLab App PIVlab. Further analysis of director and velocity
field data was performed with custom-written MatLab scripts.

###  I.8 Young’s modulus determination

The hydrogel Young’s modulus was determined in-situ, in a flow cell prepared
as explained above, filled by capillarity with the buffer solution and
gellifying agents but without motor proteins and microtubules. Cylindrical
columns were polymerized using the protocol described above and were deformed
by flowing liquids of known viscosity by capillarity. We used methylcellulose
(Sigma; M0512) solutions of concentrations between 0.18 and 0.54% (w/vol)
prepared in M2B buffer. Their viscosity was measured using a Brookfield Ametek
DVE Viscosimeter with a Yula-15 spindle. To measure the average flow speed, a
5% (vol/vol) concentration of fluorescent particles of diameter 1 $\mu$m
(Bangs Laboratories; FCDG006) were incorporated into the solution, and average
speeds where determined by tracking individual particles using ImageJ.

##  II Data Analysis

###  II.1 Force balance

The viscosity, $\eta$, and the activity parameter, $\alpha$, of the active
nematic were determined by analysing the in-plane ($xy$) resolved deflection,
$\delta_{i}$, of the hydrogel cantilevers. Invoking basic concepts in
mechanics of materials [44], the differential equation for the deflection can
be formulated in terms of the bending moment, $M$, of a cylindrical pillar,
its second moment of inertia, $I=\pi r^{4}/4$, where $r$ denotes the radius of
the beam, and the Young’s modulus of the material, $E$,

$EI\frac{d^{2}\delta_{i}}{dz^{2}}=M_{i}(z),$ (1)

where the index $i$ refers to each of the two spatial dimensions. The
cantilevers, of length $L$, were fixed on the substrate. Their bending moment
is equal to the deflection force multiplied by the distance from the position
where it is applied, $M_{i}(z)={F_{i}}^{def}z$.

By integrating Eq. 1 with the corresponding boundary conditions
$\delta_{i}(z=L)=\frac{d\delta_{i}}{dz}|_{z=L}=0$, the final relationship
between the deflection and the applied force components reads,

$F_{i}^{def}=\frac{3\delta_{i}EI}{L^{3}}=\frac{3\pi}{4}\frac{\delta_{i}Er^{4}}{L^{3}}$
(2)

expressed in terms of the geometrical parameters of a beam and its Young’s
modulus .

This deflection force is balanced with the resultant from adding the force
directly exerted by the active fluid, and those indirectly arising in the
driven contact fluids, which are hydrodynamically coupled to the AN layer. The
force acting on the submerged columns can be expressed in terms of the total
stress acting on them, $\sigma_{ij}$, as

$F_{i}=\iiint_{V_{c}}\sum_{j=1,2}\partial_{j}\sigma_{ij}dV=\iint_{S_{c}}\sum_{j=1,2}\sigma_{ij}n_{j}dS.$
(3)

Here, $V_{c}$ spans the volume of the cylinder, $S_{c}$ spans the exposed area
of the cylinder (lateral and free circular surface), and $\hat{n}$ is the unit
vector normal to the local surface. We split this computation into three
parts,

$\begin{split}F_{i}&=F_{i}^{oil}+F_{i}^{AN}+F_{i}^{water}\\\
&=\iint_{S_{o}}\sigma^{o}_{ij}n_{j}dS+\iint_{S_{AN}}\sigma^{AN}_{ij}n_{j}dS+\iint_{S_{w}}\sigma^{w}_{ij}n_{j}dS,\end{split}$
(4)

where the indices $o$, $AN$, and $w$ stand, respectively, for the oil layer in
contact with the free top cylinder surface, the cylindrical region of height
$\Delta z$ where the cylinder is submerged in the AN layer, and the layer of
aqueous subphase that supports the active layer. We justify below that, given
the large viscosity contrast between the AN layer and the aqueous phase and
the used oil, the direct contribution of the AN is clearly dominant within the
force balance.

The force exerted by the AN layer on the cylindrical region of thickness
$\Delta z$ is computer as

$\begin{split}F^{AN}_{i}&=\iint_{S_{AN}}\sum_{j=1,2}\sigma^{AN}_{ij}n_{j}dS=\oint\sum_{j=1,2}\sigma^{AN}_{ij}n_{j}dl\;\Delta
z\\\ &=\oint\sum_{j=1,2}\sigma^{AN,2d}_{ij}n_{j}dl,\end{split}$ (5)

where the rightmost integration is performed along the circular edge of the
cylinder free surface. We implicitly assume that $\sigma^{AN}_{ij}$ is
invariant across the thickness of the AN layer. As a result, we have defined
the two-dimensional AN stress tensor,
$\sigma^{AN,2d}_{ij}=\sigma^{AN}_{ij}\Delta z$, with which we implement the
force balance without the need for a precise knowledge of the thickness of the
AN layer.

The 2D stress tensor is further decomposed in terms of the dominant viscous
and active stresses,
$\sigma^{AN,2d}_{ij}=\sigma_{ij}^{vis}+\sigma_{ij}^{act}$, respectively
defined as $\sigma_{ij}^{vis}=2\eta\varepsilon_{ij}$, and
$\sigma_{ij}^{act}=-\alpha Q_{ij}$. In these last expressions,
$\varepsilon_{ij}$ denotes the strain rate tensor, and $Q_{ij}$ is the
traceless nematic tensor order parameter, both in their two-dimensional forms.

Under the hypothesis of no-slip boundary conditions on the hydrogel column
surface, the inclusion of the latter inside the active material produces a
drastic change in the velocity field next to them. To accurately compute the
strain rate tensor it is thus necessary to achieve a substantial spatial
resolution of the velocity field around the obstacles. Since this condition is
not totally guaranteed from our PIV analysis, we resorted to estimating the
viscous drag using Lamb’s equation [45],

$F_{i}^{drag}=4\pi\epsilon\eta^{3D}v_{i}\Delta z=4\pi\epsilon\eta v_{i},$ (6)

in terms of the average velocity just before printing the pillars, and the
viscous drag coefficient, $\epsilon$. For the very low Reynolds number in our
experiments the latter is estimated as [45]

$\epsilon={[0.5-\gamma-\ln\left(\frac{Re}{8}\right)]}^{-1},$ (7)

where $Re$ is the Reynolds number, and $\gamma$ is the Euler constant. In our
experiments we took $\epsilon=0.052$, corresponding to an averaged value of
the Reynolds number $Re=3.2\times 10^{-8}$. The value of $L$ is considered to
be 40 $\mu$m (see above), while $r$ varied between 5.2 and 9 $\mu$m. For each
tested condition, experiments were repeated with different cylinder radii to
achieve an average deflection between $r$ and $2r$, approximately, which
resulted in the most reliable results for $\alpha$ and $\eta$.

For each column, one can thus write the force balance separating explicitly
the distinctive linear contributions in $\eta$ and $\alpha$ as,

$F_{i}^{def}=F_{i}^{drag}+F_{i}^{act}=F_{i}^{\eta}\eta+F_{i}^{\alpha}\alpha,$
(8)

with $F_{i}^{\eta}=4\pi\varepsilon v_{i}$ and
$F_{i}^{\alpha}=-\oint\sum_{j=1,2}Q_{ij}n_{j}dl$. This equation was applied
for each column, using data for the director and velocity fields measured at
the column location right before polymerization. In addition, we assumed that
the cylindrical cantilevers reach mechanical equilibrium when they change
noticeably their initial moving direction or when they reach a steady
deflection value.

For each column we obtained two equations, corresponding to the two components
of $\vec{F}^{def}$, as shown in Eq. 8. Therefore, if we analyze $N$ columns,
we will obtain $2N$ equations for the two unknowns $\eta$ and $\alpha$. This
overdetermined system of equations was solved numerically finding the two
parameters that minimize the sum of the squares of all the errors. In the
different sets of performed experiments the number of analyzed cylindrical
cantilevers varied between 37 and 118.

In this study, hydrogel columns where polymerized in a hexagonal lattice over
the field of view, with no a priori consideration on the local configuration
of the surrounding AN. This made some of the columns unsuitable for the
analysis described above. Columns that were discarded satisfied one or more of
the following conditions:

  * •

Incorrect tracking. The contour of columns that were polymerized on top of AN
defects could not be tracked with enough precision, as those regions are
deprived of fluorescent microtubules. This makes them indistinguishable from
hydrogel columns.

  * •

Regions with aligned filaments. In highly aligned AN regions, the Q tensor is
practically uniform in space, and velocity tracking using PIV software cannot
be performed with enough precision.

  * •

Highly distorted regions. In the area of influence of more than one AN defect,
the Q tensor and the velocity change considerably over small distances. This
made force quantification in those regions too sensitive to the precise
location of cylinder boundaries, resulting in imprecise computation of the
viscous and active force components.

###  II.2 Drag force from the oil and water phases

Self-sustained motion of the AN layer leads to flows on both the contacting
oil layer and the underlying aqueous phase. Both fluid phases exerted drag
forces on the polymerized cylindrical columns. In this section, we quantify
the resulting forces and demonstrate that they can be discarded when compared
to the forces exerted by the AN layer.

The oil phase can only produce drag on the free surface of the column. This
drag force can be calculated using the classical result for drag on a disk
[45] that, assuming no-slip boundary conditions at the AN-oil interface, reads

$F_{oil}^{drag}=\frac{1}{2}(\frac{32}{3}\eta_{oil}rv).$ (9)

We compare this expression with Eq. 6 to estimate the relative relevance of
the drag by the oil layer,

$\frac{F^{drag}}{F_{oil}^{drag}}=\frac{3\pi\epsilon}{4}\frac{\eta}{\eta_{oil}r}.$
(10)

The oil used in the experiments has a viscosity of 1.24$\times 10^{-3}$ Pa s,
and cylinders are a few microns in diameter. With $\eta$ of the order of
$10^{-6}$ Pa s m (see Fig. 2D main text), and $\epsilon$ a numerical factor of
order $10^{-2}$, the numerator in Eq. 10 is around two orders of magnitude
larger than the denominator, which indicates that drag due to the oil layer
can be safely neglected.

On the other hand, for the drag force exerted by the water phase all along the
column, we take the simplifying assumption that there is a linear velocity
profile that falls to zero on the substrate where columns are attached,

$v_{water}(z)=v\frac{z}{L}.$ (11)

We, then, evaluate the deflection under the non-uniform load that results from
applying Lamb’s equation for the viscous drag to the above velocity profile.
Integrating with the same boundary conditions indicated earlier, we obtain

$\delta_{water}=\frac{8}{15}\frac{\epsilon\eta_{water}vL^{4}}{Er^{4}},$ (12)

and, by comparison with the deflection due to the AN layer,

$\frac{\delta}{\delta_{w}}=\frac{\frac{4}{3\pi}\frac{F^{drag}L^{3}}{Er^{4}}}{\frac{8}{15}\frac{\epsilon\eta_{water}vL^{4}}{Er^{4}}}=\frac{5}{2\pi}\frac{4\pi\epsilon\eta
v}{\epsilon\eta_{water}vL}=10\frac{\eta}{\eta_{water}L},$ (13)

we find the contribution of the aqueous phase to the deflection of the
hydrogel cylinders to be less than $1\%$ of the contribution by the AN layer.
Hence, the effect of the aqueous subphase can be neglected.

###  II.3 In-situ determination of the hydrogel Young’s modulus

Figure S2: Calibration of the hydrogel rigidity. (A) Columns with
$r=5.4\;\mu$m are imprinted in a standard cell filled with M2B buffer and the
precursors of the hydrogel. The circular contour of the columns is highligtd
with black dashed lines in the brightfield images (B) The unpolymerized buffer
is replaced with M2B buffer containing methylcellulose and fluorescent
microparticles. Columns bend due to viscous drag. Their contour is highlighted
with yellow dashed lines. (C) When flow stops, columns relax to their initial
shape. Seeded microparticles can be observed in this panel. White dotted lines
indicate the center of the circular cross section of the original columns at
different times. Scale bar is 50 $\mu$m.

Calibration experiments employing the same protocol for hydrogel
photopolymerization were performed for an in situ determination of its Young’s
modulus (Fig. S2). For this purpose, we measured the deflection of cylindrical
columns due to the viscous drag of an aqueous fluid of know viscosity that
flows in the cell at a known velocity.

In this case, deflection is purely due to the viscous drag from the resultant
contributions of the aqueous and oil phases present in the standard flow cell.
The aqueous phase exerts its drag all along the body of the column. For
simplicity, we will assume that a uniform average force per unit length $q$ is
applied along the height of each column. We resort once again to Lamb’s
equation to obtain $q=4\pi\epsilon\eta_{water}v$, where $v$ will be the
average flow speed in the aqueous phase. On the other hand, the contribution
from the oil phase is simply expressed in terms of the corresponding drag on
the free column surface, as given in Eq. 9. Thus, the final expression for the
total bending moment reads,

$M(z)=q\frac{z^{2}}{2}+F_{oil}^{drag}z.$ (14)

Integrating the differential equation for the deflection (Eq. 1), we arrive at
an explicit equation for the deflection written as,

$\delta=\frac{3qL^{4}+8F_{oil}^{drag}L^{3}}{6\pi Er^{4}}.$ (15)

Different experiments were carried out with aqueous methylcellulose solutions,
spanning a range of viscosities between 3.6 and 18.7 mPa s, flowing at
velocities between 5 and 27$\,\mu$m s-1. For each experiment we evaluated the
corresponding Reynolds-number-dependent drag coefficient $\epsilon$. The final
estimation for the Young’s modulus evaluated this way is $E=18.1\pm 1.5$ Pa.

###  II.4 Active nematic viscosity estimated from kinetic energy spectra

In earlier work, Martínez-Prat et al. (Nature Physics 15, 362(2019)) estimated
the AN viscosity by fitting an ad hoc model for the kinetic energy spectrum.
In that work, boundary conditions for the oil layer in contact with the AN
where different from the ones we have here. While Martínez-Prat et al. used an
open cell and, thus, no-slip boundary conditions were assumed at the oil/air
interface, here our cells are closed, so we need to invoke no-slip boundary
conditions when both the water subphase and the oil layer contact with the
confining glass plates. This leads to a modified equation for the kinetic
energy spectrum,

$E(q)=\frac{Bq{R_{*}}^{4}e^{-q^{2}{R_{*}}^{2}/2}[I_{o}\left(\frac{q^{2}{R_{*}}^{2}}{2}\right)-I_{1}\left(\frac{q^{2}{R_{*}}^{2}}{2}\right)]}{{[q+\frac{\eta_{oil}}{\eta}{\coth}\left(qH_{oil}\right)+\frac{\eta_{w}}{\eta}\coth{\left(qH_{w}\right)]}}^{2}}.$
(16)

In this equation $q$ is the wave number, $\eta_{oil}$, $\eta_{w}$, and $\eta$
are the viscosities of the oil, water, and the active nematic layer; $H_{oil}$
and $H_{w}$ are the thickness of the oil and water phases, respectively;
$R^{*}$ is the mean vortex radius of the AN flow, $I_{0}$ and $I_{1}$ are
modified Bessel functions; $B=N\omega_{v}^{2}/(32\pi^{3}\mathcal{A)}$ is a
prefactor related to the total enstrophy, where N is the average number of
vortices over a total system area $\mathcal{A}$, and $\omega_{v}$ is the
average vorticity.

As a check for the values of $\eta$ obtained in this work using direct
rheological measurements, we have computed the kinetic energy spectrum for one
experiment with an average speed around $1\;\mu$m/s and we have fitted the
data to Eq. 16 (Fig. S4), yielding an estimation for the AN viscosity
$\eta=(1.3\pm 0.1)10^{-6}$ Pa s m. This value is of the same order as the one
we obtain in this work using in-situ microrheology.

SUPPLEMENTARY FIGURES

Figure S3: (A) Calibration of the average AN speed as a function of ATP
concentration. The solid line is a fit to a Michaelis-Menten model with fitted
parameters $v_{max}=3.6\pm 0.5\mu$m/s and $K_{M}=(3.7\pm 2)10^{2}\mu$M. (B)
Measured activity and shear viscosity of the AN as a function of ATP
concentration. Error bars are the standard deviation of the mean for $v$,
confidence bands for the fitted parameters $\alpha$ and $\eta$, and estimated
confidence band from experimental observations for [ATP].  Figure S4: Kinetic
energy spectrum of the AN layer with [ATP] = 450 $\mu$M. The solid line is a
fit to a variant of the model proposed by Martínez-Prat et al. that yield an
estimation for the shear viscosity of the active material $1.3\;10^{-6}$ Pa s
m. Figure S5: (A) Active length scale, $l_{\alpha}$, estimated from the defect
density for different ATP concentrations. (B) Estimated values for the elastic
constant of the active nematic, $K$, as a function of ATP concentration
obtained from $l_{\alpha}$ and $\alpha$. Error bars are estimated confidence
band from experimental observations for [ATP], and propagated errors from
$l_{\alpha}$ and $\alpha$ for $K$. Figure S6: Maps of the force density
surrounding a $+1/2$ defect in an AN layer with ATP = $50\pm 10$ $\mu$M.
Different panels represent the total AN force density (A), the active
component of the force density (B), and the viscous component of the force
density (C). For comparison, the map of the active flow velocity is included
in panel D. The length of the arrows corresponds to the speed, as indicated by
the speed scale. The scale bar is $50\;\mu$m.

SUPPLEMENTARY MOVIES

Movie S1: A set of four rigid rings that are in-situ polymerized in the AN
cell. Scalebar, $50\;\mu$m.  Movie S2: A flexible ring that has been in-situ
polymerized in the AN cell. Scale bar, $50\;\mu$m.  Movie S3: A hybrid
combination of a rigid ring boundary with embedded soft columnar pillars, all
built around a $+1/2$ defect. Scale bar, $50\;\mu$m.  Movie S4: An array of
soft pillars distributed in a hexagonal lattice to probe the microrheology of
the AN layer. Scale bar, $50\;\mu$m.  Movie S5: Force density map during the
unbinding of a pair of complementary defects. The black ellipse is a guide for
the eye to locate the region where the unbinding event takes place. Scale bar
is 50 $\mu$m.  Movie S6: Force density map during the annihilation of a pair
of complementary defects. The black ellipse is a guide for the eye to locate
the region where the event takes place. Scale bar is 50 $\mu$m.  Movie S7:
Hybrid microfabrication of a compliant column surrounded by a rigid ring. In
both cases the ring confines a $+1/2$ defect. Top: centered column. Bottom:
off-centered column. Scale bar is 50 $\mu$m.  Movie S8: Hybrid
microfabrication of a set of compliant columns surrounded by a rigid ring that
confines a $-1/2$ defect. Scale bar is 50 $\mu$m.


# Secure Antenna Selection and Beamforming in MIMO Systems

Chongjun Ouyang, Hao Xu, Xujie Zang, and Hongwen Yang Beijing University of
Posts and Telecommunications, Beijing, 100876, China
Email: {DragonAim, Xu_Hao, zangxj<EMAIL_ADDRESS>
###### Abstract

This work proposes a novel joint design for multiuser multiple-input multiple-
output wiretap channels. The base station exploits a switching network to
connect a subset of its antennas to the available radio frequency chains. The
switching network and transmit beamformers are jointly designed to maximize
the weighted secrecy sum-rate for this setting. The principal design problem
reduces to an NP-hard mixed-integer non-linear programming. We invoke the
fractional programming technique and the penalty dual decomposition method to
develop a tractable iterative algorithm that effectively approximates the
optimal design. Our numerical investigations validate the effectiveness of the
proposed algorithm and its superior performance compared with the benchmark.

###### Index Terms:

Antenna selection, multiuser multiple-input multiple-output, physical layer
security.

##  I Introduction

Information security has always been a critical issue for wireless
communications due to the broadcast property of the wireless medium. This fact
spawned the development and application of physical layer (PHY) security,
where the transmitter exploits secrecy channel coding to ensure perfect
security, i.e., eavesdroppers cannot decipher confidential information from
wiretapped messages. A feasible method to further boost the secrecy
performance at the PHY is to exploit the spatial degrees of freedom offered by
multiple antennas [1]. However, its fully digital implementation with a
dedicated radio frequency (RF) chain at each antenna suffers from expensive
hardware costs and excessive energy consumption. To address these challenges,
numerous potential technologies have been introduced over the last few years.
Among them is antenna selection, a technique to set only a small subset of
antennas active in each coherence time [2]. This technique can alleviate the
requirement on the number of RF transceivers without significantly sacrificing
the secrecy performance [3].

The past years have seen increasingly rapid advances in antenna selection
algorithm design in multiple-input multiple-output (MIMO) wiretap channels.
Most initial efforts focused on the secrecy performance achieved by single
antenna selection; see [2] and the references therein. Extension to multiple
antenna selection settings with single legitimate and eavesdropping user
terminals (UTs) were discussed in [3, 4, 5]. These works were further
generalized to the multiuser case [6, 7]. Despite extending the basic studies
on this topic, the available literature is still restricted to some particular
scenarios. Specifically, the results in [6, 7] were based on the maximal ratio
transmission (MRT) protocol, and a joint beamforming and antenna selection
design is still lacking.

###  I-A Contributions

This paper studies the secure transmission in a multiuser MIMO wiretap channel
where the base station (BS) sets only a subset of available antennas active
for communications. We maximize the secrecy throughput by proposing a joint
design that optimizes the switching network and beamformers at the BS. The
principal design problem belongs to mixed-integer non-linear programming
(MINLP). We address this challenging problem through the following
contributions: 1) We propose a penalty dual decomposition (PDD)-based method
to tackle the non-convex joint design problem via capitalizing on the
fractional programming (FP) technique. 2) To further reduce the computational
complexity, we propose an alternative algorithm based on sequential
optimization (SO). We show through numerical experiments that this algorithm
involves reduced complexity at the expense of a minor secrecy performance
loss. Our numerical results verify the capability of the proposed approaches
to outperform the benchmark significantly.

###  I-B Notation

Throughout this paper, scalars, vectors, and matrices are denoted by non-bold,
bold lower-case, and bold upper-case letters, respectively. The notations
$[\mathbf{A}]_{i,j}$ and $[\mathbf{a}]_{i}$ denote the $(i,j)$th entry and the
$i$th entry of matrix $\mathbf{A}$ and vector $\mathbf{a}$, respectively. The
identity matrix, zero matrix, and all-one vector are represented by
$\mathbf{I}$, $\mathbf{0}$, and $\mathbf{1}$, respectively. The Hadamard
product is shown by $\odot$ and $[K]$ represents the integer set
$\\{1,\ldots,K\\}$.

##  II System Model and Problem Formulation

Consider MIMO secure transmission where one $M$-antenna BS sends messages
simultaneously to $K$ legitimate UTs in the same time/frequency resources,
while these messages are confidential to $J$ adversaries that are treated as
eavesdroppers (Eves), as depicted in Figure 1. We assume that each UT
$k\in[K]$ and each Eve $j\in[J]$ are equipped with a single antenna for
receiving. Let ${\mathbf{h}}_{k}\in{\mathbbmss{C}}^{M\times 1}$ and
${\mathbf{g}}_{j}\in{\mathbbmss{C}}^{M\times 1}$ denote the UT $k$-to-BS and
Eve $j$-to-BS channel vectors, respectively. Furthermore, we denote
${x}_{k}\in{\mathbbmss{C}}$ as the securely coded data symbol dedicated to UT
$k$, with zero mean and unit variance. The data symbol $x_{k}$ is independent
of the symbols dedicated to other UTs, i.e.,
$\mathbbmss{E}\left\\{x_{k}x_{k^{\prime}}^{\mathsf{H}}\right\\}={0}$ for
$k\neq k^{\prime}$.

Figure 1: Illustration of a multiuser MIMO wiretap channel

The system operates in the time-division duplexing (TDD) mode, where the
instantaneous channel state information (CSI) can be estimated via pilot
sequences in the uplink training phase. This work focuses on passive
eavesdropping, where the eavesdroppers are registered UTs in the system but
distrusted by the legitimate UTs. In such a scenario, the eavesdroppers will
participate in the uplink training phase; thus, the BS can estimate their CSI.
We assume that the pilots are mutually orthogonal and that the estimation
error is negligible. The BS thus learns perfectly the CSI
$\\{{\mathbf{h}}_{k}\\}_{k=1}^{K}$ and $\\{{\mathbf{g}}_{j}\\}_{j=1}^{J}$.

###  II-A Transmit Antenna Selection

The BS has $N<M$ RF chains and thus uses a switching network to select a
subset of transmit antennas. This switching network connects the selected
antennas to the available $N$ RF chains at the BS. As a result, the
observations at UT $k$ and Eve $j$ are, respectively, given by

$\displaystyle
y_{k}={\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}\sum\nolimits_{k=1}^{K}\mathbf{w}_{k}x_{k}+n_{k},$
(1)
$\displaystyle\overline{y}_{j}={\mathbf{g}}_{j}^{\mathsf{H}}{\bm{\Delta}}\sum\nolimits_{k=1}^{K}\mathbf{w}_{k}x_{k}+\overline{n}_{j},$
(2)

where ${\mathbf{w}}_{k}\in{\mathbbmss{C}}^{M\times 1}$ is the beamforming
vector associated with symbol $x_{k}$,
$n_{k}\sim{\mathcal{CN}}(0,\sigma_{k}^{2})$ and
$\overline{n}_{j}\sim{\mathcal{CN}}(0,\delta_{j}^{2})$ denote the additive
white Gaussian noises (AWGNs) with $\sigma_{k}^{2}$ and $\delta_{j}^{2}$ being
the noise powers. Moreover,
${\bm{\Delta}}={\mathsf{diag}}\\{s_{1},\ldots,s_{M}\\}$ is the antenna
selection matrix with
$[s_{1},\ldots,s_{M}]^{\mathsf{T}}\triangleq{\mathbf{s}}$ and

$\displaystyle s_{m}=\begin{cases}1&{\text{antenna}}~{}m~{}{\text{is
selected}}\\\ 0&{\text{otherwise}}\end{cases}.$ (3)

We note that ${\bm{\Delta}}{\bm{\Delta}}^{\mathsf{H}}={\bm{\Delta}}$. Assume
that each UT $k$ and each Eve $j$ have access to instantaneous CSI of their
own channels with properly designed pilot signals. The signal-to-interference-
plus-noise ratio (SINR) at UT $k$ is given by

$\displaystyle\gamma_{k}=\frac{\lvert{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k}\rvert^{2}}{\sigma_{k}^{2}+\sum_{k^{\prime}\neq
k}\lvert{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k^{\prime}}\rvert^{2}}.$
(4)

###  II-B Performance Metric: Weighted Secrecy Sum-Rate

From a worst-case design perspective, we assume that all eavesdroppers
cooperate to overhear the secure transmission and can cancel out the
interference of other legitimate UTs. Given a beamforming matrix
${\mathbf{W}}=[{\mathbf{w}}_{1},\ldots,{\mathbf{w}}_{K}]$ and a selection
matrix ${\bm{\Delta}}$, the maximum secrecy rate to UT $k\in[K]$ in the worst-
case scenario is given by

$\displaystyle{\mathcal{R}}_{k}=[\log_{2}\left(1+\gamma_{k}\right)-\log_{2}\left(1+{\overline{\gamma}}_{k}\right),0]^{+},$
(5)

where $[\cdot]^{+}\triangleq\max\\{\cdot,0\\}$ and
${\overline{\gamma}}_{k}=\sum_{j=1}^{J}\frac{1}{\delta_{j}^{2}}\lvert{\mathbf{g}}_{j}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k}\rvert^{2}$
denotes the aggregated signal-to-noise ratio (SNR) at the eavesdroppers. To
quantify the secrecy throughput of the system, we define the weighted secrecy
sum-rate (WSSR) as

$\displaystyle{\mathcal{R}}=\sum\nolimits_{k=1}^{K}w_{k}[\log_{2}\left(1+\gamma_{k}\right)-\log_{2}\left(1+{\overline{\gamma}}_{k}\right),0]^{+}$
(6)

for non-negative weights $\\{w_{k}\\}_{k=1}^{K}$ corresponding to desired
quality-of-services (QoSs) of UTs.

###  II-C Problem Formulation

Our ultimate goal is to find the system design that optimizes the secrecy
throughput. This means that we strive to jointly design the beamforming matrix
$\mathbf{W}$ and the antenna selection matrix $\bm{\Delta}$, such that the
secrecy sum-rate term $\mathcal{R}$ is maximized. Consequently, our design
problem is formulated as

$\begin{split}\max_{{\mathbf{W}},{\bm{\Delta}}}~{}&{\mathcal{R}}=\sum\nolimits_{k=1}^{K}w_{k}{\mathcal{R}}_{k}\\\
\quad{\rm{s.t.}}~{}&{\mathcal{C}}_{1}:{\mathsf{tr}}({\mathbf{W}}{\mathbf{W}}^{\mathsf{H}})=\sum\nolimits_{k=1}^{K}{\mathbf{w}}_{k}^{\mathsf{H}}{\mathbf{w}}_{k}\leq
p,\\\
&{\mathcal{C}}_{2}:s_{m}\in\\{0,1\\}~{}{\text{for}}~{}m\in[M],{\mathbf{1}}^{\mathsf{T}}{\mathbf{s}}=N,\end{split}$
(${\mathcal{P}}_{1}$)

where ${\mathcal{C}}_{1}$ represents the transmit power constraint. Problem
(${\mathcal{P}}_{1}$) is a non-convex optimization problem due to the non-
convexity of $\mathcal{R}$ with respect to $\bm{\Delta}$ and $\mathbf{W}$. In
addition, considering the discrete constraints in ${\mathcal{C}}_{2}$,
(${\mathcal{P}}_{1}$) is an MINLP problem, which is difficult to solve
directly and to derive a globally optimal solution. In the sequel, we develop
an efficient framework to approximate the optimal design via a feasible
computational complexity.

##  III Proposed Solution

By invoking the FP framework [8], we simplify problem (${\mathcal{P}}_{1}$) to
a more tractable yet equivalent form. We then handle the resulting equivalent
problem via the PDD [9] that efficiently addresses non-convex non-smooth
problems with coupling equality constraints.

###  III-A Reformulation of Problem (${\mathcal{P}}_{1}$)

The operator $[\cdot]^{+}$ in (6) makes the objective of problem
(${\mathcal{P}}_{1}$) intractable. To take off this operator, we convert
problem (${\mathcal{P}}_{1}$) into its equivalent form as follows.

######  Lemma 1.

Define ${\mathbf{b}}=[b_{1},\ldots,b_{K}]^{\mathsf{T}}$. Problem
(${\mathcal{P}}_{1}$) is equivalent to the problem defined as follows:

$\begin{split}\max_{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}}}~{}&\overline{\mathcal{R}}=\sum\nolimits_{k=1}^{K}b_{k}(\log_{2}\left(1\\!+\\!\gamma_{k}\right)-\log_{2}\left(1\\!+\\!{\overline{\gamma}}_{k}\right))\\\
\quad{\rm{s.t.}}~{}&{\mathcal{C}}_{1},{\mathcal{C}}_{2},{\mathcal{C}}_{3}:b_{k}\in[0,w_{k}],~{}{\text{for}}~{}k\in[K].\end{split}$
(${\mathcal{P}}_{2}$)

This means that the solutions $\\{{\mathbf{W}},{\bm{\Delta}}\\}$ for both the
problems are identical.

###### Proof:

Please refer to Appendix -1 for more details. ∎

We further rewrite the objective of (${\mathcal{P}}_{2}$) as

$\displaystyle\overline{\mathcal{R}}=\sum\nolimits_{k=1}^{K}b_{k}\left(\log_{2}\left(1+\gamma_{k}\right)+\log_{2}\left(\frac{1}{1+{\overline{\gamma}}_{k}}\right)\right).$
(7)

This _sum-of-functions-of-ratio_ form motivates us to exploit the FP framework
to further simplify (${\mathcal{P}}_{2}$). We start the derivations by the
following lemma:

######  Lemma 2.

Define
$g\triangleq\sum_{j=1}^{J}\frac{p}{\delta_{j}^{2}}\lVert{\mathbf{g}}_{j}\rVert^{2}$.
Problem (${\mathcal{P}}_{2}$) is equivalent to the problem defined as follows:

$\begin{split}\max_{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}}}&\hat{\mathcal{R}}=\sum\nolimits_{k=1}^{K}\\!\\!b_{k}\\!\\!\left(\\!\log_{2}\left(1\\!+\\!\gamma_{k}\right)\\!+\\!\log_{2}\left(\\!1\\!+\\!\frac{g\\!-\\!{\overline{\gamma}}_{k}}{1\\!+\\!{\overline{\gamma}}_{k}}\\!\right)\\!\\!\right)\\\
\quad{\rm{s.t.}}~{}&{\mathcal{C}}_{1},{\mathcal{C}}_{2},{\mathcal{C}}_{3},\end{split}$
(${\mathcal{P}}_{3}$)

where $g-{\overline{\gamma}}_{k}\geq 0$.

###### Proof:

Please refer to Appendix -2 for more details. ∎

We next exploit the FP framework to convert the fractional programming problem
(${\mathcal{P}}_{3}$) equivalently to a new problem which involves a more
tractable objective function.

######  Lemma 3.

Problem (${\mathcal{P}}_{3}$) is equivalent to:

$\begin{split}\max_{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}}}&{\mathcal{F}}({\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}})=\sum\nolimits_{k=1}^{K}b_{k}(f_{1}^{k}+f_{2}^{k})\\\
\quad{\rm{s.t.}}~{}&{\mathcal{C}}_{1},{\mathcal{C}}_{2},{\mathcal{C}}_{3},\end{split}$
(${\mathcal{P}}_{4}$)

where ${\bm{\alpha}}=[\alpha_{1},\ldots,\alpha_{K}]^{\mathsf{T}}$,
${\bm{\beta}}=[\beta_{1},\ldots,\beta_{K}]^{\mathsf{T}}$,

$\displaystyle
f_{1}^{k}=\log(1+\alpha_{k})-\alpha_{k}+\frac{(1+\alpha_{k})\lvert{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k}\rvert^{2}}{\sum_{i=1}^{K}\lvert{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{i}\rvert^{2}+\sigma_{k}^{2}},$
(8) $\displaystyle
f_{2}^{k}=\log(1+\beta_{k})-\beta_{k}+{(1+\beta_{k})(g-{\overline{\gamma}}_{k})}({1+g})^{-1}.$
(9)

The optimal $\alpha_{k}$ and $\beta_{k}$ are given by
$\alpha_{k}^{\star}=\gamma_{k}$ and
$\beta_{k}^{\star}=\frac{g-{\overline{\gamma}}_{k}}{1+{\overline{\gamma}}_{k}}$,
respectively.

###### Proof:

Please refer to Appendix -3 for more details. ∎

It is worth noting that $g$ is independent of
$\\{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}}\\}$. Therefore, the intractability
of
${\mathcal{F}}({\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}})$
mainly originates from the fractional terms involved in
$\\{f_{1}^{k}\\}_{k=1}^{K}$. To handle this difficulty, we introduce auxiliary
variables ${\bm{\eta}}=[\eta_{1},\ldots,\eta_{K}]^{\mathsf{T}}$ to simplify
${\mathcal{F}}({\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}})$.
Then, the following lemma can be found.

######  Lemma 4.

Problem (${\mathcal{P}}_{4}$) is equivalent to:

$\begin{split}\max_{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}}}&{\mathcal{L}}({\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}})\\!=\\!\sum\nolimits_{k=1}^{K}\\!b_{k}(g_{1}^{k}\\!+\\!f_{2}^{k})\\\
\quad{\rm{s.t.}}~{}&{\mathcal{C}}_{1},{\mathcal{C}}_{2},{\mathcal{C}}_{3},\end{split}$
(${\mathcal{P}}_{5}$)

where
$g_{1,k}=\log(1+\alpha_{k})-\alpha_{k}+(1+\alpha_{k}){\overline{g}}_{1}^{k}$
with

$\begin{split}{\overline{g}}_{1}^{k}=2\Re\\{{\eta}_{k}^{*}{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k}\\}\\!-\\!|\eta_{k}|^{2}\big{(}\sum\nolimits_{i=1}^{K}\\!|{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{i}|^{2}\\!+\\!\sigma_{k}^{2}\big{)}.\end{split}$
(10)

The optimal $\eta_{k}$ is given by
$\eta_{k}^{\star}=(\sum\nolimits_{i=1}^{K}\\!|{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{i}|^{2}\\!+\\!\sigma_{k}^{2})^{-1}{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k}$.

###### Proof:

Similar to the proof of Lemma 3. ∎

######  Remark 1.

Lemmas 1 to 4 indicate that (${\mathcal{P}}_{5}$) is a variational form of the
principal problem (${\mathcal{P}}_{1}$). This means that the solutions
$\\{{\mathbf{W}},{\bm{\Delta}}\\}$ for both the problems are identical. It is
worth noting that unlike the original form in problem (${\mathcal{P}}_{1}$),
the problem (${\mathcal{P}}_{5}$) has an objective function that is marginally
convex over each of the variable ${\mathbf{W}}$, ${\bm{\Delta}}$,
${\mathbf{b}}$, ${\bm{\alpha}}$, ${\bm{\beta}}$, and ${\bm{\eta}}$.

Furthermore, to tackle the discrete constraints in ${\mathcal{C}}_{2}$, we
define the auxiliary variables
$\overline{\mathbf{s}}=[\overline{s}_{1},\ldots,\overline{s}_{M}]^{\mathsf{T}}$
which satisfy the following constraints: ${\overline{s}}_{m}=s_{m}$ and
$s_{m}\left(1-{\overline{s}}_{m}\right)=0$. We thus can equivalently find the
solution of (${\mathcal{P}}_{1}$) by solving the following optimization:

$\begin{split}&\max_{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}},\overline{\mathbf{s}}}{\mathcal{L}}({\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}})\\\
&\quad{\rm{s.t.}}~{}{\mathcal{C}}_{1},{\mathcal{C}}_{2},{\mathcal{C}}_{3},\\\
&\qquad~{}~{}{\mathcal{C}}_{4}:{\overline{s}}_{m}=s_{m},s_{m}\left(1-{\overline{s}}_{m}\right)=0,~{}{\text{for}}~{}m\in[M].\end{split}$
(${\mathcal{P}}_{6}$)

We note that in (${\mathcal{P}}_{6}$), variables $\mathbf{s}$ and
$\overline{\mathbf{s}}$ are only constrained through
${\overline{s}}_{m}=s_{m}$, $s_{m}\left(1-{\overline{s}}_{m}\right)=0$, and
${\mathbf{1}}^{\mathsf{T}}{\mathbf{s}}=N$. To handle these equality
constraints, we resort to the PDD technique.

###  III-B The Proposed PDD-Based Algorithm

The PDD-based algorithm is characterized by an embedded double loop structure
[9]. The inner loop solves the augmented Lagrangian (AL) subproblem while the
outer loop updates the dual variables and the penalty parameters that
correspond to constraint violation. We obtain the AL problem by moving the
equality constraints as a penalty term to the objective function. By [9], the
AL problem corresponding to (${\mathcal{P}}_{6}$) is given by

$\begin{split}\max_{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}},\overline{\mathbf{s}}}{\mathcal{L}}({\mathbf{W}},{\bm{\Delta}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}})-f_{\rho}\quad{\rm{s.t.}}~{}{\mathcal{C}}_{1},{\mathcal{C}}_{3}.\end{split}$
(${\mathcal{P}}_{7}$)

where $\rho>0$ is the penalty parameter penalizing the violation of the
equality constraints, and $f_{\rho}$ is given by

$\displaystyle f_{\rho}$
$\displaystyle=\frac{1}{2\rho}\left[\left({\mathbf{1}}^{\mathsf{T}}{\mathbf{s}}\\!-\\!N\\!+\\!\rho\xi\right)^{2}+\sum\nolimits_{m=1}^{M}\left[\left(s_{m}\\!-\\!{\overline{s}}_{m}\\!+\\!\rho\mu_{m}\right)^{2}\right.\right.$
$\displaystyle\left.\left.+\left(s_{m}\left(1-{\overline{s}}_{m}\right)+\rho\lambda_{m}\right)^{2}\right]\right],$
(11)

with $\xi$, ${\bm{\mu}}=[\mu_{1},\ldots,\mu_{M}]^{\mathsf{T}}$, and
${\bm{\lambda}}=[\lambda_{1},\ldots,\lambda_{M}]^{\mathsf{T}}$ denoting the
Lagrangian dual variables associated with the equality constraints in
${\mathbf{1}}^{\mathsf{T}}{\mathbf{s}}=T$, ${\overline{s}}_{m}=s_{m}$, and
$s_{m}\left(1-{\overline{s}}_{m}\right)=0$, respectively. It is observed that
as $\rho\rightarrow 0$, the penalty term is forced to zero, i.e., equality
constraints are enforced. It is shown in [9] that updating the primal and dual
variables, as well as the penalty factor in an alternating manner, PDD
converges to a stationary-point solution. Thus, in the following, we focus on
solving problem (${\mathcal{P}}_{7}$) by invoking the block coordinate descent
(BCD) method. More specifically, we sequentially update ${\mathbf{W}}$,
$\bm{\Delta}$ (or ${\mathbf{s}}$), $\overline{\mathbf{s}}$, ${\mathbf{b}}$,
$\bm{\alpha}$, $\bm{\beta}$, and $\bm{\eta}$ in the inner loop through the
following marginal optimizations.

####  III-B1 Optimizing $\mathbf{W}$

The marginal problem for $\mathbf{W}$ reads

$\begin{split}\min_{\mathbf{W}}~{}\sum\nolimits_{k=1}^{K}({\mathbf{w}}_{k}^{\mathsf{H}}{\mathbf{A}}_{k}{\mathbf{w}}_{k}-2\Re\\{{\mathbf{w}}_{k}^{\mathsf{H}}{\mathbf{a}}_{k}\\})\quad{\rm{s.t.}}~{}{\mathcal{C}}_{1}.\end{split}$
(${\mathcal{P}}_{8}$)

where
${\mathbf{a}}_{k}=b_{k}(1+\alpha_{k}){\eta}_{k}{\bm{\Delta}}^{\mathsf{H}}{\mathbf{h}}_{k}$
and

$\begin{split}{\mathbf{A}}_{k}&=\sum\nolimits_{i=1}^{K}b_{i}(1+\alpha_{i})|\eta_{i}|^{2}{\bm{\Delta}}^{\mathsf{H}}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{\mathsf{H}}{\bm{\Delta}}\\\
&+{b_{k}(1+\beta_{k})}{({1+g})^{-1}}\sum\nolimits_{j=1}^{J}{\bm{\Delta}}^{\mathsf{H}}{\mathbf{g}}_{j}{\mathbf{g}}_{j}^{\mathsf{H}}{\bm{\Delta}}\succeq{\mathbf{0}}.\end{split}$
(12)

This is a standard convex quadratic optimization subproblem whose solution is
given by

$\displaystyle{\mathbf{w}}_{k}^{\star}=({\mathbf{A}}_{k}+\lambda{\mathbf{I}})^{-1}{\mathbf{a}}_{k},~{}{\text{for}}~{}k\in[K].$
(13)

The regularizer $\lambda$ is chosen, such that the complementarity slackness
condition, i.e.,
$\lambda({\mathsf{tr}}({\mathbf{W}}{\mathbf{W}}^{\mathsf{H}})-p)=0$, is
satisfied. If
$\sum_{k=1}^{K}{\mathbf{a}}_{k}^{\mathsf{H}}{\mathbf{A}}_{k}^{-2}{\mathbf{a}}_{k}=p$;
then, $\lambda=0$. Otherwise, we can obtain the solution of $\lambda$ from the
following identity:

$\displaystyle{\mathsf{tr}}({\mathbf{W}}{\mathbf{W}}^{\mathsf{H}})=\sum\nolimits_{k=1}^{K}{\mathbf{a}}_{k}^{\mathsf{H}}({\mathbf{A}}_{k}+\lambda{\mathbf{I}})^{-2}{\mathbf{a}}_{k}=p.$
(14)

It follows that

$\displaystyle\sum\nolimits_{k=1}^{K}\sum\nolimits_{m=1}^{M}\frac{|[{\mathbf{U}}_{{\mathbf{A}}_{k}}{\mathbf{a}}_{k}]_{m}|^{2}}{([{\bm{\Lambda}}_{{\mathbf{A}}_{k}}]_{m,m}+\lambda)^{2}}=p,$
(15)

where
${\mathbf{A}}_{k}={\mathbf{U}}_{{\mathbf{A}}_{k}}^{\mathsf{H}}{\bm{\Lambda}}_{{\mathbf{A}}_{k}}{\mathbf{U}}_{{\mathbf{A}}_{k}}$
is the eigen-decomposition of ${\mathbf{A}}_{k}$. Since
$[{\bm{\Lambda}}_{{\mathbf{A}}_{k}}]_{m,m}\geq 0$ for $m\in[M]$, the left-hand
side of (15) is a monotonously decreasing function of $\lambda\geq 0$.
Consequently, we can find $\lambda$ by solving equation (15) using the
bisection-based search method.

####  III-B2 Optimizing $\mathbf{s}$

Define ${\mathbf{W}}_{k}={\mathsf{diag}}\\{{\mathbf{w}}_{k}\\}$ for $k\in[K]$,

$\displaystyle{\mathbf{Q}}$
$\displaystyle\triangleq{\Re}\bigg{\\{}\sum\nolimits_{i=1}^{K}\sum\nolimits_{k=1}^{K}w_{k}(1+\alpha_{k})|\eta_{k}|^{2}{\mathbf{W}}_{i}^{\mathsf{H}}{\mathbf{h}}_{k}{\mathbf{h}}_{k}^{\mathsf{H}}{\mathbf{W}}_{i}$
$\displaystyle+\sum\nolimits_{k=1}^{K}\frac{b_{k}(1+\beta_{k})}{1+g}\sum\nolimits_{j=1}^{J}{\mathbf{W}}_{k}^{\mathsf{H}}{\mathbf{g}}_{j}{\mathbf{g}}_{j}^{\mathsf{H}}{\mathbf{W}}_{k}+\frac{1}{2\rho}$
$\displaystyle\times({\mathbf{1}}{\mathbf{1}}^{\mathsf{T}}+{\mathbf{I}}+{{\mathsf{diag}}\\{\left(1-{\overline{s}}_{1}\right)^{2},\cdots,\left(1-{\overline{s}}_{M}\right)^{2}\\}})\big{\\}},$
(16) $\displaystyle{\mathbf{g}}$ $\displaystyle\triangleq
2\Re\\!\left\\{{\mathbf{q}}\right\\}\\!-\\!\frac{1}{\rho}\left(\left(\rho\xi\\!-\\!N\right){\mathbf{1}}\\!+\\!\left(\rho{\bm{\mu}}\\!-\\!\overline{\mathbf{s}}\right)\\!+\\!\rho\left({\mathbf{1}}\\!-\\!\overline{\mathbf{s}}\right)\\!\odot\\!{\bm{\lambda}}\right),$
(17) $\displaystyle{\mathbf{q}}$
$\displaystyle\triangleq\sum\nolimits_{k=1}^{K}(1+\alpha_{k})\Re\\{{\eta}_{k}{\mathbf{W}}_{k}^{\mathsf{H}}{\mathbf{h}}_{k}\\},$
(18)

The marginal problem for $\mathbf{s}$ is given by

${\mathbf{s}}^{\star}=\operatorname*{argmin}\nolimits_{{\mathbf{s}}\in{\mathbbmss{R}}^{M}}\left({\mathbf{s}}^{\mathsf{T}}{\mathbf{Q}}{\mathbf{s}}-{\mathbf{s}}^{\mathsf{T}}{\mathbf{g}}\right).$
(19)

It is readily shown that ${\mathbf{Q}}\succ{\mathbf{0}}$, which means that the
problem (19) is a standard convex problem. Using its first-order optimality
condition, we obtain
${{\mathbf{s}}}^{\star}=\frac{1}{2}{\mathbf{Q}}^{-1}{\mathbf{g}}$.

####  III-B3 Optimizing $\overline{\mathbf{s}}$

The marginal optimization for $\\{\overline{s}_{m}\\}_{m=1}^{M}$ decouples
into $M$ scalar optimization problems with solution to the $m$th problem being

$\overline{s}_{m}^{\star}=\operatorname*{argmin}\nolimits_{\overline{s}_{m}\in{\mathbbmss{R}}}\left(\overline{a}_{m}\overline{s}_{m}^{2}-2\overline{b}_{m}\overline{s}_{m}\right).$
(20)

Here, $\overline{a}_{m}=1+s_{m}^{2}$ and
$\overline{b}_{m}=s_{m}+\rho\mu_{m}+s_{m}^{2}+s_{m}\rho\lambda_{m}$. The
optimization in (20) is a quadratic scalar problem whose solution is given by
${\overline{s}}_{m}^{\star}={\overline{b}_{m}}/{\overline{a}_{m}}$.

####  III-B4 Optimizing ${\mathbf{b}}$

The marginal optimization for $\\{b_{k}\\}_{k=1}^{K}$ decouples into $K$
scalar optimization problems with solution to the $k$th problem being

$b_{k}^{\star}=\operatorname*{argmax}\nolimits_{b_{k}\in[0,w_{k}]}\left(b_{k}(g_{1}^{k}\\!+\\!f_{2}^{k})\right)=w_{k}{\mathbf{1}}_{\\{g_{1}^{k}>-f_{2}^{k}\\}},$
(21)

where ${\mathbf{1}}_{\\{\cdot\\}}$ represents the indicator function.

####  III-B5 Optimizing the Auxiliary Variables

For given
$\\{{\mathbf{W}},{\mathbf{s}},\overline{\mathbf{s}},{\mathbf{b}}\\}$, the
auxiliary variables $\bm{\alpha}$, $\bm{\beta}$, and $\bm{\eta}$ can be
updated by Lemmas 3 and 4. The marginally optimal $\bm{\alpha}$, $\bm{\beta}$,
and $\bm{\eta}$ are given by $\alpha_{k}^{\star}=\gamma_{k}$,
$\beta_{k}^{\star}=\frac{g-{\overline{\gamma}}_{k}}{1+{\overline{\gamma}}_{k}}$,
and
$\eta_{k}^{\star}=(\sum\nolimits_{i=1}^{K}\\!|{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{i}|^{2}\\!+\\!\sigma_{k}^{2})^{-1}{\mathbf{h}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k}$,
respectively.

Finally, the dual variables $\\{{\xi},{\bm{\lambda}},{\bm{\mu}}\\}$ and
penalty factor $\rho$ can be updated following the PDD framework. As a
consequence, the overall algorithm for solving problem (${\mathcal{P}}_{1}$)
is summarized in Algorithm 1, where $h\triangleq\max_{\forall
m}\\{|{\mathbf{1}}\\!^{\mathsf{T}}{\mathbf{s}}\\!-\\!N|,\left|{\overline{s}}_{m}\\!-\\!s_{m}\right|,\left|s_{m}\\!\left(1\\!-\\!{\overline{s}}_{m}\right)\right|\\}$
denotes the constraint violation function. The computational complexity of the
proposed algorithm can be further characterized in terms of problem
dimensions. To this end, let $I_{\rm{out}}$ and $I_{\rm{fp}}$ denote the
numbers of iterations in the outer loop and the inner-loop FP, respectively.
It is readily shown that the complexity of marginal optimizations with respect
to ${\mathbf{W}}$, ${\mathbf{s}}$, $\overline{\mathbf{s}}$, ${\mathbf{b}}$,
$\bm{\alpha}$, $\bm{\beta}$, and $\bm{\eta}$ scale with
${\mathcal{O}}\left(M^{3}\right)$, ${\mathcal{O}}\left(M^{3}\right)$,
${\mathcal{O}}\left(M\right)$, ${\mathcal{O}}\left(KM^{2}\right)$,
${\mathcal{O}}\left(K\right)$, ${\mathcal{O}}\left(K\right)$, and
${\mathcal{O}}\left(K\right)$, respectively. Therefore, the overall complexity
of Algorithm 1 scales with
${\mathcal{O}}(I_{\text{out}}I_{\text{fp}}(2M^{3}+KM^{2}))$.

###  III-C Low-Complexity Sequential Optimization

Although the PDD-based algorithm leads to a feasible complexity, it can lead
to a computational burden in many applications. To this end, we develop an
alternative scheme based on sequential optimization approach. In this respect,
we first design the beamforming matrix using the FP-based method; then, we
select the antennas using a greedy search (GS). These two steps are
illustrated in the sequel.

####  III-C1 Transmit Precoding Design

At the first step of the SO-based method, we design the beamforming matrix
$\mathbf{W}$ by setting ${\bm{\Delta}}={\mathbf{I}}$. The resulting problem
can be solved by using the FP-based method. The detailed steps are similar to
those outlined in the previous part and are omitted here for brevity.

####  III-C2 Antenna Selection

In the next step, we optimize the antenna selection vector $\mathbf{s}$ for
the updated beamformer. The resulting problem of antenna selection is given by

$\displaystyle{\mathbf{s}}^{\star}=\operatorname*{argmax}{\mathcal{R}},~{}{\rm{s.t.}}~{}{\mathcal{C}}_{2}.$
(22)

To tackle the discrete constraint, we exploit the GS method. More details
about this method can be found in [10] and are omitted here due to page
limitations.

Based on [10], we can show that the complexity of sequential optimizations
with respect to ${\mathbf{s}}$ and $\mathbf{W}$ scale with
${\mathcal{O}}(KNM)$, and ${\mathcal{O}}(I_{\text{fp}}M^{3})$, respectively.
Hence, the complexity of the SO-based method scales with
${\mathcal{O}}(KNM+I_{\text{fp}}M^{3})$, which is lower than the PDD-based
method.

Algorithm 1 PDD-based algorithm for solving problem (${\mathcal{P}}_{1}$)

1: Initialize primary variables
$\left\\{{\mathbf{W}},{\mathbf{s}},\overline{\mathbf{s}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}}\right\\}$,
dual variables $\left\\{\xi,{\bm{\lambda}},{\bm{\mu}}\right\\}$, threshold
$\mu$, penalty factor $\rho>0$, and scaling factor $\chi\in(0,1)$

2: repeat

3: repeat

4: Update
$\left\\{{\mathbf{W}},{\mathbf{s}},\overline{\mathbf{s}},{\mathbf{b}},{\bm{\alpha}},{\bm{\beta}},{\bm{\eta}}\right\\}$
by the BCD method

5: until convergence

6: if $h<\mu$ then

7: Set $\xi=\xi+{\rho}^{-1}({\mathbf{1}}^{\mathsf{T}}{\mathbf{s}}-N)$,
${\bm{\mu}}={\bm{\mu}}+{\rho}^{-1}\left({\overline{\mathbf{s}}}-{\mathbf{s}}_{n}\right)$,
and
${\bm{\lambda}}={\bm{\lambda}}+{\rho}^{-1}{\mathbf{s}}\odot\left(1-{\overline{\mathbf{s}}}\right)$

8: else

9: Set $\rho=\chi\rho$

10: end if

11: Set $\mu=\chi h$ and $t=t+1$

12: until convergence

##  IV Numerical Results

In this section, simulation results are provided to verify the effectiveness
of the proposed algorithms. For simulations, the following parameters are used
unless explicitly mentioned otherwise: $M=24$, $K=6$, $J=4$, $w_{k}=1$,
$\log{\sigma_{k}^{2}}=\log{\delta_{j}^{2}}=-120$ dBm, $\mu=1$, $\rho=1$, and
$\chi=0.1$. We further generate the channel realizations as follows. Regarding
the large-scale fading, we assume that all the channels, i.e.,
$\\{{\mathbf{h}}_{k}\\}_{k=1}^{K}$ and $\\{{\mathbf{g}}_{j}\\}_{j=1}^{J}$,
exhibit the same path loss $-120$ dB for illustration. Meanwhile, for the
small-scale fading of the channels, we consider the standard Rayleigh fading
model. All the optimization variables are randomly initialized, and all
simulation curves are averaged over $500$ independent channel realizations.

(a) Convergence of the WSSR.

(b) Constraint violation.

Figure 2: Average convergence performances for $L=6$.

In Figure 2, we first study the convergence behavior of the proposed PDD-based
method. It can be seen from Figure 2(a) that the WSSR rapidly converges to a
stationary value. Figure 2(b) shows the constraint violation in terms of the
number of outer iterations. We observe that the constraint violation reduces
to a threshold $10^{-4}$ in less than $100$ outer iterations, meaning that the
solution has essentially met the equality constraints for problem
(${\mathcal{P}}_{1}$).

(a) WSSR vs. the power budget.

(b) WSSR vs. the RF chain number.

Figure 3: WSSR performance: (a) $L=6$, (b) $\log{p}=10$ dBm.

We next consider the following benchmark schemes for performance comparison:
1) Random scheme, in which $\mathbf{s}$ is randomly set, and ${\mathbf{W}}$ is
optimized by the FP-based method. 2) Energy-based scheme, in which the BS
antennas corresponding to the $L$ strongest channel gains are selected, and
${\mathbf{W}}$ is optimized by the FP-based method. This selection strategy
was utilized in previous works [3, 2]. 3) MRT-based scheme, in which
$\mathbf{W}$ is generated by the MRT protocol, and $\mathbf{s}$ is designed by
the GS-based method [6, 7].

Figure 3 compares the WSSR performances achieved by the proposed methods and
the benchmark methods. In Figure 3(a), the WSSR is plotted as a function of
the transmit power budget $p$ for different optimization schemes assuming
$L=6$. From this graph, we find that the proposed PDD-based method has the
best secrecy performance, followed by the SO-based scheme. As stated before,
the SO-based method involves less computational complexity than the PDD-based
scheme. The above results imply that the SO-based method is preferred for
systems with highly restricted computational capacity. What stands out in this
graph is that the methods using FP-based secure beamforming, i.e., the PDD-
based, SO-based, Energy-based, and Random schemes, are far superior to the
MRT-based scheme in terms of the WSSR, especially in the high-SNR regime. This
observation highlights the superiority of the proposed FP-based beamforming
design.

We next fix $\log{p}=10$ dBm and plot the WSSR versus the number of RF chains
for selected BS antennas. The results are shown in Figure 3(b) implying that
in all the presented cases, the increase in the number of RF chains improves
the secrecy performance. Specifically, the proposed two methods significantly
outperform the benchmarks in terms of the WSSR. It validates the advantages of
using our proposed joint design to enhance secure transmissions at the PHY.

##  V Conclusion

We proposed an iterative algorithm capitalizing on the FP and PDD techniques
for joint antenna selection and transmit beamforming in MIMO wiretap channels.
We also developed an alternative algorithm with reduced computational
complexity based on the SO scheme. Our numerical results imply that the PDD-
based method outperforms the SO-based and benchmark methods. The SO-based
algorithm closely tracks the performance of the PDD-based method while
benefiting considerably in terms of complexity.

####  -1 Proof of Lemma 1

We prove this lemma by showing that (${\mathcal{P}}_{2}$) can be equivalently
transformed to (${\mathcal{P}}_{1}$). Given
$\\{{\mathbf{W}},{\bm{\Delta}}\\}$, the optimal $b_{k}$ satisfies
$b_{k}^{\star}=w_{k}{\mathbf{1}}_{\\{\gamma_{k}>{\overline{\gamma}}_{k}\\}}$.
We note that substituting $b_{k}=b_{k}^{\star}$ into $\overline{\mathcal{R}}$
recovers the objective of (${\mathcal{P}}_{1}$):

$\displaystyle\overline{\mathcal{R}}$
$\displaystyle=\sum\nolimits_{k=1}^{K}{\mathbf{1}}_{\\{\gamma_{k}>{\overline{\gamma}}_{k}\\}}w_{k}(\log_{2}\left(1\\!+\\!\gamma_{k}\right)\\!-\\!\log_{2}\left(1\\!+\\!{\overline{\gamma}}_{k}\right))$
(23)
$\displaystyle=\sum\nolimits_{k=1}^{K}w_{k}[\log_{2}\left(1\\!+\\!\gamma_{k}\right)\\!-\\!\log_{2}\left(1\\!+\\!{\overline{\gamma}}_{k}\right),0]^{+}={\mathcal{R}}.$
(24)

The final results follow immediately.

####  -2 Proof of Lemma 2

By the Cauchy inequality, we have

$\displaystyle{\overline{\gamma}}_{k}$
$\displaystyle\leq\sum\nolimits_{j=1}^{J}{\delta_{j}^{-2}}\lVert{\mathbf{g}}_{j}\rVert^{2}\lVert{\bm{\Delta}}{\mathbf{w}}_{k}\rVert^{2}.$
(25)

We next exploit the following fact:

$\displaystyle\lVert{\bm{\Delta}}{\mathbf{w}}_{k}\rVert^{2}={\mathbf{w}}_{k}^{\mathsf{H}}{\bm{\Delta}}{\mathbf{w}}_{k}\leq{\mathbf{w}}_{k}^{\mathsf{H}}{\mathbf{I}}{\mathbf{w}}_{k}=\lVert{\mathbf{w}}_{k}\rVert^{2},$
(26)

to obtain the following relationship:

$\displaystyle{\overline{\gamma}}_{k}$
$\displaystyle\leq\sum\nolimits_{j=1}^{J}{\delta_{j}^{-2}}\lVert{\mathbf{g}}_{j}\rVert^{2}\lVert{\mathbf{w}}_{k}\rVert^{2}=g.$
(27)

It follows that

$\displaystyle\log_{2}\left(\frac{1}{1+{\overline{\gamma}}_{k}}\right)$
$\displaystyle=\log_{2}\left(\frac{1+g}{1+{\overline{\gamma}}_{k}}\right)-\log_{2}\left(1+g\right)$
(28)
$\displaystyle=\log_{2}\left(1+\frac{g-{\overline{\gamma}}_{k}}{1+{\overline{\gamma}}_{k}}\right)-\log_{2}\left(1+g\right).$
(29)

Note that $\log_{2}\left(1+g\right)$ is independent of
$\\{{\mathbf{W}},{\bm{\Delta}}\\}$. Hence, we can establish the equivalence
between (${\mathcal{P}}_{2}$) and (${\mathcal{P}}_{3}$).

####  -3 Proof of Lemma 3

Since $\mathcal{F}$ is concave over $\bm{\alpha}$ and $\bm{\beta}$ for fixed
$\\{{\mathbf{W}},{\bm{\Delta}},{\mathbf{b}}\\}$, we take its complex
derivative and solve each $\frac{\partial}{\partial\alpha_{k}}{\mathcal{F}}=0$
and $\frac{\partial}{\partial\beta_{k}}{\mathcal{F}}=0$. The optimal
$\alpha_{k}$ and $\beta_{k}$ are easily seen as $\alpha_{k}^{\star}$ and
$\beta_{k}^{\star}$, respectively. Inserting $\alpha_{k}^{\star}$ and
$\beta_{k}^{\star}$ back to $\mathcal{F}$ recovers the objective function in
(${\mathcal{P}}_{3}$), thus establishing the equivalence of these two
problems.

## References

  * [1] X. Chen, D. W. K. Ng, W. Gerstacker, and H.-H. Chen, “A survey on multiple-antenna techniques for physical layer security,” _IEEE Commun. Surveys Tuts._ , vol. 19, no. 2, pp. 1027–1053, Second Quarter, 2017.
  * [2] C. Ouyang _et al._ , “Security enhancement via antenna selection in MIMOME channels with discrete inputs,” _IEEE Trans. Commun._ , vol. 68, no. 8, pp. 5041–5055, Aug. 2020.
  * [3] S. Asaad _et al._ , “Optimal transmit antenna selection for massive MIMO wiretap channels,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 4, pp. 817–828, Apr. 2018.
  * [4] M. Tian _et al._ , “Joint beamforming design and receive antenna selection for large-scale MIMO wiretap channels,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 3, pp. 2716–2730, Mar. 2020.
  * [5] C. Ouyang, Z. Ou, L. Zhang, and H. Yang, “Optimal transmit antenna selection algorithm in massive MIMOME channels,” in _Proc. IEEE Wireless Commun. Netw. Conf. (WCNC)_ , pp. 1–6, 2019.
  * [6] A. Bereyhi _et al._ , “Securing massive MIMO systems: Secrecy for free with low-complexity architectures,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 9, pp. 5831–5845, Sep. 2021.
  * [7] A. Bereyhi _et al._ , “Iterative antenna selection for secrecy enhancement in massive MIMO wiretap channels,” in _Proc. IEEE Int. Wkshp. Signal Process. Adv. Wireless Commun. (SPAWC)_ , pp. 1–5, 2018.
  * [8] K. Shen and W. Yu, “Fractional programming for communication systems–Part I: Power control and beamforming,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2616–2630, May. 2018.
  * [9] Q. Shi and M. Hong, “Penalty dual decomposition method for nonsmooth nonconvex optimization–Part I: Algorithms and convergence analysis,” _IEEE Trans. Signal Process._ , vol. 68, pp. 4108–4122, 2020.
  * [10] M. Gharavi-Alkhansari and A. B. Gershman, “Fast antenna subset selection in MIMO systems,” _IEEE Trans. Signal Process._ , vol. 52, no. 2, pp. 339–347, Feb. 2004.