# Emergence of nuclear saturation within $\Delta$–full chiral effective field
theory

W. G. Jiang  Department of Physics, Chalmers University of Technology, SE-412
96 Göteborg, Sweden  Institut für Kernphysik and PRISMA Cluster of Excellence,
Johannes Gutenberg Universität, 55128 Mainz, Germany  C. Forssén  Department
of Physics, Chalmers University of Technology, SE-412 96 Göteborg, Sweden  T.
Djärv  Department of Physics, Chalmers University of Technology, SE-412 96
Göteborg, Sweden  Physics Division, Oak Ridge National Laboratory, Oak Ridge,
TN 37831, USA  G. Hagen  Physics Division, Oak Ridge National Laboratory, Oak
Ridge, TN 37831, USA  Department of Physics and Astronomy, University of
Tennessee, Knoxville, TN 37996, USA

###### Abstract

Nuclear saturation is a key property of low-energy nuclear physics that
depends on fine details of the nuclear interaction. Here we develop a unified
statistical framework that uses realistic nuclear forces to link the
theoretical modeling of finite nuclei and infinite nuclear matter. We
construct fast and accurate emulators for nuclear-matter observables and
employ an iterative history-matching approach to explore and reduce the
enormous parameter domain of $\Delta$-full chiral interactions. We find that
model calibration including 16O observables gives saturation predictions that
are more accurate and have a smaller variance than those that only use few-
body data.

Introduction—A key question in low-energy nuclear physics is whether it is
possible to successfully describe all systems from finite nuclei to infinite
nuclear matter using nucleons as effective degrees of freedom. Realistic
interaction models based on chiral effective field theory ($\chi$EFT) used in
combination with _ab initio_ methods, that solve the many-body problem with
controlled approximations, have the potential to deliver on that research
program Machleidt and Entem (2011); Hebeler  _et al._ (2011); Ekström  _et
al._ (2013, 2015); Jiang  _et al._ (2020); Somà  _et al._ (2020); Hu  _et al._
(2022); Maris  _et al._ (2022); Elhatisari  _et al._ (2022). However, using a
complicated interaction model and calibration data with limited independence
implies a risk for overfitting unless relevant theoretical uncertainties are
accounted for. Furthermore, the extrapolation from well-studied few-nucleon
systems to heavier nuclei and infinite nuclear matter will lead to increasing
variances in model predictions. Therefore it becomes important to study the
precision of these predictions and to explore the sensitivity to different
choices of calibration data. The emergence of nuclear saturation—represented
by a minimum in the equation-of-state for infinite symmetric nuclear matter
(SNM)—is particularly important since it affects bulk properties (such as
binding energies and radii) of finite atomic nuclei.

In the last decade significant progress has been made towards quantifying EFT
uncertainties of _ab initio_ nuclear matter predictions and identifying
possible correlations with observables in finite nuclei Hebeler  _et al._
(2011); Khan and Margueron (2012); Hagen  _et al._ (2016); Kievsky  _et al._
(2018); Drischler  _et al._ (2019); Reinhard and Nazarewicz (2016); Carson
_et al._ (2019); Drischler  _et al._ (2020a, b); Hu  _et al._ (2022); Huth
_et al._ (2022). Very recently, some of the authors of this work identified
correlations between the symmetry energy and its slope with the neutron skin
and dipole polarizability of the heavy nucleus 208Pb starting from chiral
interactions at next-to-next-to-leading order (NNLO) with explicit delta
isobars ($\Delta$) Hu  _et al._ (2022). This finding was made possible by
employing novel statistical tools such as bayesian inference von Toussaint
(2011) and history matching Vernon  _et al._ (2010, 2014, 2018); Hu  _et al._
(2022) together with accurate emulators of _ab initio_ computations of light
nuclei Ekström and Hagen (2019); König  _et al._ (2020); Wesolowski  _et al._
(2021a); Djärv  _et al._ (2022).

In this Letter we develop fast and accurate emulators of coupled-cluster
computations of pure neutron matter (PNM) and SNM Baardsen  _et al._ (2013);
Hagen  _et al._ (2014a) starting from $\Delta$-full $\chi$EFT interaction
models at NNLO van Kolck (1994); Hemmert  _et al._ (1998); Kaiser  _et al._
(1998); Krebs  _et al._ (2007); Epelbaum  _et al._ (2008); Ekström  _et al._
(2018); Jiang  _et al._ (2020). Our construction of emulators is based on the
subspace-projected coupled cluster (SPCC) method Ekström and Hagen (2019),
which we here extend with small batch voting Jiang  _et al._ (2022)—a new on-
the-fly validation approach that addresses the possible appearance of spurious
states when diagonalizing a Hamiltonian in a subspace of bi-orthogonal
coupled-cluster solutions Shavitt and Bartlett (2009). Using these emulators
we can reproduce full-space coupled-cluster computations of infinite nuclear
matter with high precision at a tiny fraction of the computational cost.

After validating our approach, we use history matching to identify a region of
the 17-dimensional parameter space of the low-energy constants (LECs) at NNLO
that give acceptable results when confronted with few-body data. This
iterative parameter search is enabled by employing emulators of few-body
systems up to 4He and allows us to collect $1.7\times 10^{6}$ non-implausible
interaction samples. Finally, we calibrate our _ab initio_ model with two
alternative LEC posterior probability density functions (PDFs) and use the
principle of importance resampling Smith and Gelfand (1992); Jiang and Forssén
(2022) to quantify the sensitivity of nuclear matter predictions to
calibration data. Here we focus on the results of this analysis, and the
emergence of nuclear saturation, while details of the small-batch voting
scheme and the history match are presented in a companion paper Jiang  _et
al._ (2022).

Method—The chiral Hamiltonian of $\Delta$NNLO is parametrized with 17 LECs
Jiang  _et al._ (2020), and following Refs. Ekström and Hagen (2019); König
_et al._ (2020); Wesolowski  _et al._ (2021a) it can be written as

$\displaystyle
H(\alpha)=h_{0}+\displaystyle\sum_{i=1}^{N_{\rm{LECs}}=17}\alpha_{i}h_{i}.$
(1)

Here, $h_{0}=t_{\rm{kin}}+v_{0}$ where $t_{\rm{kin}}$ is the kinetic energy
and $v_{0}$ represents the constant potential term without LEC-dependence.
Note that ${\vec{\alpha}}$ is a vector that denotes all LECs. In this work we
use non-local regulators with a cutoff $\Lambda=394$ MeV/$c$.

Recently, model reduction methods Melendez  _et al._ (2022) such as
eigenvector continuation (EC) Frame  _et al._ (2018); Sarkar and Lee (2021),
has proven to be both efficient and accurate for emulating the output of _ab
initio_ computations of both scattering Melendez  _et al._ (2021); Zhang and
Furnstahl (2022) and bound-state observables Ekström and Hagen (2019); König
_et al._ (2020); Wesolowski  _et al._ (2021a); Djärv  _et al._ (2022); Demol
_et al._ (2020). These methods employ the fact that the eigenvector trajectory
generated by smooth changes of the Hamiltonian matrix is well approximated by
a very low-dimensional manifold Frame  _et al._ (2018). Following Refs.
Ekström and Hagen (2019); König  _et al._ (2020) we can therefore obtain good
approximations of the ground-state of a target Hamiltonian
$H({\vec{\alpha}}_{\circledcirc})$ by projecting it on a subspace of
$N_{\rm{sub}}$ different ground-state training eigenvectors and solving the
corresponding $N_{\rm{sub}}\times N_{\rm{sub}}$-dimensional generalized
eigenvalue problem. In this work we employ emulators for 2,3H and 4He for
history matching, and the SPCC method with singles-, doubles-, and leading
order triples excitations (CCSDT-3) for subsequent model calibration on 16O
(see Ref. Hu  _et al._ (2022) for details). We also construct a new emulator
for 6Li based on Ref. Djärv  _et al._ (2022) for model checking. Here we used
$N_{\mathrm{sub}}=32$ training points reaching a relative accuracy of
$10^{-3}$ for the non-implausible parametrizations Djärv (2021).

In this Letter we use the coupled-cluster method Coester (1958); Coester and
Kümmel (1960); Čížek (1966); Kümmel  _et al._ (1978); Bishop (1991); Zeng  _et
al._ (1998); Mihaila and Heisenberg (2000); Dean and Hjorth-Jensen (2004);
Bartlett and Musiał (2007); Shavitt and Bartlett (2009); Hagen  _et al._
(2010); Binder  _et al._ (2013); Hagen  _et al._ (2014b) within the doubles
approximation (CCD) and solve for pure neutron and symmetric nuclear matter on
a cubic lattice in momentum space using periodic boundary conditions Hagen
_et al._ (2014a). The model-space has $(2n_{\rm{max}}+1)^{3}$ momentum points
with $n_{\rm{max}}=4$. We use 132 nucleons for SNM and 66 neutrons for PNM
which allows to minimize finite-size effects Gandolfi  _et al._ (2009); Hagen
_et al._ (2014a). We then use the EC-inspired SPCC method Ekström and Hagen
(2019) to construct subspace emulators at five different nuclear densities.
Since the subspace projected coupled-cluster Hamiltionian is non-Hermitian,
the variational theorem does not hold Teschl (2009). This might lead to the
appearance of “spurious” states in the subspace spectrum that are lower in
energy than the corresponding full-space coupled-cluster solution for certain
combinations of the LECs. These states appear due to the non-symmetric
treatment of the left and right bi-orthogonal coupled-cluster states Bartlett
and Musiał (2007). Therefore, we have developed a new algorithm—called small-
batch voting—which allows to identify the physical ground state using a set of
different subspace projections. This algorithm employs the strong sensitivity
of spurious states on the specific choice of basis, and the contrasting
stability of physical states. Details of this algorithm are given in the
companion paper Jiang  _et al._ (2022).

Results—We used history matching to collect a large number of $1.7\times
10^{6}$ model parametrizations that exhibit acceptable (or at least non-
implausible) fits to the calibration data. In the final wave we defined the
non-implausible volume using a rotated hyperrectangle, as illustrated in Fig.
1, that captured parameter correlations and allowed to significantly increase
the number of collected non-implausible samples. Further details of the
history match are presented in the companion paper Jiang  _et al._ (2022). We
then used a Latin hypercube design of 64 samples from the wave 5 volume for
which we performed full CCD computations for SNM and PNM at five different
densities $\rho\in\\{0.12,0.14,0.16,0.18,0.20\\}$ fm-3. Together with small-
batch voting this allowed to create SPCC emulators for $E/A$ and $E/N$ at
these densities. Finally, the nuclear matter equation-of-state (EOS) around
saturation could be obtained by interpolation using Gaussian processes C. E.
Rasmussen (2006) as described in more detail in Ref. Jiang  _et al._ (2022).

Figure 1: (Color online) Iterative parameter volume reduction via history
matching shown here for the strongly correlated LEC pairs $(c_{D},c_{E})$ and
$(C_{3S1},\tilde{C}_{3S1})$, where the former are dimensionless and the latter
are in units of $10^{4}$ GeV-2 and $10^{4}$ GeV-4, respectively. The initial
parameter domain used at the start of history-matching is represented by the
axes limits while iteratively reduced non-implausible volumes at the start of
waves 2, 3, 4, and 5 are indicated by green/dash-dotted, blue/dashed,
black/dotted, red/solid parallelograms. Note that $A=3,4$ observables (that
constrain $c_{D}$ and $c_{E}$) do not enter until the fourth wave. The final
set of $1.7\times 10^{6}$ non-implausible samples is shown as two-dimensional
histograms.

Cross-validation of emulator performance is shown in Fig. 2. We conclude that
our emulators predict nuclear matter properties with $<5\%$ precision at a
computational cost that is six orders of magnitude smaller than the full
solution.

Figure 2: (Color online) Cross-validation of the SPCC emulator with exact CCD
calculations using 50 interaction samples. Shown here from top to bottom is
the nuclear matter saturation density $(\rho_{0})$, saturation energy
($E_{0}/A)$, and symmetry energy $(S)$. Histograms of relative errors are
shown in the right column.

Having access to fast and accurate emulators, and having identified a non-
implausible region in the parameter space of our chiral interaction model, we
can proceed to study the general behaviour of nuclear matter model predictions
and possible correlations between different properties of these systems. Here
we are not interested in the usual optimization approach that results in a
single (“best fit”) model parametrization. Instead, we consider all non-
implausible samples from the history match. About 65% of these samples predict
saturation within the studied density region and are kept for further
consideration. The outputs from the nuclear matter emulators for these
$1.1\times 10^{6}$ interaction samples are shown in the upper triangle of Fig.
3. Here we have applied a density-dependent energy shift to approximately
account for triples corrections Jiang  _et al._ (2022).

We observe a strong anti-correlation between the saturation energy $E_{0}/A$
and the saturation density $\rho_{0}$ (Pearson correlation coefficient
$r=-0.90$). Similarly, the symmetry energy $S$ and its slope $L$ show a
positive correlation ($r=0.71$). These correlations have been seen in DFT
calculations Ducoin  _et al._ (2011); Santos  _et al._ (2014) and have been
indicated with small families of EFT-inspired Kievsky  _et al._ (2018) or
chiral interactions Drischler  _et al._ (2019). The comparison with empirical
nuclear matter properties reveals that the model predictions from the history
match are clustered in a region with too small $|E_{0}/A|$, $\rho_{0}$, and
$S$. This finding is consistent with previous results from various few-body
optimized interactions Coraggio  _et al._ (2014); Ekström  _et al._ (2018);
Drischler  _et al._ (2020b).

The non-implausible predictions of nuclear matter properties that is shown in
the upper triangle of Fig. 3 is a form of Bayes linear forecasting. By just
considering expectation values such as means and variances—thereby avoiding
the full probabilistic specification of uncertain quantities—we made the
analysis simpler and more technically straightforward. While this allowed us
to identify the parameter region of interest and to explore the model’s
forecasting capabilities, we now make an effort to extract a posterior
predictive distribution (PPD) with which we can make probabilistic statements.

Figure 3: (Color online) Upper triangle: nuclear matter emulator output
(saturation density $\rho_{0}$, saturation energy $E_{0}/A$, symmetry energy
$S$, symmetry energy slope $L$ and incompressibility $K$) for the non-
implausible interactions from the fifth wave of history matching. The axes
limits are the same as in the corresponding panels in the lower triangle.
Lower triangle and diagonal: PPDs for nuclear matter properties using two
different PDFs for the LECs plus error sampling. These predictions are based
either on few-body ($A=2,3,4$) calibration data (orange PPD) or the addition
of 16O to the calibration data set (blue PPD). See the Supplemental Material
sup  for a Gaussian approximation of ${\rm PPD}_{A=2,3,4,16}$. The contour
lines on the bivariate distributions denote $68\%$ and $90\%$ credible
regions. The red bands indicate empirical ranges for the observables
($E_{0}/A=-16.0\pm 0.5$ MeV, $\rho_{0}=0.16\pm 0.01$ fm-3, $S=31\pm 1$,
$L=50\pm 10$ and $K=240\pm 20$ MeV from Refs. Lattimer and Lim (2013); Bender
_et al._ (2003); Shlomo  _et al._ (2006)).

First, we revisit all remaining interaction parametrizations and keep only
those that survive an extra implausibility analysis involving selected phase
shifts in $S$ and $P$ waves up to $T_{\mathrm{lab}}=200$ MeV and that give an
unbound ${}^{1}S_{0}$ ground-state in the $np$ system. This step results in a
significant reduction to $10^{4}$ acceptable samples. Second, we introduce a
set of calibration observables $\mathcal{D}_{\mathrm{cal}}$ and a normally-
distributed likelihood,
$\mathcal{L}\left(\mathcal{D}_{\mathrm{cal}}\,|\,{\vec{\alpha}}\right)$,
assuming independent experimental, method, emulator, and model errors Jiang
_et al._ (2022). At this stage we employ the established method of
sampling/importance resampling Smith and Gelfand (1992); Jiang and Forssén
(2022) to approximately extract samples from the parameter posterior
$\mathrm{pr}\left({\vec{\alpha}}\,|\,\mathcal{D}_{\mathrm{cal}},I\right)$. In
this step we assume a uniform prior probability for all non-implausible
samples 111The prior for $c_{1,2,3,4}$ is the multivariate Gaussian resulting
from a Roy-Steiner analysis of $\pi N$ scattering data Siemens  _et al._
(2017)..

In this work we especially considered two different sets of model calibration
data: (i) $\mathcal{D}_{\mathrm{cal}}=\mathcal{D}_{A=2,3,4}$ encompassing
binding energies and point-proton radii of 2,3H and 4He plus the quadrupole
moment of the deuteron, and (ii)
$\mathcal{D}_{\mathrm{cal}}=\mathcal{D}_{A=2,3,4,16}$ where we complemented
the above set of observables with the energy and radius of 16O.

Having access to approximate parameter posteriors we first extract samples
from the model PPD defined as ${\rm
PPD}_{\mathrm{th}}=\\{{\rm\textbf{y}_{th}}({{\vec{\alpha}}}):{{\vec{\alpha}}}\sim\mathrm{pr}\left({\vec{\alpha}}\,|\,\mathcal{D}_{\mathrm{cal}},I\right)\\}$.
As a non-trivial model validation we predict the 6Li threshold energy
$S_{d}=E({}^{2}{\text{H}})+E({}^{4}{\text{He}})-E({}^{6}{\text{Li}})$ using
our EC emulators. The 90% credible intervals are $[-0.14,1.54]$ MeV and
$[1.23,1.84]$ MeV for ${\rm PPD}_{\mathrm{th}}$ with $\mathcal{D}_{A=2,3,4}$
and $\mathcal{D}_{A=2,3,4,16}$ calibration data, respectively. Both
predictions are consistent with the experimental value $S_{d}=1.474$ MeV
Tilley  _et al._ (2002), but the latter one is significantly more precise.

We then predict nuclear matter properties. For these predictions we collect
samples from the full PPD, which includes the EFT truncation error, CC method
errors, and the SPCC emulator error (estimated from cross validation shown in
Fig. 2). The density dependence and relevant cross correlation of these errors
are described by a Bayesian machine learning error model Drischler  _et al._
(2020a, b); Jiang  _et al._ (2022). Using the two different sets of
calibration data, $\mathcal{D}_{A=2,3,4}$ and $\mathcal{D}_{A=2,3,4,16}$, we
compare the resulting PPDs that we label ${\rm PPD}_{A=2,3,4}$ and ${\rm
PPD}_{A=2,3,4,16}$, respectively. These are shown in the lower triangle of
Fig. 3.

The marginal distributions on the diagonal reveal that $\rho_{0}$, $E_{0}/A$
and $S$ for ${\rm PPD}_{A=2,3,4}$ still deviate significantly from the
empirical region—as previously observed with the non-importance-weighted
samples. In other words, even by enforcing a higher accuracy for two- and few-
body systems via the data likelihood, the description of nuclear matter
properties is not improved. On the other hand, the predictions become more
accurate when we include the 16O observables. In particular, the saturation
point is much improved. We actually find that ${\rm PPD}_{A=2,3,4}$ displays a
significant asymmetry in some directions, e.g., $\rho_{0}$, which hints to a
possible bimodality. Samples from the tail region of ${\rm PPD}_{A=2,3,4}$
correspond largely to the mode of ${\rm PPD}_{A=2,3,4,16}$. This finding shows
that the emergence of saturation, represented by the position and shape of
this mode, can be altered depending on the choice of calibration observables.

Summary and outlook—In this work we constructed emulators that accurately
reproduce full-space coupled-cluster computations of nucleonic matter starting
from $\Delta$-full $\chi$EFT at NNLO while reducing the computational cost by
several orders of magnitude. The small-batch voting algorithm was developed to
remove spurious states that occured in the SPCC method. Using these tools,
together with emulators of light nuclei, we employed history matching to
identify more than one million acceptable interaction samples, and could
reveal correlations between different properties of infinite nuclear matter.
We then employed importance resampling to obtain PPDs, and studied the
sensitivity to the choice of calibration data. Broad, asymmetric marginal
distributions for the saturation energy and density were observed when making
predictions conditional on few-nucleon data only, while both accuracy and
precision were much improved when adding the 16O energy and radius to the
calibration data set. Future work using the history-matching framework and
many-body emulators can help to elucidate the information content of various
low-energy observables and the variances of _ab initio_ model predictions
across the nuclear landscape. Binding energies and radii of finite nuclei such
as 16O are obviously useful for model calibration, but are not the only
choice. Other observables such as the 3H beta-decay rate and three-nucleon
scattering cross sections should also be considered within a statistical
framework involving many-body predictions. Fortunately there is ongoing work
that shows that also these observables can be efficiently computed Miller  _et
al._ (2022a, b) or emulated Wesolowski  _et al._ (2021b); Zhang and Furnstahl
(2022).

###### Acknowledgements.

We thank Andreas Ekström and Thomas Papenbrock for useful discussions. This
work was supported by the Swedish Research Council (Grant Nos 2017-04234 and
2021-04507), the European Research Council under the European Unions Horizon
2020 research and innovation program (Grant No. 758027), and the U.S.
Department of Energy under contract DE-AC05-00OR22725 with UT-Battelle, LLC
(Oak Ridge National Laboratory). The computations and data handling were
enabled by resources provided by the Swedish National Infrastructure for
Computing (SNIC) at Chalmers Centre for Computational Science and Engineering
(C3SE), and the National Supercomputer Centre (NSC) partially funded by the
Swedish Research Council through Grant No. 2018-05973.

## References

  * Machleidt and Entem (2011) R. Machleidt and D.R. Entem, “Chiral effective field theory and nuclear forces,” Phys. Rep. 503, 1 – 75 (2011).EOS
  * Hebeler  _et al._ (2011) K. Hebeler, S. K. Bogner, R. J. Furnstahl, A. Nogga, and A. Schwenk, “Improved nuclear matter calculations from chiral low-momentum interactions,” Phys. Rev. C 83, 031301 (2011).EOS
  * Ekström  _et al._ (2013) A. Ekström, G. Baardsen, C. Forssén, G. Hagen, M. Hjorth-Jensen, G. R. Jansen, R. Machleidt, W. Nazarewicz, T. Papenbrock, J. Sarich, and S. M. Wild, “Optimized chiral nucleon-nucleon interaction at next-to-next-to-leading order,” Phys. Rev. Lett. 110, 192502 (2013).EOS
  * Ekström  _et al._ (2015) A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, “Accurate nuclear radii and binding energies from a chiral interaction,” Phys. Rev. C 91, 051301 (2015).EOS
  * Jiang  _et al._ (2020) W. G. Jiang, A. Ekström, C. Forssén, G. Hagen, G. R. Jansen, and T. Papenbrock, “Accurate bulk properties of nuclei from $A=2$ to $\infty$ from potentials with $\Delta$ isobars,” Phys. Rev. C 102, 054301 (2020), arXiv:2006.16774 [nucl-th] .EOS
  * Somà  _et al._ (2020) V. Somà, P. Navrátil, F. Raimondi, C. Barbieri, and T. Duguet, “Novel chiral hamiltonian and observables in light and medium-mass nuclei,” Phys. Rev. C 101, 014318 (2020).EOS
  * Hu  _et al._ (2022) B. S. Hu, W. G. Jiang, T. Miyagi, Z. H. Sun, A. Ekström, C. Forssén, G. Hagen, J. D. Holt, T. Papenbrock, S. R. Stroberg, and I. Vernon, “Ab initio predictions link the neutron skin of 208Pb to nuclear forces,” Nature Phys. 18, 1196–1200 (2022), arXiv:2112.01125 [nucl-th] .EOS
  * Maris  _et al._ (2022) P. Maris _et al._ (LENPIC), “Nuclear properties with semilocal momentum-space regularized chiral interactions beyond N2LO,” Phys. Rev. C 106, 064002 (2022), arXiv:2206.13303 [nucl-th] .EOS
  * Elhatisari  _et al._ (2022) Serdar Elhatisari, Lukas Bovermann, Evgeny Epelbaum, Dillon Frame, Fabian Hildenbrand, Hermann Krebs, Timo A. Lähde, Dean Lee, Ning Li, Bing-Nan Lu, Myungkuk Kim, Youngman Kim, Yuanzhuo Ma, Ulf-G. Meißner, Gautam Rupak, Shihang Shen, Young-Ho Song, and Gianluca Stellin, “Wave function matching for the quantum many-body problem,” arXiv e-prints , arXiv:2210.17488 (2022), arXiv:2210.17488 [nucl-th] .EOS
  * Khan and Margueron (2012) E. Khan and J. Margueron, “Constraining the nuclear equation of state at subsaturation densities,” Phys. Rev. Lett. 109, 092501 (2012), arXiv:1204.0399 [nucl-th] .EOS
  * Hagen  _et al._ (2016) G. Hagen, A. Ekström, C. Forssén, G. R. Jansen, W. Nazarewicz, T. Papenbrock, K. A. Wendt, S. Bacca, N. Barnea, B. Carlsson, C. Drischler, K. Hebeler, M. Hjorth-Jensen, M. Miorelli, G. Orlandini, A. Schwenk, and J. Simonis, “Neutron and weak-charge distributions of the 48Ca nucleus,” Nature Phys. 12, 186 (2016).EOS
  * Kievsky  _et al._ (2018) A. Kievsky, M. Viviani, D. Logoteta, I. Bombaci, and L. Girlanda, “Correlations imposed by the unitary limit between few-nucleon systems, nuclear matter and neutron stars,” Phys. Rev. Lett. 121, 072701 (2018), arXiv:1806.02636 [nucl-th] .EOS
  * Drischler  _et al._ (2019) C. Drischler, K. Hebeler, and A. Schwenk, “Chiral interactions up to next-to-next-to-next-to-leading order and nuclear saturation,” Phys. Rev. Lett. 122, 042501 (2019).EOS
  * Reinhard and Nazarewicz (2016) P. G. Reinhard and W. Nazarewicz, “Nuclear charge and neutron radii and nuclear matter: Trend analysis in Skyrme density-functional-theory approach,” Phys. Rev. C 93, 051303 (2016), arXiv:1601.06324 [nucl-th] .EOS
  * Carson  _et al._ (2019) Z. Carson, A. W. Steiner, and K. Yagi, “Constraining nuclear matter parameters with GW170817,” Phys. Rev. D 99, 043010 (2019), arXiv:1812.08910 [gr-qc] .EOS
  * Drischler  _et al._ (2020a) C Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, “How Well Do We Know the Neutron-Matter Equation of State at the Densities Inside Neutron Stars? A Bayesian Approach with Correlated Uncertainties,” Phys. Rev. Lett. 125, 202702 (2020a), 2004.07232 .EOS
  * Drischler  _et al._ (2020b) C. Drischler, J. A. Melendez, R. J. Furnstahl, and D. R. Phillips, “Quantifying uncertainties and correlations in the nuclear-matter equation of state,” Phys. Rev. C 102, 054315 (2020b), arXiv:2004.07805 [nucl-th] .EOS
  * Huth  _et al._ (2022) Sabrina Huth, Peter T. H. Pang, Ingo Tews, Tim Dietrich, Arnaud Le Fèvre, Achim Schwenk, Wolfgang Trautmann, Kshitij Agarwal, Mattia Bulla, Michael W. Coughlin, and Chris Van Den Broeck, “Constraining neutron-star matter with microscopic and macroscopic collisions,” Nature 606, 276–280 (2022).EOS
  * von Toussaint (2011) Udo von Toussaint, “Bayesian inference in physics,” Rev. Mod. Phys. 83, 943–999 (2011).EOS
  * Vernon  _et al._ (2010) I. Vernon, M. Goldstein, and R. G. Bower, “Galaxy formation: a bayesian uncertainty analysis,” Bayesian Anal. 5, 619–669 (2010).EOS
  * Vernon  _et al._ (2014) I. Vernon, M. Goldstein, and R. Bower, “Galaxy formation: Bayesian history matching for the observable universe,” Statist. Sci. 29, 81–90 (2014).EOS
  * Vernon  _et al._ (2018) I. Vernon, J. Liu, M. Goldstein, J. Rowe, J. Topping, and K. Lindsey, “Bayesian uncertainty analysis for complex systems biology models: emulation, global parameter searches and evaluation of gene functions,” BMC Systems Biology 12, 1 (2018).EOS
  * Ekström and Hagen (2019) A. Ekström and G. Hagen, “Global sensitivity analysis of bulk properties of an atomic nucleus,” Phys. Rev. Lett. 123, 252501 (2019).EOS
  * König  _et al._ (2020) S. König, A. Ekström, K. Hebeler, D. Lee, and A. Schwenk, “Eigenvector Continuation as an Efficient and Accurate Emulator for Uncertainty Quantification,” Phys. Lett. B 810, 135814 (2020), arXiv:1909.08446 [nucl-th] .EOS
  * Wesolowski  _et al._ (2021a) S. Wesolowski, I. Svensson, A. Ekström, C. Forssén, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, “Rigorous constraints on three-nucleon forces in chiral effective field theory from fast and accurate calculations of few-body observables,” Phys. Rev. C 104, 064001 (2021a).EOS
  * Djärv  _et al._ (2022) T. Djärv, A. Ekström, C. Forssén, and H. T. Johansson, “Bayesian predictions for A=6 nuclei using eigenvector continuation emulators,” Phys. Rev. C 105, 014005 (2022), arXiv:2108.13313 [nucl-th] .EOS
  * Baardsen  _et al._ (2013) G. Baardsen, A. Ekström, G. Hagen, and M. Hjorth-Jensen, “Coupled-cluster studies of infinite nuclear matter,” Phys. Rev. C 88, 054312 (2013).EOS
  * Hagen  _et al._ (2014a) G. Hagen, T. Papenbrock, A. Ekström, K. A. Wendt, G. Baardsen, S. Gandolfi, M. Hjorth-Jensen, and C. J. Horowitz, “Coupled-cluster calculations of nucleonic matter,” Phys. Rev. C 89, 014319 (2014a).EOS
  * van Kolck (1994) U. van Kolck, “Few-nucleon forces from chiral Lagrangians,” Phys. Rev. C 49, 2932–2941 (1994).EOS
  * Hemmert  _et al._ (1998) T. R. Hemmert, B. R. Holstein, and J. Kambor, “Heavy baryon chiral perturbation theory with light deltas,” J. Phys. G 24, 1831 (1998).EOS
  * Kaiser  _et al._ (1998) N. Kaiser, S. Gerstendörfer, and W. Weise, “Peripheral nn-scattering: role of delta-excitation, correlated two-pion and vector meson exchange,” Nucl. Phys. A 637, 395 – 420 (1998).EOS
  * Krebs  _et al._ (2007) H. Krebs, E. Epelbaum, and U. G. Meißner, “Nuclear forces with $\mathrm{\Delta}$ excitations up to next-to-next-to-leading order, part i: Peripheral nucleon-nucleon waves,” Eur. Phys. J. A 32, 127–137 (2007).EOS
  * Epelbaum  _et al._ (2008) E. Epelbaum, H. Krebs, and U.-G. Meißner, “$\Delta$-excitations and the three-nucleon force,” Nucl. Phys. A 806, 65 – 78 (2008).EOS
  * Ekström  _et al._ (2018) A. Ekström, G. Hagen, T. D. Morris, T. Papenbrock, and P. D. Schwartz, “$\Delta$ isobars and nuclear saturation,” Phys. Rev. C 97, 024332 (2018), arXiv:1707.09028 [nucl-th] .EOS
  * Jiang  _et al._ (2022) W. G. Jiang, C. Forssén, T. Djärv, and G. Hagen, “Emulating _ab initio_ computations of infinite nucleonic matter,” (2022), in preparation.EOS
  * Shavitt and Bartlett (2009) I. Shavitt and R. J. Bartlett,  _Many-body Methods in Chemistry and Physics_ (Cambridge University Press, Cambridge UK, 2009).EOS
  * Smith and Gelfand (1992) A. F. M. Smith and A. E. Gelfand, “Bayesian statistics without tears: A sampling-resampling perspective,” Am. Stat. 46, 84–88 (1992).EOS
  * Jiang and Forssén (2022) Weiguang Jiang and Christian Forssén, “Bayesian probability updates using sampling/importance resampling: Applications in nuclear theory,” Front. in Phys. 10, 1058809 (2022), arXiv:2210.02507 [nucl-th] .EOS
  * Melendez  _et al._ (2022) J. A. Melendez, C. Drischler, R. J. Furnstahl, A. J. Garcia, and Xilin Zhang, “Model reduction methods for nuclear emulators,” J. Phys. G 49, 102001 (2022), arXiv:2203.05528 [nucl-th] .EOS
  * Frame  _et al._ (2018) D. Frame, R. He, I. Ipsen, D. Lee, D. Lee, and E. Rrapaj, “Eigenvector continuation with subspace learning,” Phys. Rev. Lett. 121, 032501 (2018).EOS
  * Sarkar and Lee (2021) A. Sarkar and D. Lee, “Convergence of Eigenvector Continuation,” Phys. Rev. Lett. 126, 032501 (2021), arXiv:2004.07651 [nucl-th] .EOS
  * Melendez  _et al._ (2021) J. A. Melendez, C. Drischler, A. J. Garcia, R. J. Furnstahl, and Xilin Zhang, “Fast & accurate emulation of two-body scattering observables without wave functions,” Phys. Lett. B 821, 136608 (2021), arXiv:2106.15608 [nucl-th] .EOS
  * Zhang and Furnstahl (2022) Xilin Zhang and R. J. Furnstahl, “Fast emulation of quantum three-body scattering,” Phys. Rev. C 105, 064004 (2022), arXiv:2110.04269 [nucl-th] .EOS
  * Demol  _et al._ (2020) P. Demol, T. Duguet, A. Ekström, M. Frosini, K. Hebeler, S. König, D. Lee, A. Schwenk, V. Somà, and A. Tichai, “Improved many-body expansions from eigenvector continuation,” Phys. Rev. C 101, 041302 (2020), arXiv:1911.12578 [nucl-th] .EOS
  * Djärv (2021) T. Djärv,  _JupiterNCSM: A Pantheon of Nuclear Physics_ , Ph.D. thesis, Chalmers (2021), ISBN: 978-91-7905-552-3.EOS
  * Coester (1958) F. Coester, “Bound states of a many-particle system,” Nucl. Phys. 7, 421 – 424 (1958).EOS
  * Coester and Kümmel (1960) F. Coester and H. Kümmel, “Short-range correlations in nuclear wave functions,” Nucl. Phys. 17, 477 – 485 (1960).EOS
  * Čížek (1966) Jiri Čížek, “On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell-Type Expansion Using Quantum-Field Theoretical Methods,” J. Chem. Phys. 45, 4256–4266 (1966).EOS
  * Kümmel  _et al._ (1978) H. Kümmel, K. H. Lührmann, and J. G. Zabolitzky, “Many-fermion theory in expS- (or coupled cluster) form,” Phys. Rep. 36, 1 – 63 (1978).EOS
  * Bishop (1991) R. F. Bishop, “An overview of coupled cluster theory and its applications in physics,” Theor. Chim. Acta 80, 95–148 (1991).EOS
  * Zeng  _et al._ (1998) Chen Zeng, D. J. J. Farnell, and R. F. Bishop, “An Efficient Implementation of High-Order Coupled-Cluster Techniques Applied to Quantum Magnets,” J. Stat. Phys. 90, 327–361 (1998).EOS
  * Mihaila and Heisenberg (2000) Bogdan Mihaila and Jochen H. Heisenberg, “Microscopic Calculation of the Inclusive Electron Scattering Structure Function in ${}^{16}O$,” Phys. Rev. Lett. 84, 1403–1406 (2000).EOS
  * Dean and Hjorth-Jensen (2004) D. J. Dean and M. Hjorth-Jensen, “Coupled-cluster approach to nuclear physics,” Phys. Rev. C 69, 054320 (2004).EOS
  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, “Coupled-cluster theory in quantum chemistry,” Rev. Mod. Phys. 79, 291–352 (2007).EOS
  * Hagen  _et al._ (2010) G. Hagen, T. Papenbrock, D. J. Dean, and M. Hjorth-Jensen, “Ab initio coupled-cluster approach to nuclear structure with modern nucleon-nucleon interactions,” Phys. Rev. C 82, 034330 (2010).EOS
  * Binder  _et al._ (2013) Sven Binder, Piotr Piecuch, Angelo Calci, Joachim Langhammer, Petr Navrátil, and Robert Roth, “Extension of coupled-cluster theory with a noniterative treatment of connected triply excited clusters to three-body hamiltonians,” Phys. Rev. C 88, 054319 (2013).EOS
  * Hagen  _et al._ (2014b) G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, “Coupled-cluster computations of atomic nuclei,” Rep. Prog. Phys. 77, 096302 (2014b).EOS
  * Gandolfi  _et al._ (2009) S. Gandolfi, A. Yu. Illarionov, K. E. Schmidt, F. Pederiva, and S. Fantoni, “Quantum monte carlo calculation of the equation of state of neutron matter,” Phys. Rev. C 79, 054005 (2009).EOS
  * Teschl (2009) G. Teschl, “Mathematical methods in quantum mechanics,” Graduate Studies in Mathematics 99, 106 (2009).EOS
  * C. E. Rasmussen (2006) C. K. I. Williams C. E. Rasmussen,  _Gaussian processes for machine learning_ (The MIT Press, 2006).EOS
  * Ducoin  _et al._ (2011) Camille Ducoin, Jerome Margueron, Constanca Providencia, and Isaac Vidana, “Core-crust transition in neutron stars: predictivity of density developments,” Phys. Rev. C 83, 045810 (2011), arXiv:1102.1283 [nucl-th] .EOS
  * Santos  _et al._ (2014) B. M. Santos, M. Dutra, O. Lourenço, and A. Delfino, “Correlations between the nuclear matter symmetry energy, its slope, and curvature from a nonrelativistic solvable approach and beyond,” Phys. Rev. C 90, 035203 (2014), arXiv:1409.4986 [nucl-th] .EOS
  * Coraggio  _et al._ (2014) L. Coraggio, J. W. Holt, N. Itaco, R. Machleidt, L. E. Marcucci, and F. Sammarruca, “Nuclear-matter equation of state with consistent two- and three-body perturbative chiral interactions,” Phys. Rev. C 89, 044321 (2014), arXiv:1402.0965 [nucl-th] .EOS
  * (64) See Supplemental Material at [URL will be inserted by publisher] for numerical parameters of a Gaussian approximation to the PPD.
  * Lattimer and Lim (2013) James M. Lattimer and Yeunhwan Lim, “Constraining the symmetry parameters of the nuclear interaction,” Ap. J. 771, 51 (2013).EOS
  * Bender  _et al._ (2003) Michael Bender, Paul-Henri Heenen, and Paul-Gerhard Reinhard, “Self-consistent mean-field models for nuclear structure,” Rev. Mod. Phys. 75, 121–180 (2003).EOS
  * Shlomo  _et al._ (2006) S Shlomo, VM Kolomietz, and G Colo, “Deducing the nuclear-matter incompressibility coefficient from data on isoscalar compression modes,” Eur. Phys. J. A 30, 23–30 (2006).EOS
  * Note (1) The prior for $c_{1,2,3,4}$ is the multivariate Gaussian resulting from a Roy-Steiner analysis of $\pi N$ scattering data Siemens  _et al._ (2017).
  * Tilley  _et al._ (2002) D. R. Tilley, C. M. Cheves, J. L. Godwin, G. M. Hale, H. M. Hofmann, J. H. Kelley, C. G. Sheu, and H. R. Weller, “Energy levels of light nuclei A=5, A=6, A=7,” Nucl. Phys. A 708, 3–163 (2002).EOS
  * Miller  _et al._ (2022a) Sean B. S. Miller, Andreas Ekström, and Kai Hebeler, “Neutron-deuteron scattering cross sections with chiral NN interactions using wave-packet continuum discretization,” Phys. Rev. C 106, 024001 (2022a), arXiv:2201.09600 [nucl-th] .EOS
  * Miller  _et al._ (2022b) Sean B. S. Miller, Andreas Ekström, and Christian Forssén, “Posterior predictive distributions of neutron-deuteron cross sections,” (2022b), arXiv:2209.06501 [nucl-th] .EOS
  * Wesolowski  _et al._ (2021b) S. Wesolowski, I. Svensson, A. Ekström, C. Forssén, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, “Rigorous constraints on three-nucleon forces in chiral effective field theory from fast and accurate calculations of few-body observables,” Phys. Rev. C 104, 064001 (2021b), arXiv:2104.04441 [nucl-th] .EOS
  * Siemens  _et al._ (2017) D. Siemens, J. Ruiz de Elvira, E. Epelbaum, M. Hoferichter, H. Krebs, B. Kubis, and U. G. Meißner, “Reconciling threshold and subthreshold expansions for pion–nucleon scattering,” Phys. Lett. B 770, 27–34 (2017), arXiv:1610.08978 [nucl-th] .EOS


# Emergent dynamics due to chemo-hydrodynamic self-interactions in active
polymers

Manoj Kumar<EMAIL_ADDRESS>Simons Centre for the Study of Living
Machines, National Centre for Biological Sciences, Tata Institute of
Fundamental Research, Bangalore, India  Aniruddh Murali  Simons Centre for the
Study of Living Machines, National Centre for Biological Sciences, Tata
Institute of Fundamental Research, Bangalore, India  Arvin Gopal Subramaniam
Department of Physics, Indian Institute of Technology, Chennai, India  Rajesh
Singh<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology, Chennai, India  Shashi Thutupalli<EMAIL_ADDRESS>Simons
Centre for the Study of Living Machines, National Centre for Biological
Sciences, Tata Institute of Fundamental Research, Bangalore, India
International Centre for Theoretical Sciences, Tata Institute of Fundamental
Research, Bangalore, India

###### Abstract

We create freely-jointed active polymers using self-propelled droplets as
monomeric units, and show that self-shaping chemo-hydrodynamic interactions
within the polymer result in novel dissipative structures and steady-states.
Our experiments reveal that the interactions between the monomeric droplets
give rise to ballistic propulsion of the active polymers and is associated
with rigidity and stereotypical shapes of the polymers. These traits,
quantified by the curvatures and speeds of the active polymers, vary
systematically with the chain length. Using simulations of a minimal model, we
establish that the emergent propulsion and rigidity are a generic consequence
of quasi two-dimensional confinement and auto-repulsive chemical interactions
between the freely jointed active droplets. Finally, we tune the interplay
between the chemical and hydrodynamic fields to experimentally demonstrate
oscillatory dynamics of the rigid polymer propulsion. Altogether, our work
highlights the possible first steps towards synthetic self-morphic active
matter.

Understanding and tuning the orchestration between mechanics and other guiding
fields, especially self-generated ones such as chemistry, is a fundamental
challenge in the field of active matter Ramaswamy (2010); Marchetti  _et al._
(2013); Cates and Tailleur (2015); Bechinger  _et al._ (2016); Saha  _et al._
(2014), both in the context of living systems Howard  _et al._ (2001); Berg
(1975); Gross  _et al._ (2017) and synthetic non-living emulations Howse  _et
al._ (2007); Paxton  _et al._ (2004); Thutupalli  _et al._ (2011); Hokmabad
_et al._ (2022); Meredith  _et al._ (2020); Ren  _et al._ (2020). While there
has been an increasing focus of experimental work on the emergent consequences
of chemo-mechanical self-coupling in synthetic active matter systems comprised
of point-like unit active particles Aubret  _et al._ (2018); Hokmabad  _et
al._ (2022); Meredith  _et al._ (2020), the experimental realisation and study
of such dynamics in extended active objects such as molecules Löwen (2018),
polymers Jayaraman  _et al._ (2012); Laskar and Adhikari (2015); Winkler  _et
al._ (2017); Prathyusha  _et al._ (2022) and sheets Manna  _et al._ (2022)
remains elusive.

Several approaches have attempted the construction of active assemblies –
especially linear chains – using monomer units such as colloids, to serve as
idealised experimental realizations of multi-body “polymeric” systems. It must
be noted that the assemblies have mostly been achieved and kept out of
equilibrium using external electric fields, magnetic fields, or light energy
Zhang and Granick (2016); Vutukuri  _et al._ (2017); Nishiguchi  _et al._
(2018); Snezhko and Aranson (2011); Dreyfus  _et al._ (2005); Stenhammar  _et
al._ (2016). While these systems are out-of-equilibrium, the dynamics is
predominantly due to external fields, in contrast to active and living systems
in which drive is at the scale of the individual units. This difference in the
external versus internal driving is crucial – since the control fields are
imposed externally, the interactions between the monomers are not self-
responsive, leading to a fundamental difference in the emergent dynamics
Ramaswamy (2010); Marchetti  _et al._ (2013); Gross  _et al._ (2017).
Chemically linked catalyst-coated colloids, that are actuated locally, have
also been used to construct linear polymers Biswas  _et al._ (2017, 2021) —
however, the “monomeric” units are neither self-propelled nor are
orientationally free. Emulsion droplets, due to their internal fluidity offer
a promising monomeric unit — in the past, passive linear assemblies that are
flexible or freely jointed have been made using emulsion droplets using sticky
DNA linkers Zhang  _et al._ (2017); McMullen  _et al._ (2018, 2022); Zhang
_et al._ (2018).

In this study, we construct polymeric chains of freely jointed self-propelling
active droplets. The propulsion mechanism of the monomers, _i.e._ droplets,
creates external chemical and hydrodynamic fields, causing self-shaping
interactions within the polymer. We measure these chemical and hydrodynamic
fields and then quantify, as a function of the polymer lengths, the emergent
rigidity, shapes and ballistic self-propulsion of the active polymers in quasi
two-dimensional confinement. We show that a model which includes only chemical
interactions between the monomers captures all these activity-driven
conformational features quantitatively, thereby establishing the minimal
requirements for the emergent self-organization: (i) the monomer droplets of
the chain propel due to (gradients of) their self-generated chemical field and
(ii) the structures are confined to move in quasi two-dimensions. Finally, by
tuning the coupling between the chemical and hydrodynamic fields, we
demonstrate oscillatory gaits of the polymers. Altogether, we envision these
as first steps towards a kind of synthetic active matter capable of emergent
self-morphic dynamics.

###  .1 Freely jointed active polymers of self-propelling emulsion droplets

The self-propelled droplets that we use as monomers are comprised of oil,
slowly dissolving into an external supramicellar aqueous solution of ionic
surfactants Peddireddy  _et al._ (2012); Peddireddy (2014). The dissolution of
the droplets results in the spontaneous development of self-sustaining
gradients of surfactant coverage around the droplets. These gradients give
rise to Marangoni stresses causing the droplets to propel Herminghaus  _et
al._ (2014) (Fig. 1A). As such, these self-propelled droplets may be viewed as
a physico-chemical realization of the so-called squirmer model for
microswimmers, a sphere with a prescribed surface slip velocity that exchanges
momentum with the fluid in which it is immersed Blake (1971); Ehlers  _et al._
(1996); Ishikawa  _et al._ (2006). The qualitative and quantitative nature of
this slip velocity governs the near- and far-field hydrodynamic flow
perturbation around the squirmer and eventually its self-propulsion.

The hydrodynamic flow fields around the oil droplet microswimmer (Fig. 1B,
Supplementary Video SV1) depend in tunable ways on geometric and, as we also
show later, chemical conditions Peddireddy (2014); Hokmabad  _et al._ (2021);
Dwivedi  _et al._ (2021); Ramesh  _et al._ (2022). In addition to the
hydrodynamic flows, these microswimmers leave a trail of chemical fields
Hokmabad  _et al._ (2022) comprised of oil-filled surfactant micelles formed
by the transfer of oil molecules into empty surfactant micelles, which can be
visualised using an oil-soluble fluorescent dye (Fig. 1C, Supplementary Video
SV1). The trail persists because the oil-filled micelles take longer to
diffuse than the surrounding oil-free surfactant micelles and causes a
remodeling of the environment around the droplets, leading to time delayed
negative chemotatic self-interactions Hokmabad  _et al._ (2022).

To create tethered assemblages, we use biotin-streptavidin chemistry to form
freely jointed chains of the active droplets (Fig. 1D, Materials and Methods).
Briefly, we use a hybrid (surfactant-lipid) monolayer of surfactant and
biotinylated lipids that coat the droplets spontaneously (SI Fig. S1A, B).
Using these biotinylated lipids in conjunction with the well known specific
interaction of biotin with streptavidin, we “polymerise” the droplets to form
linear chains via controlled assembly and incubation (details in Materials and
Methods and Supplementary Information, SI Fig. S1B). This assembly process
robustly results in linear assemblies of 5CB oil emulsion droplets such as
dimers, trimers, and so on up to decamers (Fig. 1E). Our protocol rarely
yields chains longer than $N=10$ without branching, and longest linear chains
we obtain consist of $N=13$ droplets (Fig. 1E and SI Fig. S1D). The droplets
comprising the polymers are freely-jointed, as seen from the distribution of
the internal bond angles (Fig. 1F), only restricted by the steric interaction
between the droplets. Furthermore, since these are liquid droplets, there are
also convective flows that are generated inside these self-propelled oil
droplets Thutupalli  _et al._ (2018) which can cause dynamic reorientation of
the propulsion direction (this is in strong contrast with asymmetric colloids
such as Janus colloids, whose orientation is fixed by design). We activate the
chains by immersing them in 25 wt% SDS surfactant solution. This causes the
monomers within the droplet to propel (Fig. 1G) and thereby the entire polymer
is set into motion – this propulsion is fully three-dimensional and for ease
of visualisation we only focused on confined settings (using Hele-Shaw like
geometries) in this paper (SI Fig. S1E). Furthermore, since these are liquid
droplets, there are also convective flows that are generated inside these
self-propelled oil droplets Thutupalli  _et al._ (2018) which can cause
dynamic reorientation of the propulsion direction; these fully flexible active
polymer dynamics are evident in the partially three-dimensional random motion
and orientation of the polymer chains (Fig. 1G and Supplementary Video SV2).

Figure 1: Self-propelled droplets are used as monomers to assemble freely-
jointed active polymers. (A) A monomeric liquid crystal droplet in a
surfactant micellar solution, swims due to a spontaneously generated
surfactant gradient at its interface. (B) The propulsion of the nematic liquid
crystal 5CB droplets results in characteristic hydrodynamic flow fields (left
half: experimental measurements and right half: simulations of a squirmer
model) and (C) chemical fields visualised using oil-soluble fluorescent (Nile
red) dye mixed with 5CB. Inset: the symmetry breaking associated with the
self-propulsion and the defect structure within the nematic droplet is
apparent in the bright field and cross-polarised images of the self-propelling
droplets. (D) Schema for assembling linear polymers of the active droplets
using biotin–streptavidin chemistry. (E) Linear, chemically linked, inactive
N-meric chains of 5CB emulsion droplets, comprised of increasing numbers of
monomer units. (F) The bond angles of the polymer exhibit a uniform
distribution i.e. the polymers are freely-jointed. (G) When activated, i.e.
immersed into surfactant solutions of high enough concentration, the polymers
exhibit active motion. The image shows three time points of the polymer,
overlaid with the center of mass trajectory, moving in a weakly confined
setting. (All scale bars represent $50~{}\mu m$).

###  .2 Emergent rigidity and ballistic propulsion of the polymers in quasi
two-dimensional confinement

When the active polymers are strongly confined to a quasi two-dimensional
setting _i.e._ a Hele-Shaw cell, the chains no longer exhibit flexible motion
but rather become rigid and adopt stereotypic C-like shapes (Fig. 2A, B,
Supplementary Videos SV3). We quantify the curvature of the steady state
chain, via the Monge representation Helfrich (1973), which reveals a
systematic dependence, i.e. a straightening of the polymer, with increase in
the polymer length, suggesting stiffening of the polymers (Fig 2C). To test
this further, we compute the “positional-” and ”orientational rigidity” of the
chains. The positional rigidity is quantified by the distribution of the
positions of the monomers from center of mass of the chain along both parallel
and perpendicular directions (Fig. 2C inset and Fig. 2D, SI Fig. S1F). The
orientational rigidity is computer by computing the correlation between the
orientations of the monomers with respect direction orthogonal to the long
axis of the polymer (Fig. 2C inset and Fig. 2E). In weak confinements
($h/2b=1.4$), the histograms for the positional (Fig. 2D) and orientational
(Fig. 2E) rigidities exhibit a broad distribution for droplet positions
indicating a freely-jointed or flexible nature of the chain. However, in
strong confinement $(h/2b=1)$, we see distinct peaks for the positions and
orientations of the monomers within the chains, clearly indicating an
emergence of rigidity of the C-shape configurations in strong confinement
(Fig. 2D–E). This emergent rigidity of the active chains is associated with
ballistic propulsion of the polymers in a direction orthogonal to the length
of the polymer. As with the curvature, the propulsion is sensitive to the
number of monomeric units comprising the polymer (Fig. 2F–G, Supplementary
Videos SV4) with longer polymers moving faster than shorter ones — it is
noteworthy that while the speeds of the polymers less than that of the
monomeric units, they asymptote to the instantaneous speed of a monomeric
droplet with increasing polymer length (Fig. 2H). Altogether, we find that in
confined settings, the active polymers exhibit emergent rigidity, straighten
with increasing polymer lengths and propel ballistically in a direction
orthogonal to the long axis of the polymers(Supplementary Videos SV4).

Figure 2: Emergent ballistic propulsion and rigidity of the polymer in strong
quasi two-dimensional confinement. A and B The active chains are flexible in
weak confinement ($\textit{h/2b}=1$) and become rigid in strong confinement
($\textit{h/2b}=1.4$). Panels A and B are the time-dependent configurations
(left to right) of the active chains in their weak and strong 2D confinements
(scale bar: $50~{}\mu m$). (C) Dimensionless mean curvature (obtained by
multiplying curvature with the droplet radius) of different polymer chain
lengths decreases with chain length ($N$) (Supplementary Experimental Section,
SI Fig. S2 and Fig. S3) . D and E are the distributions of the positional and
orientational rigidity of the active polymer chain measured in weak
confinement ($\textit{h/2b}=1$) and in strong confinement
($\textit{h/2b}=1.4$). The positional rigidity measured from the chain centre
of the mass in a direction parallel and orthogonal to the direction of motion
of the centre of mass. The orientational rigidity of the chain measured as
$e_{i}$*$n_{CM}$ where $e_{i}$ is the orientation of the unit vector $n_{i}$
and $n_{CM}$ is the unit vector of the center of mass. (F) Trajectories of
different polymer chain lengths are shown in different colors, follow the (G)
for the color representations (Supplementary Experimental Section, SI Fig. S2
and Fig. S3). (G) Mean squared displacement (MSD) profiles of experimental
trajectories (scale bar: 1000 $\mu m$) for different length chains ($N=1$ to
$N=7$), represented in dimensionless form by dividing the MSD with $b^{2}$ and
time with $1~{}s$ where $b$ is the radius of the droplet (Supplementary Video
SV4). Dashed line indicates the ballistic $\sim t^{2}$ limit. (H) Normalised
speed of different polymer chain lengths ($N$) which increases with the chain
length. The speed of a polymer chain is normalised with respect to the monomer
speed, $v_{s}$ (Supplementary Experimental Section, SI Fig. S3).

###  .3 Self-shaping chemo-hydrodynamic interactions within the polymer

Figure 3: Self-shaping hydrodynamic and chemical fields of active polymer
chains. (A) Experimentally measured hydrodynamic flow fields (Supplementary
Section PIV, FlowTrace), SI Fig. S(4C, D, E), from left to right: dimer
($N=2$), trimer ($N=3$), and octamer ($N=8$). (B) Experimentally measured
chemical fields (Supplementary Section CF, SI Fig. S4 A, B), from left to
right: dimer, trimer, and octamer (scale bar: $50~{}\mu m$). All error bars
indicate standard deviations.

To guide our understanding of these emergent traits, we start with
quantitatively mapping the self-shaping hydrodynamic and chemical interaction
fields within the polymer (exemplary hydrodynamic and chemical fields are
shown in Figs. 3A, B and associated Supplementary Videos SV5–SV10). The flow
fields were measured using standard microPIV techniques with fluorescent
tracer particles of diameter 0.5 $\mu m$ pre-mixed in the surrounding
continuous medium. On the other hand, the chemical fields were visualized
using an oil-soluble fluorescent dye (Supplementary Experimental Section).
When the oil droplets dissolve via micellar solubilization, micelles filled
with dye-doped oil are formed, enabling the visualisation of the field of
filled-oil-micelles (“chemical field”) around the monomer droplets in the
chain (Fig. 3B). The hydrodynamic flow fields around the chains are not
obtained from a simple superposition of the field around a monomer (Fig. 1B)
underscoring the non-linearity of the chemo-hydrodynamic coupling between the
monomers within the chains (Supplementary Information). Further, there is a
qualitative difference in the fields of a dimer from that of the other N-mers:
the dimer attains a metastable symmetric configuration during which it does
not propel while the steady-state hydrodynamic and chemical fields of all
other N-mers are asymmetric (see SI Fig. S6 A–H, Supplementary Videos
SV5–SV10). The symmetries of the chemical and hydrodynamic flow fields follow
each other and it is remarkable to note that steady state chemical fields are
also established in front of the chains in the direction of propulsion — the
polymers propel even through the negatively chemotactic fields that they self-
generate – this together with the force-balance resulting from the
orientations of the swimmers could account for the reduction in the speed of
the polymers. Despite the advection-diffusion coupled chemo-hydrodynamic
effects, we next ask if chemical interactions between the monomeric droplets
(in a freely jointed chain) alone can rationalize the experimental
observations. Indeed, such a rationalization might be anticipated a priori,
given that such approaches successfully describe the dynamics of individual
droplets Hokmabad  _et al._ (2022) and also that in confinement (two
dimensions), chemical fields decay slower than hydrodynamic fields Kanso and
Michelin (2019).

###  .4 Chemical interactions alone account for the emergent rigidity and
ballistic propulsion

We develop a minimal model for the active polymers in which the chemical field
around a monomer droplet is modeled by considering the $i$th monomer as a
point source of the chemical field centered at $\bm{R}_{i}$, which self-
propels with speed $v_{s}$ along the direction $\bm{e}_{i}$. The velocity and
direction of the monomer can change due to chemical interactions. The position
$\bm{R}_{i}$ and orientation, given by the unit vector $\bm{e}_{i}$, of the
$i$th monomer is updated using the following kinematic equations:

$\displaystyle\frac{d\bm{R}_{i}}{dt}={\bm{V}}_{i},\qquad\frac{d\bm{e}_{i}}{dt}={\bm{\Omega}}_{i}\times\bm{e}_{i}.$
(1)

Here, the translational velocity $\bm{V}_{i}$ and angular velocity
$\mathbf{\Omega}_{i}$ of the $i$th monomer are given as:

$\displaystyle{\bm{V}}_{i}=v_{s}\bm{e}_{i}+\chi_{t}\,\bm{\mathcal{J}}_{i}+\mu\bm{F}_{i},\qquad\bm{\Omega}_{i}=\chi_{r}\left(\bm{e}_{i}\times\bm{\mathcal{J}}_{i}\right).$
(2)

Figure 4: Simulations of a minimal model including only chemical interactions
capture the active polymer dynamics. (A) Snapshots of active chains from
simulations are shown in four columns, for chain sizes of $N=2$, $N=3$, $N=4$
and $N=8$, respectively (from left to right). Black arrows indicate the
orientation $\bm{e}_{i}$ of the particles. We show snapshots for the initial-
state ($t=0$), transient-state ($t=9s$), and steady states ($t=38s$ and
$93$s). The pseudo-color plot of the chemical field has also been overlaid on
the plots for configuration at $t=93$s. (B) Shows the saturating angles of the
chain in the steady state, with $\beta_{in}$ and $\beta_{out}$ retaining their
definitions as in Fig. 3. (C) shows the MSD for different chain lengths, along
with the ballistic limit of $\sim t^{2}$ as reference. (D) contains the speed
of polymers in steady state, normalised by the speed of a single particle
$v_{s}$ in simulations, as a function of the number of monomers in the polymer
($N$). (E) shows the dimensionless curvature of the active polymer in
simulations as a function of the $N$. For Figures F, H, and I, we display the
positional and angular dynamics for the $N=8$ chain. (F) is the distribution
of the angle between monomer inside the chain ($\beta_{in}$) and on the edge
of the chain ($\beta_{out}$) for the trajectory shown in panel (A). (G) is a
schematic of distance from the centre of mass in directions parallel
$r^{\parallel}-r_{CM}^{\parallel}$ and perpendicular
$r^{\perp}-r_{CM}^{\perp}$ to the motion of centre of mass, with the $r_{CM}$
coordinates denoted by dotted red lines. The distribution of these distances
are shown in (F) and (I) respectively for the dynamics shown in (A).

In the above equation, $\mu$ is the mobility of a monomer, while the chemical
interactions between the monomers are contained in the vector
$\bm{\mathcal{J}}_{i}=-\left(\bm{\nabla}c\right)_{\bm{r}=\bm{R}_{i}}$, where
$c$ is the concentration of filled micelles. The concentration field of the
filled micelles is obtained by considering each emulsion droplet as a point
source of the chemical field (explicit expressions of $c$ and
$\bm{\mathcal{J}}_{i}$ are given in the Supplementary Information). The
constants $\chi_{t}$ and $\chi_{r}$ take positive values. $\chi_{t}>0$ implies
a repulsive chemical interaction between the monomers while they are held
together by the attractive spring potential described below. The constant
$\chi_{r}>0$ implies that the monomers rotate away from each other in the
freely joined chain — such reorientation is consistent with the negative auto-
chemotactic behavior of the monomer droplets Hokmabad  _et al._ (2021), unlike
the chemical interactions of purely phoretic colloids. In the experiments, the
droplets are freely joined using biotin-streptavidin chemistry, which is
captured in the model by an attractive spring potential. The resulting spring
force on the $i$th monomer is given as:
$\bm{F}_{i}=-\left(\bm{\nabla}U\right)_{\bm{r}=\bm{R}_{i}}$, where
$U=\sum_{i=1}^{N-1}U^{C}(\mathbf{R}_{i},\mathbf{R}_{i+1})$ and
$U^{C}=k\left(r_{ij}-r_{0}\right)^{2}$ is spring potential of stiffness $k$
and natural length $r_{0}$ which holds the chain together. Here,
$r_{ij}=|\bm{R}_{i}-\bm{R}_{j}|$ and the spring’s natural length is equal to
the diameter of the monomer used in the experiment, $r_{0}=2b$. In our
experiment, the typical active force [$\mathcal{O}(6\pi\eta bv_{s})\sim
10^{-11}N$], with $\eta$ being the viscosity of the solvent is dominant to the
typical Brownian forces [$\mathcal{O}(k_{B}T/b)\sim 10^{-16}N$], and is thus,
ignored in our model. We simulate the model by integrating the above equations
numerically (simulation details and a complete list of parameters are given in
the Supplementary Information).

Snapshots from simulations of the above model are shown in Fig. 4A, where the
polymer structures are confined to move in two dimensions. This setting
corresponds to the strong confinement of our experiments. We find from our
simulations that the active polymer chain exhibits an emergent rigidity and
propels ballistically in C-configurations in the steady state, as is the case
with experiments as shown in Fig. 2 (Supplementary Video SV11). We also
measure the saturating angles between the monomers in the steady state; see
Fig. 4B. The MSD (mean-squared displacement) is also shown in Fig. 4C scales
as $\varpropto t^{2}$, which is in agreement with the experimentally measured
MSD (see Fig. 2G). It should be noted that a dimer stops propelling in the
steady state once the monomers comprising it point away from each other (i.e.
$\bm{e}_{1}\cdot\bm{e}_{2}=-1$). We find that the speed of the chain increases
with the chain length $N$ (which is the number of monomers in the chain), as
shown in Fig. 4D. Finally, we also show that the curvature decreases with the
chain length $N$; see Fig. 4E (see SI for details on computation of the
curvature). Thus, we show that the experimental observations in Fig. 2 are
reproduced from simulations (see Fig. 4) of our minimal model. We can now
focus on a particular chain, namely the $N=8$ case (i.e the dynamics shown in
Fig. 3A), right), and further quantify the positional and angular rigidity of
our system across its dynamical evolution; these are shown in Figures 4F,G, H
and I. Fig. 4F shows the angle between monomer inside the chain ($\beta_{in}$)
and on the edge of the chain ($\beta_{out}$) become fixed as the polymer
reaches a steady state. Here, peaks of the histogram correspond to the steady
state $\beta$s as in Fig. 4B, whereas the tail of the distribution corresponds
to the initial transient $\beta$s in the simulation before the chain acquires
rigidity. Fig. 4G shows the schematic of distance from the centre of mass in
directions parallel $r^{\parallel}-r_{CM}^{\parallel}$ and perpendicular
$r^{\perp}-r_{CM}^{\perp}$ to the motion of centre of mass, whose
distributions are shown in Fig. 4H and Fig. 4I, respectively. These
distributions are for the entire dynamical evolution of the chain, see SI for
corresponding distributions in the steady-state - which correspond to the
peaks of these figures. Overall, we see that the positions and hence internal
angles of the chain acquire fixed values as the system reaches its steady
state.

It is to be noted that the chemical fields generated in our model are not an
exact match to those as seen in the experiments – for instance, the dipolar
field for $N=2$ is not reproduced (compare Fig. 4A, top row with Fig. 3B).
Additionally, the dependence of the speed and curvature of the active chain on
the number of monomers is only a semi-quantitative match in that the
qualitative behaviour is reproduced but the exact numbers differ (compare the
y-axis of Fig. 2H with Fig. 4D and Fig. 2C with Fig. 4E). These differences
arise from (presumably) neglecting hydrodynamic flow fields in the
simulations, which advect the chemicals. Nevertheless, our model captures the
essential features of the emergent dynamical traits seen in the experiments,
indicating that chemical interactions between the monomers are responsible for
the emergent rigidity of the active polymers with hydrodynamics potentially
giving rise only to higher order quantitative corrections.

### Metastable configurations exhibiting spontaneous polymer rotations

To emphasise the versatility of the model, we now show that the chemical
interactions can also capture the dynamics of other metastable shapes of the
active polymers. As discussed earlier, the chains adopt a stable C-shape self-
propelling configuration. These shapes can be destabilised in two ways: (i) by
increasing the height of confinement such that the droplets forming the chain
can sample the full three dimensional space (e.g. Fig. 1F, Fig. 2A and
Supplementary Video SV2); and (ii) by direct-collision or interactions with
surrounding assemblies (e.g. Supplementary Video SV12) . In the experiments,
we often see such destabilisation (in quasi-2D settings due to the
interactions of the chain with surrounding assemblies) and this causes the the
chains to transition from the C-configurations to a metastable S-shaped
symmetric configuration (Fig. 5 along with Supplementary Videos SV13 and
SV14).

Figure 5: Persistent rotation in a metastable S-shaped configuration and its
transition to a stable C-shaped configuration. Longer active polymers
transiently adopt S-shaped configurations with associated characteristic (A)
hydrodynamic flow fields and (B) symmetry broken chiral chemical fields
(Supplementary Videos SV13, SV14). (C) In this metastable configuration the
polymers undergo persistent rotation shown by the snapshots and transition
spontaneously into the stable C-shaped configuration. Top: The top panel (left
to right) shows time snapshots of cross-polarised microscopy images of active
polymer chain configurations. The instantaneous propulsion directions of
droplets in the chain can be visualised with the help of a nematic director
pointing in the propulsion direction. The initial time ($t=0s$ and $40s$)
frames show that active chain rotation in transient S-configuration (indicated
with curved and red color arrows). At $t\sim 168s$, the active chain
completely transitions to a stable C-configuration and the arrow (white)
pointing in the chain translational direction (self-propulsion). Note that the
transition time (from S to C-configuration) can be different for different
chain lengths. Bottom: The above phenomena are reproduced in our simulations.
We see the rotation (curved red arrows) in the S-configuration ($t=0s$ and
$t=40s$), and a purely linearly translating (long black arrow) rigid chain a
$t=168s$.

In the S-configuration, the active chains exhibit a broken rotation symmetry
in the hydrodynamic and chemical fields and undergo spontaneous rotational
motion (Fig. 5A, B). In these configurations, the chains cannot escape long
term interaction with the self-generated chemical fields. As a consequence,
the S-configuration is only transiently stable, switching back to the
C-configuration via a series of monomer reorganizations (Fig. 5C, top panel,
Supplementary Videos SV15 and SV16). Simulations of our model, match this
dynamical behaviour (both the polymer rotation and spontaneous transition to
the C-shape) and show that, indeed, the dynamics of the S-configuration is due
to interactions with the self-generated chemical field (Fig. 5C, bottom
panel).

Our results so far suggest that the emergent dynamics of the flexible active
polymers — rigidity, stereoscopic shapes and propulsion — can be qualitatively
(and to some extent quantitatively) captured by considering the effects of
confinement and chemical interactions between the droplets. However, the
effects of the coupling between the chemical and hydrodynamic fields are
apparent in our system due to their mutual shaping of each other. Such
coupling can lead to feedback and time-delay effects which in some scenarios
can give rise to time-dependent steady-states such as oscillations beyond what
we have considered so far. Hints of such behavior can already be seen from the
dynamics of monomers and dimers (Fig. 2F). We now experimentally explore such
regimes of chemo-hydrodynamic coupling via tuning of the monomer propulsion.

### Tuning the chemo-hydrodynamic fields of the monomer self-propelled droplet

Figure 6: Tuning the chemo-hydrodynamic fields of a monomer. The schematics
showing swimming monomer droplets without additional oil-filled micelles
($\Phi=0$) (A) and with additional oil-filled micelles to the swimming droplet
environment ($\Phi=0.7$) (B). (C) Shows the speed of the monomer droplet
decreasing with increasing concentration of oil-filled micelles. (D)
Experimental trajectories of monomer droplets at a different concentration of
oil-filled micelles added to its environment. (E) The contractile
(quadrupolar) nature of the steady-state hydrodynamic field of monomer in a
chemically tuned environment with oil-filled micellar solution ($\Phi=0.7$)
and (F) is the corresponding modified steady-state chemical field. (G) Shows
the normalised chemical field intensity behind the swimmer, along the line
normal to the direction of motion of swimming monomer droplets.

.

The propulsion of the monomer can be tuned by modulating the transport of
fresh surfactant micelles to the droplet interface Dwivedi  _et al._ (2021);
Ramesh  _et al._ (2022, 2023). Specifically, such modulation affects the
distribution of the surfactant gradient, and thereby the slip velocity, on the
surface of the droplet — altogether, this qualitatively and quantitatively
changes the chemical and hydrodynamic flow fields around a self-propelling
monomer. We achieved such tuning by controlling the ratios between fresh
(“unfilled”) and “filled” surfactant micelles in the aqueous medium in which
the swimmers are embedded (Fig. 6A). The filled micelles are prepared by pre-
dissolving 5CB oil droplets in a 25 wt% SDS solution to a saturation. We then
adjust the ratio $\Phi$, the total fraction of the filled micelles in the
solution, using which we obtain a quantitative tuning of the monomer speed
(Fig. 6C) and its dynamics (Fig. 6D). It can be seen that with increasing
$\Phi$, the speed of the monomer reduces sharply. This reduction in speed is
associated with qualitative changes in the flow and chemical fields around the
monomer (Fig. 6E–F and SI Fig. S9). Since our focus is to demonstrate non-
trivial effects due to the coupling between the chemical and hydrodynamic
fields, we focus on the case with $\Phi=0.7$ which results in a flow field
that is reminiscent of a “puller” type of squirmer (Fig. 6E). Our rationale is
as follows — due to the reduced speed of the droplet in this configuration,
the associated chemical field is able to effectively diffuse further, and form
a wider trail behind the droplet (Fig. 6F–G). In addition, due to the nature
of the flow field, the filled micelles are transported to the front of the
droplet advectively thus changing time-scales associated with the self-
interaction of the monomers with their own chemical trails.

### Time varying chemo-hydrodynamic fields and the emergence of polymer
oscillations

We now turn our attention to the dynamics of the chains in this condition i.e.
$\Phi=0.7$. There are now striking differences in the chemical and
hydrodynamic flow-fields around the polymers (Fig. 7, in comparison with Fig.
3). Specifically, clear asymmetries in the fields appear in a time-snapshot
and these represent time-dependent dynamics (Supplementary Video SV17–SV18) –
these asymmetries are also reflected in the orientations of the monomers
within the chains. Generally, these represent complex dynamics of active
polymers and here, we specifically focus only on time-periodic motion
(oscillations). For this, we illustratively focus on the case of trimers —
oscillatory gaits of the trimers can readily be seen from the trajectories
(Fig. 8A) Supplementary Videos SV19). By measuring the time-varying chemical
fields around the trimer, we rationalise that these oscillations are caused
due to the repulsive autochemotactic nature of the droplets. It can be seen
that the interplay between the flow and chemical fields leads to an asymmetric
build-up of the chemical field in front of the chain (Fig. 8A, lower row) due
to which the monomer self-propulsion direction undergoes a re-orientation. The
periodicity of these oscillations can be seen clearly from the quantification
of the intensity and phase angle measurements $\phi$, which periodically
change with time on either side of the trimer (Fig. 8B) — the sharp
transitions are suggestive of threshold concentrations of the chemical field
that cause sharp re-orientations of the monomers.

Figure 7: Self-shaping hydrodynamic and chemical fields of active polymer
chains in a chemically tuned environment. (A) Experimentally measured
hydrodynamic flow fields (Supplementary Section PIV, FlowTrace), SI Fig. S4(C,
D, E), from left to right: dimer ($N=2$), trimer ($N=3$), and nonamer ($N=9$).
(B) Experimentally measured chemical fields (Supplementary Section CF, SI Fig.
S4A, B), from left to right: dimer, trimer, and nonamer (scale bar: $50~{}\mu
m$).

.

## Discussion

Using chemo-hydrodynamically active monomer droplets, we constructed freely
jointed polymers; not only are the joints between the monomers fully flexible
but the self-propulsion direction of each monomer can also freely evolve. The
self-interactions within the polymer give rise to emergent self-organization,
particularly self-propulsion — stable translation and metastable persistent
rotations — of the polymers in rigid configurations. We quantitatively mapped
these interactions and using a simple model based on chemical interactions
alone, identified minimal, generic conditions in which such novel self-
organization can arise. It must be noted that the active chains in our
experiments self-propel in a direction normal to their body axis, which is the
direction of maximum drag for a slender body; this is in contrast to natural
objects such as microbes and polymers which typically propel along their major
body axis _i.e._ the direction of minimum drag.

The “dry” chemical active matter model we have considered here is in line with
similar recent successful descriptions for the scattering dynamics of monomer
droplets due to their repulsive auto-chemotaxis Hokmabad  _et al._ (2022);
Kranz  _et al._ (2016). However, the effects of the hydrodynamic coupling are
apparent in our system due to their shaping of the chemical fields and vice-
versa. While such coupling does not qualitatively effect the emergent rigidity
and self-propulsion that we have reported here, feedback and time-delay
effects from such coupling can give rise to steady-states such as
oscillations. We have experimentally explored that such states may indeed be
achieved via tuning of the monomer propulsion Ramesh  _et al._ (2022); Dwivedi
_et al._ (2021) — developing the full theoretical framework for such dynamics
involving the coupling of the hydrodynamics with the chemical fields remains a
future challenge.

Figure 8: Oscillatory dynamics of active chain in a chemically tuned
environment. (A) Oscillatory trajectories of a trimer and the time -dependent
asymmetry in the corresponding chemical field in a strong two-dimensional
confinement and chemically tuned environment. (B) Time-dependent periodic
changes are observed in chemical field intensities measured along the droplet
orientations. The periodic changes in the droplet orientations are measured
for the droplet at the edge in reference to the orientations to the central
droplet.

.

The emergent dynamics we have reported are in quasi two-dimensional settings.
As we have seen, relaxation of this criterion results in a reduction of the
rigidity and affects the stability of the stereotypical polymer
configurations. For example, even when the height of the Hele-Shaw cell
$h>2b$, there is room for extra conformational degrees of freedom for the
freely jointed chain. The fully three-dimensional dynamics of active polymers
with chemical and hydrodynamic interactions should hold rich possibilities for
future exploration. Further, while we have only explored short linear
assemblies here, the extension to longer polymers and higher dimensional
assemblies such as sheets Manna  _et al._ (2022) with more tunable bonds
(“multiflavored assemblies”) Zhang  _et al._ (2018), complex shapes and
thereby a richer dynamical repertoire is possible. Such multicomponent
monomers with attractive, repulsive and even non-reciprocal interactions
Meredith  _et al._ (2020) may be used to create functional self-morphing
assemblies. Finally, emulsion droplets, as such, offer many possible
modifications in both controlling the internal chemistry of the droplets to
tune their activity and in controlling the flexibility of the assemblies
because of the mobile surface of the droplets in comparison to colloidal
systems. Therefore, these systems could also be step in the direction of
designing the multicomponent and multistimuli active delivery systems Downs
_et al._ (2020).

## Materials and Methods

###  .5 Materials

We used 5CB (4-cyno-4’-pentylbiphenyl) liquid crystal procured from Frinton
Laboratories, Inc. and SDS (sodium dodecyl sulphate) surfactant procured from
Sigma Aldrich. DOPC (1,2-dioleoyl-sn-glcero-3-phosphocoline), Liss-Rhod-PE
(1,2-dioleoyl-sn-glycero-phosphoethanolamine-N-(lissamine B sulfonyl)
(ammonium salt) and biotinyl-cap-PE (1,2-dioleoyl-sn-glycero-
phosphoethanolamine-N-(cap biotinyl) (sodium salt) were purchased from Avanti
Polar. Alexa Fluor®488 conjugated streptavidin was procured from Sigma-
Aldrich. FluoSpheres™ carboxylate-modified Microspheres, yellow-green
fluorescent (505/515 $nm$) and Nile red dye were purchased from Invitrogen (by
ThermoFisher Scientific). All chemicals were used as received.

###  .6 Preparation of (Lipid (DOPC) - surfactant (SDS)) micellar solution

All lipid stocks were prepared in chloroform: DOPC at 25 $mg/mL$, Biotinyl-
cap-PE at 10 $\mu g/mL$, and Liss-Rhod-PE at 10 $\mu g/mL$. We aliquoted 25
$\mu L$ of DOPC from the stock, 2 $\mu L$ of Biot-Cap-PE, and 10 $\mu L$ of
Liss-Rhod-PE in a 2 $mL$ effendorf tube to prepare the aqueous phase. The
lipids were vacuum dried overnight and then nitrogen dried before being mixed
with the aqueous phase. To the dried lipids containing eppendorf, we added 1
$mL$ of 0.125 (w/v)% SDS solution prepared with Milli-Q water. The eppendorf
was then left at room temperature (T = 25$\degree$C) for lipid hydration. A
probe sonicator (VC 750 (750W), stepped microtip diameter: 3 $mm$) with an
on/off pulse of 1.5 s/1.5 s and 32% of the maximum sonicator amplitude was
used to mix the components. The tube was placed in an ice bath during the
sonication steps (sonication for 30 s with a 1 min gap) to avoid overheating
of the sample.

###  .7 Droplet production

We used an oil-in-water (O/W) emulsion system to make the droplets. 5CB liquid
crystal droplets (oil phase) were stabilised using micellar solution
(surfactant and lipids) in the aqueous phase. When 5CB LC (the oil phase) is
injected through a microfluidic channel and encounters the aqueous phase (an
aqueous solution of lipids and surfactant, SI Fig. S1A), adsorption of free
lipid-surfactant results in stable 5CB droplets. The monodisperse droplets
generated by the microfluidic device were collected and stored in 0.25% SDS
solution, where they remained stable for months (SI Fig. S1C, left panel). The
presence of lipids in the monolayer of lipids and surfactant stabilising the
5CB emulsion droplets was confirmed using a red channel of the fluorescence
microscopy (SI Fig. S1C, middle panel), streptavidin functionalization of
biotinylated 5CB droplets was confirmed using a green channel of the
fluorescence microscopy (SI Fig. S1C, right panel) . Details of the
fabrication of microfluidic devices can be found in the Supplementary
Information, SI Fig. S1A, and the experimental section (DF).

###  .8 Preparation of droplet assemblies

We used a simple method to prepare the 5CB droplet assemblies: we split the
sample into two populations. Only biotinylated 5CB droplets were used in
population ”I” (SI Fig. S1B, a) and streptavidin functionalised biotinylated
5CB droplets in population ”II” (SI Fig. S1B, b, c, right panel). To remove
free biotinylated lipids and streptavidin from the bulk, both populations (I
and II) were washed 2-3 times with 0.25% SDS solution. The population ”I” is
then centrifuged at 6000 rpm for 30 s to settle down all the droplets at the
bottom of the tube and remove the solvent from the top. On top of the settled
fraction of the droplets, we gently added 200 $\mu L$ of fresh 0.25% SDS
solution. The tube was then filled with 5 $\mu L$ of a 50 $mM$ NaCl salt
solution prepared in 0.25% SDS solution, and it was centrifuged for 30 s (at
6000 rpm). We then gently added a 20 $\mu L$ volume of 5CB biotinylated
droplets functionalized with streptavidin from population ”II” to population
”I” at the bottom of the tube (SI Fig. S1B). The tube is then centrifuged at
12000 rpm for 35 min in an eppendorf centrifuge (Eppendorf centrifuge model
5418) at T = 15$\degree$C. The assemblies were studied inside a quasi two-
dimensional flow cell (SI Fig. S1E); the details of the flow cell fabrication
are discussed in the experimental section of the supporting information
(FlowCell).

###  .9 Imaging (Bright Field and Fluorescence Microscopy)

We used both multi-channel epifluorescence and brightfield microscopy to
visualize 5CB emulsion droplets and their linear assemblies. The assemblies
were imaged on an Olympus IX81 microscope with a 4x, 10x, and 20x UPLSAPO
objective, and images were captured with a Photometrics Prime camera. For the
fluorescence images, we used a CoolLED PE-4000 lamp. A 16-bit multichannel
image of 2048 x 2048 pixels was taken, consisting of a green channel
(excitation: 460 $nm$, dichroic: Sem-rock’s Quad Band), a red channel
(excitation: 550 $nm$, dichroic: Sem-rock’s Quad Band), and a bright field
channel. We recorded the images and videos using Olympus cellSens software. We
used a Leica bright field microscope (model M205FA) equipped with Leica-
DFC9000GT camera to record longer time videos of the linear assemblies. Videos
are recorded with a microscope zoom of 2, an exposure time of 0.02 s, and a 1
fps frame rate. We captured a 16-bit image of 2048 x 2048 pixels using a
Leica-DFC9000GT-VSC06748 camera.

## Data Availability Statement

No additional data was used besides the results described in the text.
Additional summary statistics of the data plotted may be available upon
reasonable request to the corresponding authors.

## Code Availability Statement

All code used to produce these results will be made available upon reasonable
request to the corresponding authors.

## References

  * Ramaswamy (2010) Sriram Ramaswamy, “The mechanics and statistics of active matter,” Annual Reviews of Condensed Matter Physics 1, 323–345 (2010).
  * Marchetti  _et al._ (2013) M Cristina Marchetti, Jean-François Joanny, Sriram Ramaswamy, Tanniemola B Liverpool, Jacques Prost, Madan Rao, and R Aditi Simha, “Hydrodynamics of soft active matter,” Reviews of Modern Physics 85, 1143 (2013).
  * Cates and Tailleur (2015) Michael E Cates and Julien Tailleur, “Motility-induced phase separation,” Annual Reviews of Condensed Matter Physics 6, 219–244 (2015).
  * Bechinger  _et al._ (2016) Clemens Bechinger, Roberto Di Leonardo, Hartmut Löwen, Charles Reichhardt, Giorgio Volpe, and Giovanni Volpe, “Active particles in complex and crowded environments,” Reviews of Modern Physics 88, 045006 (2016).
  * Saha  _et al._ (2014) Suropriya Saha, Ramin Golestanian, and Sriram Ramaswamy, “Clusters, asters, and collective oscillations in chemotactic colloids,” Physical Review E 89, 062316 (2014).
  * Howard  _et al._ (2001) Jonathon Howard _et al._ ,  _Mechanics of motor proteins and the cytoskeleton_, Vol. 743 (2001).
  * Berg (1975) Howard C Berg, “Chemotaxis in bacteria,” Annual review of biophysics and bioengineering 4, 119–136 (1975).
  * Gross  _et al._ (2017) Peter Gross, K Vijay Kumar, and Stephan W Grill, “How active mechanics and regulatory biochemistry combine to form patterns in development,” Annual review of biophysics 46, 337–356 (2017).
  * Howse  _et al._ (2007) Jonathan R Howse, Richard AL Jones, Anthony J Ryan, Tim Gough, Reza Vafabakhsh, and Ramin Golestanian, “Self-motile colloidal particles: from directed propulsion to random walk,” Physical review letters 99, 048102 (2007).
  * Paxton  _et al._ (2004) Walter F Paxton, Kevin C Kistler, Christine C Olmeda, Ayusman Sen, Sarah K St. Angelo, Yanyan Cao, Thomas E Mallouk, Paul E Lammert, and Vincent H Crespi, “Catalytic nanomotors: autonomous movement of striped nanorods,” Journal of the American Chemical Society 126, 13424–13431 (2004).
  * Thutupalli  _et al._ (2011) Shashi Thutupalli, Ralf Seemann, and Stephan Herminghaus, “Swarming behavior of simple model squirmers,” New Journal of Physics 13, 073021 (2011).
  * Hokmabad  _et al._ (2022) Babak Vajdi Hokmabad, Jaime Agudo-Canalejo, Suropriya Saha, Ramin Golestanian, and Corinna C Maass, “Chemotactic self-caging in active emulsions,” Proceedings of the National Academy of Sciences 119, e2122269119 (2022).
  * Meredith  _et al._ (2020) Caleb H Meredith, Pepijn G Moerman, Jan Groenewold, Yu-Jen Chiu, Willem K Kegel, Alfons van Blaaderen, and Lauren D Zarzar, “Predator–prey interactions between droplets driven by non-reciprocal oil exchange,” Nature Chemistry 12, 1136–1142 (2020).
  * Ren  _et al._ (2020) Lin Ren, Liyuan Wang, Qingyu Gao, Rui Teng, Ziyang Xu, Jing Wang, Changwei Pan, and Irving R Epstein, “Programmed locomotion of an active gel driven by spiral waves,” Angewandte Chemie 132, 7172–7178 (2020).
  * Aubret  _et al._ (2018) Antoine Aubret, Mena Youssef, Stefano Sacanna, and Jérémie Palacci, “Targeted assembly and synchronization of self-spinning microgears,” Nature Physics 14, 1114–1118 (2018).
  * Löwen (2018) Hartmut Löwen, “Active colloidal molecules,” Europhysics Letters 121, 58001 (2018).
  * Jayaraman  _et al._ (2012) Gayathri Jayaraman, Sanoop Ramachandran, Somdeb Ghose, Abhrajit Laskar, M Saad Bhamla, PB Sunil Kumar, and R Adhikari, “Autonomous motility of active filaments due to spontaneous flow-symmetry breaking,” Physical review letters 109, 158302 (2012).
  * Laskar and Adhikari (2015) Abhrajit Laskar and R Adhikari, “Brownian microhydrodynamics of active filaments,” Soft Matter 11, 9073–9085 (2015).
  * Winkler  _et al._ (2017) Roland G Winkler, Jens Elgeti, and Gerhard Gompper, “Active polymers—emergent conformational and dynamical properties: A brief review,” Journal of the Physical Society of Japan 86, 101014 (2017).
  * Prathyusha  _et al._ (2022) KR Prathyusha, Falko Ziebert, and Ramin Golestanian, “Emergent conformational properties of end-tailored transversely propelling polymers,” Soft Matter 18, 2928–2935 (2022).
  * Manna  _et al._ (2022) Raj Kumar Manna, Abhrajit Laskar, Oleg E Shklyaev, and Anna C Balazs, “Harnessing the power of chemically active sheets in solution,” Nature Reviews Physics 4, 125–137 (2022).
  * Zhang and Granick (2016) Jie Zhang and Steve Granick, “Natural selection in the colloid world: active chiral spirals,” Faraday Discussions 191, 35–46 (2016).
  * Vutukuri  _et al._ (2017) Hanumantha Rao Vutukuri, Bram Bet, René Van Roij, Marjolein Dijkstra, and Wilhelm TS Huck, “Rational design and dynamics of self-propelled colloidal bead chains: from rotators to flagella,” Scientific reports 7, 16758 (2017).
  * Nishiguchi  _et al._ (2018) Daiki Nishiguchi, Junichiro Iwasawa, Hong-Ren Jiang, and Masaki Sano, “Flagellar dynamics of chains of active janus particles fueled by an ac electric field,” New Journal of Physics 20, 015002 (2018).
  * Snezhko and Aranson (2011) Alexey Snezhko and Igor S Aranson, “Magnetic manipulation of self-assembled colloidal asters,” Nature Materials 10, 698–703 (2011).
  * Dreyfus  _et al._ (2005) Rémi Dreyfus, Jean Baudry, Marcus L Roper, Marc Fermigier, Howard A Stone, and Jérôme Bibette, “Microscopic artificial swimmers,” Nature 437, 862–865 (2005).
  * Stenhammar  _et al._ (2016) Joakim Stenhammar, Raphael Wittkowski, Davide Marenduzzo, and Michael E Cates, “Light-induced self-assembly of active rectification devices,” Science Advances 2, e1501850 (2016).
  * Biswas  _et al._ (2017) Bipul Biswas, Raj Kumar Manna, Abhrajit Laskar, PB Sunil Kumar, Ronojoy Adhikari, and Guruswamy Kumaraswamy, “Linking catalyst-coated isotropic colloids into “active” flexible chains enhances their diffusivity,” ACS Nano 11, 10025–10031 (2017).
  * Biswas  _et al._ (2021) Bipul Biswas, Debarshi Mitra, Fayis Kp, Suresh Bhat, Apratim Chatterji, and Guruswamy Kumaraswamy, “Rigidity dictates spontaneous helix formation of thermoresponsive colloidal chains in poor solvent,” ACS Nano 15, 19702–19711 (2021).
  * Zhang  _et al._ (2017) Yin Zhang, Angus McMullen, Lea-Laetitia Pontani, Xiaojin He, Ruojie Sha, Nadrian C Seeman, Jasna Brujic, and Paul M Chaikin, “Sequential self-assembly of dna functionalized droplets,” Nature Communications 8, 21 (2017).
  * McMullen  _et al._ (2018) Angus McMullen, Miranda Holmes-Cerfon, Francesco Sciortino, Alexander Y Grosberg, and Jasna Brujic, “Freely jointed polymers made of droplets,” Physical Review Letters 121, 138002 (2018).
  * McMullen  _et al._ (2022) Angus McMullen, Maitane Muñoz Basagoiti, Zorana Zeravcic, and Jasna Brujic, “Self-assembly of emulsion droplets through programmable folding,” Nature 610, 502–506 (2022).
  * Zhang  _et al._ (2018) Yin Zhang, Xiaojin He, Rebecca Zhuo, Ruojie Sha, Jasna Brujic, Nadrian C Seeman, and Paul M Chaikin, “Multivalent, multiflavored droplets by design,” Proceedings of the National Academy of Sciences 115, 9086–9091 (2018).
  * Peddireddy  _et al._ (2012) Karthik Peddireddy, Pramoda Kumar, Shashi Thutupalli, Stephan Herminghaus, and Christian Bahr, “Solubilization of thermotropic liquid crystal compounds in aqueous surfactant solutions,” Langmuir 28, 12426–12431 (2012).
  * Peddireddy (2014) Karthik Reddy Peddireddy, “Liquid crystals in aqueous ionic surfactant solutions: Interfacial instabilities & optical applications,” PhD Dissertation, Göttingen University (2014).
  * Herminghaus  _et al._ (2014) Stephan Herminghaus, Corinna C Maass, Carsten Krüger, Shashi Thutupalli, Lucas Goehring, and Christian Bahr, “Interfacial mechanisms in active emulsions,” Soft matter 10, 7008–7022 (2014).
  * Blake (1971) John R Blake, “A spherical envelope approach to ciliary propulsion,” Journal of Fluid Mechanics 46, 199–208 (1971).
  * Ehlers  _et al._ (1996) KM Ehlers, HC Berg, and R Montgomery, “Do synechococcus swim using traveling surface waves,” Proceedings of the National Academy of Sciences 93, 8340–8344 (1996).
  * Ishikawa  _et al._ (2006) Takuji Ishikawa, MP Simmonds, and Timothy J Pedley, “Hydrodynamic interaction of two swimming model micro-organisms,” Journal of Fluid Mechanics 568, 119–160 (2006).
  * Hokmabad  _et al._ (2021) Babak Vajdi Hokmabad, Ranabir Dey, Maziyar Jalaal, Devaditya Mohanty, Madina Almukambetova, Kyle A Baldwin, Detlef Lohse, and Corinna C Maass, “Emergence of bimodal motility in active droplets,” Physical Review X 11, 011043 (2021).
  * Dwivedi  _et al._ (2021) Prateek Dwivedi, Bishwa Ranjan Si, Dipin Pillai, and Rahul Mangal, “Solute induced jittery motion of self-propelled droplets,” Physics of Fluids 33, 022103 (2021).
  * Ramesh  _et al._ (2022) Prashanth Ramesh, Babak Vajdi Hokmabad, Myriam Rahalia, Arnold JTM Mathijssen, Dmitri O Pushkin, and Corinna C Maass, “Interfacial activity dynamics of confined active droplets,” arXiv preprint arXiv:2202.08630 (2022).
  * Thutupalli  _et al._ (2018) Shashi Thutupalli, Delphine Geyer, Rajesh Singh, Ronojoy Adhikari, and Howard A. Stone, “Flow-induced phase separation of active particles is controlled by boundary conditions,” Proceedings of the National Academy of Sciences 115, 5403–5408 (2018), https://www.pnas.org/content/115/21/5403.full.pdf .
  * Helfrich (1973) Wolfgang Helfrich, “Elastic properties of lipid bilayers: theory and possible experiments,” Zeitschrift für Naturforschung c 28, 693–703 (1973).
  * Kanso and Michelin (2019) Eva Kanso and Sébastien Michelin, “Phoretic and hydrodynamic interactions of weakly confined autophoretic particles,” The Journal of Chemical Physics 150, 044902 (2019).
  * Ramesh  _et al._ (2023) Prashanth Ramesh, Yibo Chen, Petra Räder, Svenja Morsbach, Maziyar Jalaal, and Corinna C Maass, “Arrested on heating: controlling the motility of active droplets by temperature,” arXiv preprint arXiv:2303.13442 (2023).
  * Kranz  _et al._ (2016) W Till Kranz, Anatolij Gelimson, Kun Zhao, Gerard CL Wong, and Ramin Golestanian, “Effective dynamics of microorganisms that interact with their own trail,” Physical Review Letters 117, 038101 (2016).
  * Downs  _et al._ (2020) Florence G Downs, David J Lunn, Michael J Booth, Joshua B Sauer, William J Ramsay, R George Klemperer, Craig J Hawker, and Hagan Bayley, “Multi-responsive hydrogel structures from patterned droplet networks,” Nature chemistry 12, 363–371 (2020).

## Acknowledgements

We acknowledge discussions with Sriram Ramaswamy, Ignacio Pagonabarraga, and
Ronojoy Adhikari; and thank Abhrajit Laskar for his help with the numerical
computation of the flow fields. We acknowledge support from the Department of
Atomic Energy (India), under project no. RTI4006, the Simons Foundation (Grant
No. 287975), the Human Frontier Science Program and the Max Planck Society
through a Max-Planck-Partner-Group. Work funded in part by the Indian
Institute of Technology, Madras, India through their seed and initiation
grants and Start-up Research Grant, SERB, India, to RS.

## Author Contributions

ST, RS, and MK conceived and designed the study. MK and AM performed the
experiments and data analysis with help from ST. RS and AGS designed the
theoretical model and carried out the simulations. MK, RS, AGS, and ST wrote
the paper.

## Competing Interests Statement

The authors declare no competing interests.