# Shakebot: A Low-cost, Open-source Robotic Shake Table for Earthquake
Research and Education

Zhiang Chen, Devin Keating, Yash Shethwala, Aravind Adhith Pandian
Saravanakumaran,
Ramon Arrowsmith, Albert Kottke, Christine Wittich, Jnaneshwar Das

###### Abstract

Shake tables provide a critical tool for simulating earthquake events and
testing the response of structures to seismic forces. However, existing shake
tables are either expensive or proprietary. This paper presents the design and
implementation of a low-cost, open-source shake table named Shakebot for
earthquake engineering research and education, built using Robot Operating
System (ROS) and robotic concepts. The Shakebot adapts affordable and high-
accuracy components from 3D printers, particularly a closed-loop stepper motor
for actuation and a toothed belt for transmission. The stepper motor enables
the bed to reach a maximum horizontal acceleration of 11.8 $m/s^{2}$ (1.2
$g$), and velocity of 0.5 $m/s$, with a 2 $kg$ specimen. The Shakebot is
equipped with an accelerometer and a high frame-rate camera for bed motion
estimation. The low cost and easy use make the Shakebot accessible to a wide
range of users, including students, educators, and researchers in low-resource
settings. An important application of the Shakebot is to examine the dynamics
of precariously balanced rocks (PBRs), which are negative indicators of
earthquakes in nature. Our earlier research built a virtual shake robot in
simulation for the PBR study. The Shakebot provides an approach to validate
the simulation through physical experiments. The ROS-based perception and
motion software facilitates the code transition from our virtual shake robot
to the physical Shakebot. The reuse of the control programs ensures that the
implemented ground motions are consistent for both the simulation and physical
experiments, which is critical to validate our simulation experiments.

##  I Introduction

For decades, researchers have used shake tables to simulate the effects of
earthquakes on structures, systems, and materials, which provides a critical
approach to earthquake research and education [1, 2]. Simulated earthquake
conditions can help civil engineers identify potential weaknesses in the
design and suggest ways to improve the structure’s resilience [3, 4]. In a
controlled laboratory setting, shake tables can also be used to study the
behavior of fragile geological features, which aids seismologists in finding
upper-bound constraints on ground motions [5]. By simulating earthquakes in
the classroom, students can better understand the forces involved and the
factors that influence the behavior of structures [6], inspiring the next
generation of engineers and seismologists.

(a) Shakebot CAD model

(b) Shakebot experimental setup

Figure 1: Shakebot CAD model and actual experimental setup. Stepper motor
actuates pedestal through toothed belt and pulley transmission. Control box
includes power supply, stepper motor driver, Raspberry Pi, emergency push
button, and touchscreen. For scale, the control box has a dimension of
390$\times$290$\times$160 $mm$.

An important application of shake tables is to examine the dynamics of
precariously balanced rocks (PBRs), which are free-standing rocks balanced on,
but not attached to sub-horizontal pedestals [7, 8, 9]. Because an earthquake
ground motion of sufficient amplitude and duration may topple PBRs, the
fragile configurations of PBRs provide upper-bound constraints on ground
motions, refining probabilistic seismic hazard analysis [10, 11]. The dynamics
of the PBR overturning process are non-linear [12, 5]. Purvance et al. [5]
used a shake table to study the PBR overturning dynamics. However,
difficulties can arise when using conventional shake tables to repeat
overturning experiments. For example, when analyzing PBR fragility, shake
experiments would involve accurate and repeated positioning of a large boulder
on the shake table to record the overturning responses given ground motions.
These boulders are often so heavy that the use of a crane or forklift is
needed, and repositioning them precisely is difficult (e.g., [13]),
complicating experiment repeatability.

Table I: Shake tables and their specifications. | Open source | Cost (USD) | Payload ($kg$) | Vel ($m/s$) | Acc ($g$) | Software | Programmable | Transmission
---|---|---|---|---|---|---|---|---
Quanser [14] | $\times$ | 30k | 7.5 | 0.4 | 2.5 | NA | ✓ | ball screw
CUSHAKE [15] | $\times$ | 45k | 3500 | 0.01 | 1 | DEPSIM | ✓ | threaded rod
Kelvin [16] | $\times$ | 1.3k | NA | NA | NA | NA | $\times$ | reciprocating arm
SARSAR [17] | ✓ | 1.8k | 200 (0.4$g$) | 0.35 | 1.5 | NA | ✓ | ball screw
Shao and Enyart [18] | $\times$ | NA | 228 | NA | 4 | NA | ✓ | hydraulic actuators
Kınay [19] | ✓ | NA | 80 | 0.25 | 2 | Scilab | ✓ | ball screw
QuakeLogic [20] | $\times$ | 15k | 50 | 0.5 | 1 | TESTLAB | ✓ | ball screw
Shakebot (ours) | ✓ | 1.3k | 2 | 0.5 | 1.2 | ROS | ✓ | toothed belt
| Note: Shakebot cost does not include camera or logistics. The cost of SARSAR
is estimated based on available items in the United States.
---
The Specifications of payload, velocity, and acceleration are the maximum
values. SARSAR has 200 $kg$ maximum payload at 0.4 $g$ acceleration.

In an effort to reduce variability, explore a broad parameter space of ground
motions, and automate PBR overturning analyses, our previous study developed a
virtual shake robot in simulation [21]. The virtual shake robot utilized core
technologies in robotics, such as Robot Operating System (ROS), Gazebo
simulation toolbox, and Bullet physics engine. Although the virtual shake
robot allows us to conduct overturning experiments in simulation, the results
need to be validated through physical experiments. To complete this
validation, we designed and built a low-cost, small-scale shake table, called
Shakebot (Fig. 1).

The Shakebot provides a reverse method for simulation validation. Physics
engines have been used for robotic development. However, when robotic systems
are transitioned from simulations to the real world, simulation gaps may cause
critical failures (e.g., in the arena of reinforcement learning [22, 23]).
When using a physics engine for scientific studies, we also need to address
the simulation gap to calibrate and scale-up the simulation experiments. Using
the Shakebot, we test the overturning dynamics of down-scaled, 3D-printed
PBRs. We build the same PBR models in simulation and validate the overturning
dynamics using our virtual shake robot. Following this approach, we are able
to compare the overturning patterns obtained from the Shakebot and the virtual
shake robot, further allowing us to quantify the uncertainty in simulation.
The Shakebot serves as a validation tool for simulation studies, aiding in
learning more about the real-time dynamics of the PBRs and potentially
calibrating physics engine parameters to reduce the simulation gaps.

The Shakebot demonstrates interdisciplinary research in the use of automated
robotic systems for geoscience studies, and enhances the impact of robotics in
the natural sciences. Besides its contribution to seismology, the key
contributions of this study include:

  * •

Shakebot is the first tabletop shake table developed based on robotic concepts
(perception and motion systems) and tools (ROS). Because of using ROS, the
control systems are reused from our previous virtual shake robot, ensuring
that the ground motions (the motion of the bed) are consistent in physical and
simulation experiments. This consistency is vital for simulation validation.

  * •

The hardware mechanism is simplified by design. The Shakebot adopts a closed-
loop stepper motor with a toothed belt and pulley transmission mechanism. Most
of the mechanical parts and components are off-the-shelf, and only a few parts
are 3D-printed. The simplified design lowers costs and expedites prototyping.
Additionally, for safety anticipation, we employ a double-emergency mechanism
to allow system override in case of malfunction.

  * •

Compared with traditional tabletop shake tables that only use accelerometers,
our Shakebot also leverages a top-down view camera to estimate ground motion.
The perception system provides an option of fusing accelerations and camera-
based displacements for ground velocity estimation.

  * •

We have open-sourced our hardware and software to promote robotic applications
in seismic research and education: https://github.com/DREAMS-
lab/asu_shake_table. Despite focusing on PBR dynamics in this study, the
Shakebot provides an affordable and accessible platform for other earthquake
research and education.

The paper is structured as follows. In Sec. II, we discuss commercial and
academic shake tables available for civil engineering and seismology, and
highlight a few of their capabilities and limitations. From Sec. III to Sec.
VI, we describe the mechanical design, the perception system, the motion
system, and the calibration system. In Sec. VII, we discuss the results of our
experiments, followed by our conclusions and future directions in Sec. VIII.

##  II Related Work

Full- or large-scale shake tables and tabletop shake tables have different
applications for research and education. Full-scale shake tables are designed
to simulate earthquakes on full-scale buildings, bridges, and other large
structures [24]. Full-scale shake tables, which usually incorporate
acceleration and displacement sensors for ground motion estimations, are
expensive to build and maintain. The high cost limits their availability in
research and education. Full-scale shake tables may also involve safety
considerations in operation and damage prevention.

Tabletop shake tables, on the other hand, are typically used to test smaller
models and structure components. Because of their low cost, these tables are
also useful in educational settings to teach students about the behavior of
structures under earthquake conditions. TABLE I compares existing tabletop
shake tables and the Shakebot. Existing tabletop shake tables are
comparatively expensive [14, 15, 20] or only allow a fixed ground motion
pattern [16]. None of the existing tabletop shake tables have open-sourced
their hardware. Another limitation of existing shake tables is that they
solely employ accelerometers to measure ground motions.

##  III System Description

In this section, we discuss the mechanical design, motor selection strategy,
and safety mechanism. The system is designed to satisfy the ground motion
requirements—maximum dynamic acceleration of 11.8 $m/s^{2}$ (1.2 $g$), maximum
velocity of 0.5 $m/s$, and maximum displacement of 0.45 $m$, with maximum
specimen (e.g., PBR) payload mass of 2 $kg$.

###  III-A Hardware Mechanism

The two main components of the Shakebot are a chassis and a control box, as
shown in Fig. 1. The chassis is constructed of an extruded aluminum T-slot
frame that supports two linear shafts (625 $mm$ length for each shaft). A
carriage is mounted on the linear shafts with ball bearings, which reduce the
friction between the linear shafts and the carriage. A flatbed (pedestal) is
attached to the carriage to hold specimens. A stepper motor (i.e., NEMA 34)
actuates the carriage through a toothed belt and pulley transmission. The
stepper motor driver (CL86T), contained in the control box, drives the stepper
motor. The control box also houses the power supply (S-350-60), an emergency
stop button, a touchscreen for the user interface, and a Raspberry Pi that
handles trajectory generation and low-level system control.

###  III-B Motor Selection

To select a stepper motor, we considered the stepper motor torque at its
maximum angular velocity to satisfy the required ground motion acceleration
and velocity. Because stepper motor torque decreases with velocity, if the
torque at the maximum velocity satisfies the acceleration requirement, the
stepper motor can also achieve so at any smaller velocities. We selected Nema
34HS31, a closed-loop stepper motor (with an encoder) with a torque output of
1.56 $Nm$ at the maximum speed of 1200 RPM. The maximum translational force is
calculated from the torque,

$\displaystyle F=ma=\frac{\tau}{r}$ (1)

where $F$ is the translational force, $r$ is the outside radius of the toothed
pulley, $m$ is the maximum payload mass (including carriage and specimen), $a$
is the maximum translational acceleration, and $\tau$ is the stepper motor
torque. We select a toothed pulley with an outside radius of 25.91 $mm$. The
translational force at the maximum speed is 60.21 $N$. Considering a maximum
total payload mass $m$ of 4 $kg$ (including a 2 $kg$ carriage and 2 $kg$
specimen), the corresponding acceleration produced from the translational
motion is 15.1 $m/s^{2}$ (1.54 $g$). The ratio of the calculated acceleration
to the required acceleration is 1.28, which is adequate empirically to
compensate for the friction and energy losses in transmission. Additionally,
at the maximum stepper motor speed of 1200 RPM, the translational speed can
reach 0.52 $m/s$, which also satisfies our design requirement for PBR studies.
Note that we used the outside radius on the specifications of the toothed
pulley for the calculation here. The outside radius value is appropriate for
the stepper motor selection. We conducted a calibration process that uses a
camera and several fiducial markers to estimate this radius precisely for
accurate control. The calibration process is discussed in Sec. VI-B.

The stepper motor driver reads signals from the encoder attached to the
stepper motor, forming a closed-loop system to control the stepper motor. The
microstepping technology in the stepper motor driver allows a higher step
resolution (1.8°/2000 per step), which achieves smooth motions at low speeds
(as low as 0.001 $m/s$ in our implementation).

###  III-C Safety Mechanism

The design incorporates two levels of safety features. The first level
consists of a pair of limit emergency switches situated at either end of the
carriage’s range of motion, as shown in Fig. 2. These switches link to the
emergency stop function on the stepper motor driver. Once either of these
switches is triggered, the stepper motor is disabled. The second level of
protection is supplied by an emergency push button in the control box, which
disconnects the power supply from the stepper motor driver. The emergency push
button can be manually activated by the Shakebot’s operator. Fig. 2 shows a
top-down view of the Shakebot.

Figure 2: Shakebot from top-down view illustrating double emergency mechanism.

##  IV Perception System

The perception system estimates the status of the bed using a camera and an
accelerometer. Four fiducial markers are affixed to the corners of the bed
(Fig. 1(b) and Fig. 2). A top-down view camera (FLIR Chameleon®3 Color Camera
CM3-U3-13Y3C-CS 1/2") captures the position and orientation of the bed using
the ROS package apriltag_ros [25, 26]. Additionally, an accelerometer (Wit
Motion HWT905-TTL) attached underneath the bed measures acceleration (ground
motion acceleration). We directly obtain ground motion displacement and
acceleration from the fiducial marker detection and accelerometer,
respectively. We additionally estimate the ground motion velocities by fusing
the information from the fiducial marker detection and accelerometer.

###  IV-A Displacement Estimation

The camera estimates ground motion displacement by detecting the fiducial
markers at each corner of the bed. The apriltag_ros package provides the pose
of each visible marker in the camera coordinates. When the bed moves from one
position to another, we obtain the relative pose of a marker attached to the
bed,

$H_{i}^{j}=H_{c}^{j}\cdot\ H_{i}^{c}=[H_{j}^{c}]^{-1}\cdot\ H_{i}^{c}$ (2)

where $i$ and $j$ indicate marker positions at two time stamps, $c$ indicates
the camera coordinates, $H_{a}^{b}$ is the transformation matrix for position
$a$ with respect to position $b$. $H_{j}^{c}$ and $H_{i}^{c}$ are obtained
using the apriltag_ros package.

We extract the translation vector from $H_{i}^{j}$ as the relative
displacement of the fiducial marker at two positions. We estimate ground
motion displacement by averaging the relative displacements from all visible
fiducial marker detections. When a fiducial marker is obstructed by a
specimen, the ground motion displacement estimation only relies upon the
remaining visible fiducial markers. Note that the relative displacement from
apriltag_ros is inaccurate despite intrinsic camera calibration. To address
this issue, perception calibration is performed (Sec. VI-A).

###  IV-B Acceleration Estimation

The accelerometer installed beneath the bed measures the ground motion
acceleration vector in three dimensions. However, the bed’s motion is limited
to one degree-of-freedom prismatic movement. To obtain the acceleration along
the bed movement direction, we align the accelerometer’s x-axis with the bed’s
direction of motion. However, the manual alignment is not perfect. We
additionally calculate the absolute value (Euclidean distance) of the
acceleration vector, because the absolute acceleration value is the same value
of the acceleration vector along the bed movement direction. To eliminate
high-frequency vibration, we apply a low-pass filter to the raw acceleration
data.

###  IV-C Velocity Estimation

We estimate ground motion velocity by fusing the derived velocities from the
fiducial marker detection and accelerometer data. Fig. 3 illustrates the
velocity fusion process. During a shake experiment, the timestamped ground
motion displacements and accelerations are recorded. After the shake
experiment, we implement numerical derivation to obtain velocity $v_{d}(t)$
from the ground motion displacements d(t). From the ground motion
accelerations $a(t)$, we conduct numerical integration to obtain velocity
$v_{a}(t)$.

Figure 3: Velocity fusion workflow. Estimated velocity is regression of
velocities from displacement and acceleration.

From the recordings of fiducial markers and acceleremeter, we obtain two sets
of derived velocities: $\mathbf{v_{d}}=\\{v_{d}(t_{1}),...,v_{d}(t_{M})\\}$
and $\mathbf{v_{a}}=\\{v_{a}(t_{1}),...,v_{a}(t_{N})\\}$. We phrase the ground
motion velocity estimation as a regression problem with an objective function
that minimizes the square of errors between estimated and derived velocities,

$\arg\min_{\theta}||\hat{v}(t;\theta)-\mathbf{v_{s}}(t)||$ (3)

where $\hat{v}(t;\theta)$ is a parametric ground motion velocity function, and
$\mathbf{v_{s}}=\mathbf{v_{d}}\cup\mathbf{v_{a}}$ is a set of velocities
derived from measurements. In this ground motion velocity estimation setting,
$\\{v_{d}(t)\\}$ and $\\{v_{a}(t)\\}$ need not be synchronized, simplifying
the electronic hardware configuration of our Shakebot. The fiducial marker
detection and accelerometer have different frame rates. The parametric
regression function $\hat{v}(t;\theta)$ has many options. We use a 6rd order
polynomial function for demonstration in Sec. VII.

##  V Motion System

The Shakebot has a hierarchical motion system and a user interface (UI) to
facilitate shake experiments. The hierarchical motion system provides two
options to emulate ground motions: single-pulse cosine displacement motion and
realistic ground motion from seismogram acceleration recordings. The UI
assists users in repeating the shake experiment where single-pulse cosine
displacement motions are deployed with different ground motion parameters.

Figure 4: Motion system supports two types of ground motions: single-pulse
cosine displacement and realistic ground motion from seismogram acceleration
recordings. Single-pulse cosine displacement is defined by peak ground
velocity (PGV) and peak ground acceleration (PGA).

###  V-A Motor Control

The Shakebot motion is enabled by a hierarchical motion planner in the
Raspberry Pi, which outputs PWM signals to the closed-loop stepper motor
driver. Fig. 4 shows the hierarchical motion system. The first ground motion
option is a single-pulse cosine displacement,

$d(t)=A-A\cos(2\pi ft),$ (4)

where $d(t)$ is the ground displacement function, $A$ is the amplitude, $f$ is
the frequency, and $t\in[0,1/f]$ is time. $A$ and $f$ are derived from ground
motion parameters,

$\displaystyle f=\frac{1}{2\pi\kappa}$ (5) $\displaystyle A=\frac{\alpha
g}{4\pi^{2}f^{2}}$ (6)

where $\alpha$ is the peak ground acceleration (PGA), $\kappa$ is the ratio of
peak ground velocity to peak ground acceleration (PGV/PGA), and $g$ is the
gravitational acceleration. As shown in Fig. 4, the motion interpreter
converts (PGV, PGA) to ($A,f$) using Eqs. 5-6. With ($A,f$), the trajectory
planner takes derivative of the cosine displacement function (Eq. 4) to obtain
ground velocity function,

$v(t)=2\pi Af\sin(2\pi ft)$ (7)

where $v(t)$ is the ground velocity function, and $t\in[0,1/f]$ is time. From
the ground velocity function, we uniformly sample velocity points $\\{v\\}$ as
the input of the velocity controller. The sampling frequency is a user-defined
parameter (usually 200 Hz). Note that the velocity commands $\\{v\\}$ are
desired translational velocities of the bed. The velocity controller converts
the translation velocity to the angular velocity of the stepper motor,

$\omega(t)=\frac{v(t)}{r}$ (8)

where $\omega(t)$ is angular velocity for the stepper motor, and $r$ is the
radius of the toothed pulley. The angular velocity commands are further
converted to PWM signals for the input of the stepper motor driver. The
stepper motor driver, stepper motor, and encoder attached to the stepper motor
form a closed-loop control to execute the PWM signals.

Besides the cosine ground displacement function, the Shakebot also supports
ground motion from observed seismograms. As shown in Fig. 4, we first
implement a low-pass filter to remove the high-frequency noise in the raw
acceleration data. Numerical integration computes velocities from the
accelerations. Then we have a high-pass filter to remove the low-frequency
noise in the velocities, because the low-frequency noise in velocities may
cause accumulated displacement errors. The output of the high-pass filter is a
set of velocities $\\{v\\}$, which have the same format as the output from the
trajectory planner. Because both the trajectory planner and high-pass filter
had the same output format, we share a velocity controller to process the
desired velocity commands.

###  V-B Control User Interface

The UI is a portal through which users conduct shake experiments to study PBR
overturning responses. Once initialized, the UI prompts the user to place a
PBR on the bed. Then the user needs to input (PGV/PGA, PGA), which is
converted to (PGV, PGA). The motion module (Fig. 4) executes a single-pulse
cosine displacement ground motion based on the (PGV, PGA). After the ground
motion is completed, the UI prompts the user to input the overturning response
of the PBR (i.e., being toppled or remaining balanced). The overturning
responses with ground motion parameters and perception estimates are saved in
a CSV file. Based on the user’s input, the UI exits or continues to conduct a
new overturning experiment.

##  VI System Calibration

The Shakebot needs to be calibrated before conducting the shake experiments.
We first calibrate the camera intrinsics and fiducial marker detection, and
this calibrated perception system is used to facilitate velocity controller
calibration.

###  VI-A Perception System Calibration

The perception system requires intrinsic camera calibration and fiducial
marker detection scaling. We used the camera_calibration ROS package [27] to
obtain the intrinsic camera parameters including scale factor, focal length,
pixel dimension, and distortion. Despite the intrinsic camera calibration, the
translation estimation from apriltag_ros was found to be inaccurate in our
implementation. To calibrate this fiducial marker detection process, we
estimate a multiplicative factor between the fiducial marker displacement and
manually measured displacement,

$D=\sigma d$ (9)

where $D$ is the manually measured displacement, $\sigma$ is the
multiplication factor, and $d$ is the displacement from fiducial marker
detection.

To better estimate $\sigma$, we have an overdetermined system by sampling a
number of displacement points along the linear shaft $\\{d\\}$,

$\bm{D}=\sigma\bm{d}$ (10)

where $\bm{D}=[D_{1},D_{2},\cdots,D_{N}]^{T}$ is a vector of manually measured
displacements, and $\bm{d}=[d_{1},d_{2},\cdots,d_{N}]^{T}$ is a vector of
fiducial marker displacements. The Moore–Penrose inverse is applied to solve
the overdetermined system,

$\sigma=\bm{D}\bm{d}^{T}[\bm{d}\bm{d}^{T}]^{-1}$ (11)

We utilize the limit switches at the ends of the linear shaft (i.e., bed
position calibration switches shown in Fig. 2) to facilitate the measurement
of bed displacement, $\bm{D}$. To begin the perception system calibration, the
bed moves to the left end until the left bed position calibration switch is
triggered. At the same time, the fiducial marker detection at this position is
recorded. The bed then moves to the right end, triggering the right
calibration switch, and the fiducial marker detection is recorded again. We
input the manually measured length between the two calibration switches.
During the perception calibration, we have a position controller that rotates
the motor shaft by one angular step given one PWM signal. As the bed moves
from one end to the other, we count and record the motor steps, from which we
calculate the bed displacement per angular step (translational step
resolution). Then, the bed randomly moves along the linear shaft to sample
data points, such as step counts (from the left calibration switch) and
fiducial marker detections. Using the step counts and translation step
resolution, we calculate the actual bed displacements. The samples of the
actual bed displacements and fiducial marker displacements are used in Eq. 11
to estimate $\sigma$.

###  VI-B Velocity Controller Calibration

The velocity controller calibration aims to reduce the uncertainties from the
toothed belt transmission. Because of the closed-loop stepper motor, we assume
that the error between the actual angular velocity and the desired angular
velocity of the stepper motor is diminished. Given a desired translational
velocity of the bed, Eq. 8 calculates the corresponding desired angular
velocity. However, the toothed pulley radius $r$ is close but not equal to the
outside radius of the pulley. The uncertainty in $r$ comes from factors such
as the tightness of the toothed belt, manufacturing imprecision, and tooth
geometry. We correct for such uncertainty using a multiplicative factor,

$\omega^{\prime}(t)=\frac{\gamma v(t)}{r}$ (12)

where $\gamma$ is a constant that captures the first-order systematic
uncertainty in the toothed belt transmission. The velocity controller
calibration needs to estimate a $\gamma$ that minimizes the error between the
desired and actual translational velocities. Because $\gamma$ is a factor that
calibrates velocity, we denote $\gamma v(t)$ as calibrated velocity. Instead
of directly measuring translational velocity in real time, we conduct half-
cosine displacement motions (i.e., $t\in[0,1/(2f)]$ in Eq. 4) and estimate
$\gamma$ from displacement differences,

$D^{\prime}=\gamma D$ (13)

where $D$ is the desired half-cosine displacement, and $D^{\prime}$ is the
actual half-cosine displacement measured from the fiducial marker detection.
Similarly, we randomly sample the half-cosine displacement motions and use the
Moore–Penrose inverse to estimate $\gamma$. Using displacement difference to
estimate velocity difference has the advantage of easy implementation, because
directly measuring translation velocity is challenging.

We briefly demonstrate that $\gamma$ estimated from the displacements is the
same multiplicative factor for velocity calibration. Based on Eq. 13, we have
calibrated displacement function,

$D^{\prime}(t)=\gamma A-\gamma A\cos(2\pi ft).$ (14)

Note that $\gamma$ here is estimated from the displacement differences
described above. We take the derivative of $D^{\prime}(t)$ and obtain the
calibrated velocity function,

$v^{\prime}(t)=\dot{D^{\prime}}(t)=2\pi\gamma Af\sin(2\pi ft).$ (15)

We substitute Eq. 7 into Eq. 15,

$v^{\prime}(t)=\gamma v(t).$ (16)

Therefore, $\gamma$ estimated from displacements is also the factor that
calibrates velocity.

##  VII Experiments and Discussion

We conducted experiments to demonstrate an application of the Shakebot in PBR
fragility studies and to test its performance in validating simulation. These
experiments focused on a PBR at the Double Rock site in coastal Central
California [11]. The PBR, as shown in Fig. 5, was mapped by an unpiloted
aerial vehicle and constructed using structure from motion and Poisson
reconstruction algorithms [21].

(a) Double Rock site in coastal Central California

(b) PBR extracted from Double Rock site

Figure 5: Precariously Balanced Rock (PBR) at Double Rock site in coastal
Central California. The yaw orientation indicates the relative direction of
single-pulse cosine displacement motion.

The first experiment demonstrates the Shakebot application in studying the
fragility anisotropy of PBRs. We down-scaled the Double Rock PBR from 151.0 cm
to 12.0 cm in height and 3D-printed it with PLA material. The resulting PBR,
denoted as 3D-printed PBR(a), had a total mass of 105g. Because of the
asymmetric geometry of 3D-printed PBR(a), the overturning responses were
expected to vary with different ground motion directions. To investigate this
fragility anisotropy idea, we placed the 3D-printed PBR(a) with two initial
orientations (yaw angles of 0° and 270°, as shown in Fig. 5) on the bed. Using
the UI, we obtained the response diagrams from a set of ground motions, as
shown in Fig. 6. The response diagrams indicated that 3D-printed PBR(a)
oriented at a yaw angle of 0° was more fragile than at yaw angle of 270°. This
result is consistent with previous studies showing that PBRs with smaller
minimal contact angles along the motion direction are more fragile [5, 28].
Addtionally, the resulting boundary curves closely resembled the simulation
results of the Double Rock PBR with original shape [21].

(a) Yaw 0°

(b) Yaw 270°

Figure 6: Overturning response diagrams of the Double Rock PBR from
orientations of yaw 0° and yaw 270°. The red and blue dots represent
overturning responses of being toppled and balanced after single-pulse cosine
ground motions, respectively. The curves indicate the boundaries of the
overturning responses.

In the second experiment, we used the Shakebot to evaluate the overturning
simulation. The Double Rock PBR was down-scaled to 12.8 cm in height and
3D-printed with the PETG material. The resulting PBR, denoted as 3D-printed
PBR(b), had a total mass of 404 g. We imported the 3D model of the down-scaled
PBR to Gazebo and conducted overturning experiments using the virtual shake
robot [21]. To examine the simulation results, we used the Shakebot to
overturn 3D-printed PBR(b). In both simulation and real-world experiments, the
PBR had a yaw angle of 0°. As shown in Fig. 7, the overturning boundary curves
from the real-world and simulation experiments are similar. Future work should
explore the causes of their difference. For example, we observed that Gazebo
reduced the precision of configuration parameters when passing those
parameters to the Bullet physics engine. The data from real-world experiments
will aid in calibrating and improving the simulation experiments.

Figure 7: Overturning boundary curves from Shakebot and virtual shake robot.
Reprinted from [21]. Reprinted with permission.

Fig. 8 shows the estimation results of acceleration, displacement, and
velocity from the input PGA of 0.98 $m/s^{2}$ and PGV/PGA of 0.1 $s$.
According to the PGA and PGV/PGA, the desired PGV and peak ground displacement
are 0.1 $m/s$ and 0.01 $m$, respectively. The acceleration estimation and
displacement estimation can be improved with better measurement equipment such
as a better-quality accelerator and motion capture system. Our velocity
estimation method used the 6rd order polynomial function to fuse data from
afford sensors, resulting in a PGV estimate of 0.12 $m/s$. The resulting
velocity estimation is better than individual velocities derived from the
accelerometer and fiducial marker detection. Future work should quantify
uncertainties in displacement, velocity, and acceleration estimates.
Investigations on other regression models (e.g., Gaussian process regression)
and real-time state estimation may also improve the ground motion estimation.

Because validating ground motion estimation requires future investment in
expensive measurement equipment, we assumed that the actual ground motions
were the same as the desired ground motions during the overturning experiment.
This assumption is reasonable given our hierarchical control system and
calibration processes. However, future work should verify this assumption. We
also anticipate that the PBR response status registration can be automated by
an object detection system and PBR placement can be implemented by a robotic
arm. By automating the overturning experiments, we can save on human labor and
ensure accurate PBR placement for reliable results.

Figure 8: Velocity estimation from regression of acceleration and
displacement.

##  VIII Conclusion

We introduced a low-cost, open-source shake table named Shakebot for
earthquake research and education. Leveraging concepts and tools from
robotics, we simplified and advanced the hardware and software compared to the
existing shake tables. We presented the technical details of the perception
and motion systems, which are important for the open-source community to
upgrade or customize the Shakebot. Experiments showed the applications of the
Shakebot for PBR studies. Our work demonstrates an interdisciplinary method of
using robotics for natural sciences, providing a step towards robotics-enabled
automated geoscience [21]. The Shakebot provides an affordable platform for
education, outreach, and other ground motion experiments.

## ACKNOWLEDGEMENTS

This research was supported by the Pacific Gas and Electric Company and
Southern California Earthquake Center (Contribution No. 21137). SCEC is funded
by NSF Cooperative Agreement EAR-1600087 and USGS Cooperative Agreement
G17AC00047. Thank you to Chris Madugo for his support and advice in the
research.

## References

  * [1] Lydik S Jacobsen  “Vibration research at Stanford university”  In _Bulletin of the Seismological Society of America_ 19.1 The Seismological Society of America, 1929, pp. 1–27
  * [2] Xilin Lu, Gongkang Fu, Weixing Shi and Wensheng Lu  “Shake table model testing and its application”  In _The Structural Design of Tall and Special Buildings_ 17.1 Wiley Online Library, 2008, pp. 181–201
  * [3] Aleksandra Bogdanovic, Zoran Rakicevic and Ehsan Noroozinejad Farsangi  “Shake table tests and numerical investigation of a resilient damping device for seismic response control of building structures”  In _Structural Control and Health Monitoring_ 26.11 Wiley Online Library, 2019, pp. e2443
  * [4] Zoran Rakicevic, Aleksandra Bogdanovic, Ehsan Noroozinejad Farsangi and Abbas Sivandi-Pour  “A hybrid seismic isolation system toward more resilient structures: Shaking table experiment and fragility analysis”  In _Journal of Building Engineering_ 38 Elsevier, 2021, pp. 102194
  * [5] Matthew D Purvance, Abdolrasool Anooshehpoor and James N Brune  “Freestanding block overturning fragilities: Numerical simulation and experimental validation”  In _Earthquake Engineering & Structural Dynamics_ 37.5 Wiley Online Library, 2008, pp. 791–808
  * [6] Ahmed Elgamal, Michael Fraser and Flora McMartin  “On-line educational shake table experiments”  In _Journal of Professional Issues in Engineering Education and Practice_ 131.1 American Society of Civil Engineers, 2005, pp. 41–49
  * [7] James N Brune  “Precariously balanced rocks and ground-motion maps for southern California”  In _Bulletin of the Seismological Society of America_ 86.1A The Seismological Society of America, 1996, pp. 43–54
  * [8] Baoping Shi, Abdolrasool Anooshehpoor, Yuehua Zeng and James N Brune  “Rocking and overturning of precariously balanced rocks by earthquakes”  In _Bulletin of the Seismological Society of America_ 86.5 The Seismological Society of America, 1996, pp. 1364–1371
  * [9] Abdolrasool Anooshehpoor, James N Brune and Yuehua Zeng  “Methodology for obtaining constraints on ground motion from precariously balanced rocks”  In _Bulletin of the Seismological Society of America_ 94.1 Seismological Society of America, 2004, pp. 285–303
  * [10] John G Anderson, Glenn P Biasi and James N Brune  “Precarious rocks: Providing upper limits on past ground shaking from earthquakes”  In _Earthquake hazard, risk and disasters_ Elsevier, 2014, pp. 377–403
  * [11] AH Rood et al.  “Earthquake hazard uncertainties improved using precariously balanced rocks”  In _AGU Advances_ 1.4 Wiley Online Library, 2020, pp. e2020AV000182
  * [12] Chik-Sing Yim, Anil K. Chopra and Joseph Penzien  “Rocking response of rigid blocks to earthquakes”  In _Earthquake Engineering & Structural Dynamics_ 8.6, 1980, pp. 565–587
  * [13] M Khalid Saifullah and Christine E Wittich  “Seismic response of two freestanding statue-pedestal systems during the 2014 South Napa Earthquake”  In _Journal of Earthquake Engineering_ 26.10 Taylor & Francis, 2022, pp. 5086–5108
  * [14] Quanser  “Quanser Shake Table II”  URL: https://www.quanser.com/products/shake-table-ii/
  * [15] T Baran, AK Tanrikulu, C Dundar and AH Tanrikulu  “CONSTRUCTION AND PERFORMANCE TEST OF A LOW-COST SHAKE TABLE”  In _Experimental Techniques_ 35.4 Blackwell Publishing Ltd Oxford, UK, 2011, pp. 8–16
  * [16] Kelvin Educational  “Kelvin advanced earthquake simulator”  URL: https://kelvin.com/kelvin-e-q-machine-2021/
  * [17] E. Damcı and Ç. Şekerci  “Development of a Low-Cost Single-Axis Shake Table Based on Arduino”  In _Experimental Techniques_ 43.2, 2019, pp. 179–198
  * [18] X Shao and G Enyart  “Development of a versatile hybrid testing system for seismic experimentation”  In _Experimental Techniques_ 38.6 Springer, 2014, pp. 44–60
  * [19] Gökçe Kınay  “Construction and control of a desktop earthquake simulator”, 2006
  * [20] QUAKELOGIC  “Testbox-shaketable”  URL: https://quakelogic.net/actuators-shaketables/testbox_shaketable
  * [21] Zhiang Chen  “Automated Geoscience with Robotics and Machine Learning: A New Hammer of Rock Detection, Mapping, and Dynamics Analysis”, 2022
  * [22] Wenshuai Zhao, Jorge Peña Queralta and Tomi Westerlund  “Sim-to-real transfer in deep reinforcement learning for robotics: a survey”  In _2020 IEEE Symposium Series on Computational Intelligence (SSCI)_ , 2020, pp. 737–744  IEEE
  * [23] Zihan Ding and Hao Dong  “Challenges of reinforcement learning”  In _Deep Reinforcement Learning_ Springer, 2020, pp. 249–272
  * [24] Sofia Gavridou et al.  “Shake-table test of a full-scale 4-story precast concrete building. I: Overview and experimental results”  In _Journal of Structural Engineering_ 143.6 American Society of Civil Engineers, 2017, pp. 04017034
  * [25] Danylo Malyuta  “Guidance, Navigation, Control and Mission Logic for Quadrotor Full-cycle Autonomy”, 2017
  * [26] John Wang and Edwin Olson  “AprilTag 2: Efficient and robust fiducial detection”  In _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ IEEE, 2016, pp. 4193–4198
  * [27] ROS.org  “camera_calibration”  URL: http://wiki.ros.org/camera_calibration
  * [28] David E Haddad et al.  “Estimating Two-dimensional Static Stabilities and Geomorphic Settings of Precariously Balanced Rocks from Unconstrained Digital Photographs”  In _Geosphere_ 8.5 Geological Society of America, 2012, pp. 1042–1053


# Ion Correlation Induced Non-monotonic Height Change and Microphase
Separation of Polyelectrolyte Brushes

Chao Duan  Department of Chemical and Biomolecular Engineering, University of
California Berkeley, CA 94720, USA  Nikhil R. Agrawal  Department of Chemical
and Biomolecular Engineering, University of California Berkeley, CA 94720, USA
Rui Wang<EMAIL_ADDRESS>Department of Chemical and Biomolecular
Engineering, University of California Berkeley, CA 94720, USA  Materials
Sciences Division, Lawrence Berkeley National Lab, Berkeley, CA 94720, USA

(August 28, 2024)

###### Abstract

Modeling ion correlations in inhomogeneous polymers and soft matters with
spatially varying ionic strength or dielectric permittivity remains a great
challenge. Here, we develop a new theory which systematically incorporates
electrostatic fluctuations into the self-consistent field theory for polymers.
Applied to polyelectrolyte brushes, the theory predicts that ion correlations
induce non-monotonic change of the brush height: collapse followed by
reexpansion. The scaling analysis elucidates the competition between the
repulsive osmotic pressure due to translational entropy and the attraction
induced by ion correlations. We also clarify the absence of causal
relationship between the brush collapse-reexpansion and the inversion of the
surface electrostatic potential. Furthermore, strong ion correlations can
trigger microphase separation, either in the lateral direction as pinned
micelles or in the normal direction as oscillatory layers. Our theoretical
predictions are in good agreement with the experimental results reported in
the literature.

Ion-containing polymers and soft matters with spatially varying ionic strength
or dielectric permittivity are ubiquitous [1, 2, 3, 4]. One example is
polyelectrolyte (PE) brushes grafted on substrates which regulate surface
properties such as transport, wettability, adhesion, antifouling and
lubrication [5, 6, 7, 8, 9, 10]. PE brushes are very sensitive to the addition
of salt ions [11, 12]. The primary role of monovalent ions, screening and
translational entropy, is well captured by the mean-field Poisson-Boltzmann
theory [13, 14] . However, brush morphology and mechanical behaviors in
multivalent salt solutions are quite non-trivial [15, 16]. Using neutron
reflectivity and surface forces apparatus, Tirrell group found that the onset
of brush collapse is reduced by 3-6 orders of magnitudes in multivalent ions
compared to monovalent ions [17]. It was observed that brush height changes
non-monotonically: collapse followed by reexpansion as salt concentration
further increases. Similar non-monotonic size change was reported in the study
of single PE conformation [18]. AFM images showed that trivalent ions induce
lateral structural inhomogeneities of PE brushes, which were further confirmed
by molecular simulations [19, 20]. Furthermore, Yu et al. also reported that
the friction force between two opposing PE brushes increases dramatically in
multivalent salt solutions, diminishing the lubricity [21]. These phenomena
cannot be even qualitatively explained by mean-field electrostatics.

An obvious factor missing in the mean-field theories is the electrostatic
correlation [22, 23, 24]. Despite many efforts, it remains a great challenge
to model the coupling between electrostatic correlation and chain conformation
particularly in inhomogenous systems. Muthukumar developed a variational
method based on a two-state model for counterion binding [25, 26]. This theory
pre-integrates the salt degrees of freedom, resulting in a screened
interactions between charged monomers, which ignores the feedback on the
correlations between ions. Applying this method to PE brushes, Li et al.
explained the sharp collapse in multivalent solutions but failed to capture
the subsequent reexpansion [27, 28]. Shen and Wang addressed the electrostatic
fluctuations coupled with chain connectivity through quantifying the self
energy of a PE [29, 30]. Recently, Sing and Qin treated the correlation via
cluster expansions [31]. However, all these three approaches invoked
preassumed functional form of chain conformation, limiting their application
to PEs with various structures like globule and pearl-necklace [32, 33]. These
theories are also not applicable to inhomogeneous systems. Furthermore, a
hybrid liquid state integral equation and self-consistent field theory was
developed to model correlations in inhomogenous polymer blends and block
copolymers [34, 35]. Local charge neutrality is assumed, which ignores long-
range components of electrostatics. This assumption fails wherever a double
layer with local charge-separation exist, e.g. PE adsorption and brushes [36,
37, 38, 39, 15, 16, 17, 19, 21, 20, 18].

We develop a new theory which systematically incorporates the variational
Gaussian renormalized fluctuation theory for electrostatics [40, 41, 42, 43,
44] into the self-consistent field theory (SCFT) for polymers [33, 45, 46].
The coupling between ion correlation and chain conformation can be fully
captured without any preassumption on the PE structure. Our theory is
particularly capable of modeling inhomogeneous PEs with spatially varying
ionic strength or dielectric permittivity and can be performed in the high-
dimensional calculation. Applied to PE brushes, the theory predicts the non-
monotonical salt concentration dependence of brush height and lateral
microphase separation, in good agreement with experiments.

We consider a general system consisting of $n_{P}$ negatively charged PEs each
with $N_{P}$ monomers in an electrolyte solution which contains $n_{S}$
solvent molecules, $n_{C}$ monovalent counterions ($z_{C}=+1$) and $n_{\pm}$
salt ions with valency $z_{\pm}$. The system is taken as a semicanonical
ensemble: the number of PEs and counterions are fixed whereas solvents and
salt ions are connected with a bulk salt solution of ion concentration
$\rho^{b}_{\pm}$ that maintains the chemical potentials of $\mu_{S}$ and
$\mu_{\pm}$. To more accurately describe the charge interactions, the ionic
charge is modeled by a finite spread function $h_{K}({\bf r},{\bf
r}^{\prime})$ ($K=P,C,\pm$) instead of using the point-charge model[40, 42,
41, 43, 44]. For convenience, we take $h_{K}({\bf r},{\bf r}^{\prime})$ to be
Gaussian which retains the Born radius $a_{K}$. The smeared charge model for
PE widely used in existing theories fails to capture any correlations between
charged monomers as well as between them and mobile ions. We adopt a discrete
Gaussian chain model where each monomer is a charged particle with valency
$z_{P}$. The number of dissociated counterions is thus
$n_{C}=-z_{P}n_{P}N_{P}$. The partition function is

$\displaystyle\mathcal{Z}$
$\displaystyle=\frac{1}{n_{P}!v_{P}^{n_{P}N_{P}}}\prod_{i=1}^{n_{P}}\int\hat{D}\\{{\bf
R}_{i}\\}\frac{1}{n_{C}!v_{C}^{n_{C}}}\prod_{j=1}^{n_{C}}\int d{\bf r}_{C,j}$
$\displaystyle\cdot\sum_{\gamma=S,\pm}\sum_{n_{\gamma}=0}^{\infty}\frac{e^{\mu_{\gamma}n_{\gamma}}}{n_{\gamma}!v_{\gamma}^{n_{\gamma}}}\prod_{k=1}^{n_{\gamma}}\int
d{\bf r}_{\gamma,k}{\rm exp}\left(-H\right)$ (1)

where $v_{P}$ and $v_{S}$ are the volume of the PE monomers and small
molecules, respectively. We assume $v_{P}=v_{S}=v$. $\int\hat{D}\\{{\bf
R}_{i}\\}$ denotes integration over all configurations weighted by the
Gaussian-chain statistics [45]. The integrals in Eq. 1 are under the
incompressibility constraint $\sum_{K}\hat{\rho}_{K}v_{K}=1$ for all species
where $\hat{\rho}_{K}$ is the instantaneous number density of species $K$. The
Hamiltonian

$\displaystyle H=\chi v\int d{\bf r}\hat{\rho}_{P}({\bf r})\hat{\rho}_{S}({\bf
r})+\frac{1}{2}\int d{\bf r}d{\bf r}^{\prime}\hat{\rho}_{e}({\bf r})C({\bf
r},{\bf r}^{\prime})\hat{\rho}_{e}({\bf r}^{\prime})$ (2)

includes hydrophobic interaction between polymer and solvent in terms of
Flory-Huggins $\chi$ parameter and the long-range electrostatic interaction
between all charged species. $\hat{\rho}_{e}({\bf r})=\sum_{K}\int d{\bf
r}^{\prime}z_{K}h_{K}({\bf r}-{\bf r}^{\prime})\hat{\rho}_{K}({\bf r})$ with
$K=P,C,\pm$ is the net charge density, and $C({\bf r},{\bf r}^{\prime})$ is
the Coulomb operator.

Following the standard field theoretical approach [45], the interacting system
can be decoupled into noninteracting chains and ions in fluctuating fields
(see the detailed derivation in Sec. I of the Supplemental Materials). For the
electric field, instead of saddle-point approximation, we use a non-
perturbative variational approach to account for fluctuation by introducing a
general Gaussian reference action [40, 42, 41, 43, 44, 29, 30]. The
consequence of this fluctuation is the self-energy of charged species. For
small mobile counter ions and salt ions, the self energy $u^{0}_{K}({\bf r})$
can be individually evaluated from their own charge spread [42, 41, 43, 44].
In contrast, the charge of a polymer spreads the pervaded volume of the entire
chain [29, 30]. The polymer self energy, when evaluated from its constituting
monomer $i$, includes both the individual contribution $u^{0}_{P}({\bf r})$
and an extra $u^{ex}_{P}({\bf r};i)$ attributed by all other intra-chain
monomers.

The key results of our theory are the self consistent equations that capture
the coupling between charge correlation (manifested by the self-energy
$u^{0}_{K}({\bf r})$ and $u^{ex}_{P}({\bf r};i)$) and chain conformation
(described by the propagator $q({\bf r};i)$ and its complementary
$q^{\dagger}({\bf r};i)$):

$\displaystyle u^{0}_{K}({\bf r})$ $\displaystyle=\frac{z^{2}_{K}}{2}\int
d{\bf r^{\prime}}d{\bf r^{\prime\prime}}h_{K}({\bf r},{\bf r}^{\prime})G({\bf
r}^{\prime},{\bf r}^{\prime\prime})h_{K}({\bf r^{\prime\prime}},{\bf r})$ (3a)
$\displaystyle u^{ex}_{P}({\bf
r};i)=\frac{z^{2}_{P}}{2}\biggl{[}\sum^{i-1}_{j=1}\frac{\int d{\bf
r^{\prime}}q_{0}({\bf r}^{\prime};j)e^{V^{0}_{P}({\bf r}^{\prime})}G({\bf
r}^{\prime},{\bf r})g^{\dagger}_{0}({\bf r}^{\prime},{\bf r};j,i)}{q_{0}({\bf
r};i)}$ $\displaystyle+\sum^{N_{P}}_{j=i+1}\frac{\int d{\bf
r^{\prime}}g_{0}({\bf r},{\bf r}^{\prime};i,j)G({\bf r},{\bf
r}^{\prime})e^{V^{0}_{P}({\bf r}^{\prime})}q_{0}^{\dagger}({\bf
r}^{\prime};j)}{q_{0}^{\dagger}({\bf r};i)}\biggl{]}$ (3b) $\displaystyle
q({\bf r};i)$ $\displaystyle=e^{-V_{P}({\bf r};i)}\int d{\bf
r^{\prime}}\Phi({\bf r}-{\bf r^{\prime}})q({\bf r^{\prime}};i-1)$ (4a)
$\displaystyle q^{\dagger}({\bf r};i)$ $\displaystyle=e^{-V_{P}({\bf
r};i)}\int d{\bf r^{\prime}}\Phi({\bf r}-{\bf r^{\prime}})q^{\dagger}({\bf
r^{\prime}};i+1)$ (4b)

$G({\bf r},{\bf r}^{\prime})$ in Eq. 3 is the charge correlation function:

$\displaystyle-\nabla_{\bf r}\cdot[\epsilon$ $\displaystyle({\bf
r})\nabla_{\bf r}G({\bf r},{\bf r}^{\prime})]+2I_{0}({\bf r})G({\bf r},{\bf
r}^{\prime})$ $\displaystyle+\int d{\bf r^{\prime\prime}}2I_{ex}({\bf r},{\bf
r^{\prime\prime}})G({\bf r^{\prime\prime}},{\bf r}^{\prime})=\delta({\bf
r}-{\bf r}^{\prime})$ (5)

where $\epsilon({\bf r})$ is the scaled permittivity and $I_{0}({\bf
r})=\sum_{K}z^{2}_{K}\rho_{K}({\bf r})/2$ is the local component of ionic
strength. Different from the commonly used Green function, the ionic strength
of PE systems has an extra non-local contribution $I_{ex}({\bf r},{\bf
r^{\prime}})$ from the long-range charge spread of polymer (see Sec. I in the
Supplemental Material). $\Phi({\bf r}-{\bf r^{\prime}})$ in Eq. 4 is the
Gaussian bond transition probability [45]. $V_{P}({\bf r};i)=V^{0}_{P}({\bf
r})+u^{ex}_{P}({\bf r};i)$ is the total interaction field experienced by
monomer $i$, where $V^{0}_{P}({\bf r})=\omega_{P}+z_{P}\psi+u^{0}_{P}$.
$\omega_{P}=\chi\phi_{S}+\xi-(\partial\epsilon/\partial\rho_{P})(\nabla\psi)^{2}/2$
is the field conjugate to monomer density with $\xi$ reinforcing the
incompressibility. The electrostatic potential $\psi$ satisfies the Poisson
equation $-\nabla\cdot[\epsilon({\bf r})\nabla\psi({\bf
r})]=\sum_{K}z_{K}\rho_{K}({\bf r})$. The monomer density $\rho_{P}({\bf
r})=(n_{P}/Q_{P})\sum_{i=1}^{N}q({\bf r};i)e^{V_{P}({\bf
r};i)}q^{\dagger}({\bf r};i)$ with $Q_{P}=\int d{\bf r}q({\bf r};N_{P})$.
Counterion density $\rho_{C}({\bf r})=n_{C}{\rm
exp}(-z_{C}\psi-u^{0}_{C}-v_{C}\xi/v)/Q_{C}$ with $Q_{C}=\int d{\bf r}{\rm
exp}(-z_{C}\psi-u^{0}_{C}-v_{C}\xi/v)$. Salt ion densities $\rho_{\pm}({\bf
r})=\lambda_{\pm}{\rm exp}(-z_{\pm}\psi-u^{0}_{\pm}-v_{\pm}\xi/v)$ with the
fugacity $\lambda_{\pm}=e^{\mu_{\pm}}/v_{\pm}$ determined by the bulk
concentration $\rho^{b}_{\pm}$.

When evaluating $u^{ex}_{P}({\bf r};i)$ in Eq. 3b, we expand the Hamiltonian
based on a reference state where charged monomers are treated isolatedly that
neglects charge correlations from other intra-chain monomers. The chain
propagator $q_{0}({\bf r};j)$ and intra-chain correlation function $g_{0}({\bf
r},{\bf r}^{\prime};i,j)$ are thus subject to the same recurrence relation
(Eq. 4) as $q$ but with interaction field $V^{0}_{P}$ corresponding to the
reference state.

For a given chain conformation, the correlation between charged species can be
calculated from Eq. 5, leading to an update of the self energy (Eq. 3). This
will in turn feedback to the polymer density and chain conformation through
Eq. 4. The coupling between ion correlations and chain conformation is hence
self-consistently achieved via this iteration. The numerical details are
provided in Sec. II of the Supplemental Material [47, 48].

The theory is general for all homogeneous and inhomogeneous PE systems. Here
we apply it to PE brushes where one end of each chain is tethered to a planar
surface. Without loss of generality, we consider negatively charged brushes
immersed in a $z_{+}$:1 salt solution with the bulk salt concentration
$\rho_{b}$. To focus on the ion correlation effects, we assume $\Theta$
solvent ($\chi$=0.5) and neglect the excluded volume of ions
($v_{C}=v_{\pm}=0$). The dielectric constant is uniformly set to be 80. The
contribution of intra-chain correlation to self-energy $u^{ex}_{P}$ is found
less important than the individual contribution of a charged monomer
$u^{0}_{P}$ under the condition of high ionic strengths (either high salt
concentration or high PE density) [29]. For systems we are interested in,
$u^{ex}_{P}$ is also neglected.

Figure 1: The effect of ion correlations on the density profile of PE brushes
in the normal direction in comparison with the mean-field results (dashed
lines). $\phi_{P}(z)$ is the polymer volume fraction. $N=100$, $b=1.0$nm,
$v=1.0$nm3, $z_{P}=z_{-}=-1$, the Born radius
$a_{P}=a_{C}=a_{\pm}=2.5\mathring{\rm A}$, and the grafting density
$\sigma=0.1$nm-2. Bulk ionic strength $I_{b}=0.3$M for all cases. Figure 2:
Non-monotonic change of brush height on salt concentration. $h$ is normalized
by its the salt-free value $h_{0}$. Lines represent our theoretical
predictions and scattering symbols represent experimental data adopted from
Ref. [17] with Na+ ($\bigcirc$), Ba2+ ($\square$) and Y3+ ($\triangle$). The
horizontal dashed line is for the charge neutral brushes. Critical salt
concentrations $\rho^{*}_{b}$, $\rho^{**}_{b}$, and $\rho^{***}_{b}$ are
identified on the $z_{+}=2$ curve. $\sigma$ in the calculations has the same
value $0.1$nm-2 as the experiments.

Electrostatic correlation has dramatic impacts on the brush morphology. Figure
1 plots polymer density profiles in the presence of cations with different
valency at the same bulk ionic strength. The results based on mean-field
electrostatics are also provided for comparison, which predict swollen brushes
in all cases and insensitivity of the profile on $z_{+}$. In the stark
contrast, with the inclusion of ion correlation, our theory predict that
brushes are highly collapsed in multivalent salt solutions. The brush height
is even smaller than the charge-neutral case and becomes more compact as
$z_{+}$ increases. This collapse originates from the net attractive force
between charged monomers induced by electrostatic correlation.

Our theory fully captures the non-monotonic height change of PE brush in
response to multivalent salt. Figure 2 plots the brush height $h$ as a
function of $\rho_{b}$, where $h$ is calculated from the Gibbs dividing
surface. For $z_{+}=1$, $h$ shows three regimes consistent with the mean-field
results: an “osmotic regime” at low $\rho_{b}$ where brushes are highly
swollen with a constant $h$, a “salted regime” at intermediate $\rho_{b}$
where the swollen brushes shrink due to salt screening, and a “neutral regime”
at high $\rho_{b}$ where the charges on brushes are almost screened out and
$h$ returns to the neutral level [49]. However, $h$ for $z_{+}=2$ and $3$ is
qualitatively different from the mean-field result. A new “collapsed” regime
appears after the salted regime, where $h$ drops below the neutral brush
level, reaches the minimum and reexpands to reenter the neutral regime. Our
theoretical predictions are in good agreement with the experimental data
reported by Tirrell group [17].

We perform the scaling analysis to elucidate the insights of the non-monotonic
height change (see Sec. III in the Supplemental Materials). It is built upon
the force balance in the brush region: $f_{\rm ela}+f_{\rm int}+f_{\rm ion}=0$
[49]. $f_{\rm ela}\approx-\sigma h/Nb^{2}$ is the elastic stress due to chain
deformation. $f_{\rm int}\approx w_{2}\rho^{2}_{P}+w_{3}\rho^{3}_{P}$ is the
non-electrostatic interaction between monomers. $\rho_{P}=\sigma N/h$; $w_{2}$
and $w_{3}$ are the second and third virial coefficients, respectively.
$f_{\rm ion}$ is the electrostatic osmotic pressure including ion
transnational entropy and ion correlation. In the dilute limit, $f_{\rm
ion}=\sum_{K=\pm}(\rho_{K}-\rho^{b}_{K})-\pi/24(\kappa^{3}-\kappa_{b}^{3})$
with $\kappa$ the inverse of the Debye length [40]. We identify the critical
salt concentrations which locate the boundaries between adjacent regimes (see
Fig. 2). The onset of the brush shrink $\rho^{*}_{b}$ occurs when the ion
correlation in the brush becomes comparable to the translational entropy (i.e.
$\kappa^{3}\approx\rho_{+}$), leading to the scaling of $\rho^{*}_{b}$ on the
cation valency as

$\rho^{*}_{b}\sim\begin{cases}1/\left[z_{+}(z_{+}+1)\right]^{3}&\text{(with
ion correlations)}\\\ 1/\left[z_{+}(z_{+}+1)\right]&\text{(mean-
field)}\end{cases}$ (1)

Compared with the mean-field result, it clearly indicates a much stronger
$z_{+}$ dependence of $\rho^{*}_{b}$ in the presence of ion correlations. This
scaling explains the experimental observation that the onset of brush shrink
occurs at the salt concentration a few orders of magnitudes lower for
multivalent ions compared to monovalent ions [17].

Ion correlations get stronger as $\rho_{b}$ increases. At the critical
concentration $\rho^{**}_{b}$, the correlation-induced attraction exactly
compensates the repulsive pressure from the translational entropy of ions,
which leads to $f_{\rm ion}=0$ and hence neutral brush apparently. In the
salted regime and thereafter, we work out that $f_{\rm
ion}\approx(\rho^{2}_{P}/z_{+}\rho_{b})\Gamma$, where $\Gamma$ can be
understood as the effective relative charge on the brush after counterion
adsorption. $\Gamma$ satisfies

$\displaystyle\frac{(z_{+}+1)\Gamma-1}{\Gamma-1}$
$\displaystyle=-\frac{z_{+}(z_{+}+1)\kappa_{b}}{16\pi\epsilon}\frac{du}{d(a\kappa_{b})}$
(7)

$u/(8\pi\epsilon a)$ is the self energy of monovalent ion with Born radius $a$
in the bulk salt solution [40]. $\rho^{**}_{b}$ is thus determined at
$\Gamma=0$ through Eq. 7. It should be noted that Eq. 7 is universal for all
charge spread models. If $h$ takes the Gaussian form particularly,
$u=1-(a\kappa_{b}){\rm exp}(a^{2}\kappa^{2}_{b}/\pi){\rm
erfc}(a\kappa_{b}/\sqrt{\pi})$ [42, 41, 43, 44]. The right hand side of Eq. 7
reflects the correction from ion correlation, which vanishes at the mean-field
level leading to $\Gamma=1/(z_{+}+1)>0$. It shows that overcharging cannot
happen in the mean-field electrostatics.

For $\rho_{b}>\rho^{**}_{b}$, the effective charge on the brushes change sign,
i.e. $\Gamma<0$. Counterions over compensate the bare backbone charges on
polymers. $f_{\rm ion}$ turns from repulsive to attractive in this regime,
which leads to the transition from swollen to collapsed brushes. However,
overcharging will not continue increasing; instead, it gets suppressed at high
salt concentrations because the strength of ion correlations in bulk also
increase, which reduces the driving-force for counterions to migrate to brush
region from bulk. At the critical concentration $\rho^{\rm{valley}}_{b}$ where
$\partial\Gamma/\partial\rho_{b}=0$, overcharging as well as the impact of ion
correlation gets maximized. PE brushes are compacted to the lowest height,
which is followed by a reexpansion as $\rho_{b}$ further increases. At the
high salt concentration $\rho^{***}_{b}$, the screening length reduces to a
monomer size, i.e. $\kappa_{b}b\approx 1$. The charges on brushes can be
negligible, and the brushes behave neutral again. Using the scaling theory, we
obtain the brush height in the four regimes as

$\frac{h}{Nb}\approx\begin{cases}1&(\rho_{b}<\rho^{*}_{b})\\\
{\left(\sigma\Gamma/z_{+}b\right)}^{1/3}\rho_{b}^{-1/3}&(\rho^{*}_{b}<\rho_{b}<\rho^{**}_{b})\\\
-(\sigma w_{3}z_{+}/\Gamma
b)\rho_{b}&(\rho^{**}_{b}<\rho_{b}<\rho^{***}_{b})\\\ (\sigma
w_{2}/b)^{1/3}&(\rho_{b}>\rho^{***}_{b})\\\ \end{cases}$ (2)

Figure 3: The effective relative charge on the brushes $\Gamma$ and the
surface electrostatic potential $\psi_{S}$ (normalized by $\sigma Nv/h$) vs.
$\rho_{b}$ for cations with different valency.

The relationship between chain reexpansion and charge inversion remains under
debate [18]. The latter in practice is often manifested as the reversal of the
surface electrostatic potential $\psi_{S}$ in electroosmotic flow or
electrophoritic mobility [50, 51, 44]. In the scaling theory, $\psi_{S}$ can
be approximated by the Donnan potential in the brushes:

$\displaystyle\psi_{S}$ $\displaystyle\approx-\frac{\sigma
N}{z^{2}_{+}(z_{+}+1)h\rho_{b}}\left[(z^{2}_{+}-1)\Gamma+1\right]$ (9)

It indicates the absence of charge inversion for $z_{+}=1$. For multivalent
cation, the charge inversion occurs at $\Gamma=-1/(z^{2}_{+}-1)$, which is
different from the criteria of either the collapse point $\rho^{**}_{b}$ or
the reexpansion point $\rho^{\rm{valley}}_{b}$. This is clearly illustrated by
the salt concentration dependence of $\Gamma$ and $\psi_{S}$ in Fig. 3. For
$z_{+}=2$, $\Gamma$ changes sign whereas $\psi_{S}$ remains negative: the
brush collapse-reexpansion occurs without charge inversion. For $z_{+}=3$,
both $\Gamma$ and $\psi_{S}$ change sign; however, they are out of sync and
happen at different $\rho_{b}$. There are hence no causal relationship between
the brush collapse-reexpansion and charge inversion, although they are both
induced by ion correlation. Our theoretical finding is consistent with the
simulation by Hsiao and Luijten on the conformation and mobility of a single
PE [18].

Ion correlations can have more profound impacts on the brush morphology
besides the apparent height change. At a moderate grafting density
$\sigma=0.03$nm-2, our SCFT calculation shows that brushes in the presence of
trivalent ions undergo microphase-separation and form lateral inhomogeneity.
As shown in Figs. 4a and 4b, nearby brushes are pinned to micelles due to
correlation-induced attraction, which further assemble into hexagonal pattern.
We found that the cylindrical bundles proposed by previous scaling theory is
unstable [52]. Our theoretical predictions are in good agreement with the
observations from AFM images and simulations [19, 21, 20, 53].

The micro-phase separated morphologies of PE brushes are versatile. At high
grafting densities (e.g. $\sigma=0.1$nm-2), lateral inhomogeneities disappear;
the segregation can be induced in the normal direction instead, as shown in
Figs. 4c and 4d. Oscillatory structure appears inside the collapsed layer if
ion correlations are enhanced by reducing the radius of the trivalent cation
from $2.5\mathring{\rm A}$ to $1.5\mathring{\rm A}$. Similar ion-correlation
induced oscillation has also been reported in the trivalent salt solutions
near a charged surface [44].

Figure 4: Microphase separation in the morhphology of PE brushes triggered by
trivalent ions. Lateral inhomogeneity of pinned micelles at moderate grafting
density $\sigma=0.03$nm-2 presented by a) 3D isosurface plot and b) 2D
visualization of the polymer density in the $xy$ plane at $z=0$. c) and d)
plot oscillatory layers formed in the normal direction at high grafting
density $\sigma=0.1$nm-2 as cation radius $a_{+}$ decreases from
$2.5\mathring{\rm A}$ to $1.5\mathring{\rm A}$. $\rho_{b}=1$mM.

In this letter, we develop a new theory which systematically incorporates
electrostatic fluctuations into the SCFT for polymers. Compared to existing
work, our theory fully captures the coupling between ion correlations and
polymer conformation without any preassumption on the chain structure, and is
particularly powerful in modeling inhomogeneous ion-containing polymers where
ionic strength or dielectric permittivity is spatially varying. Applied to PE
brushes, the theory predicts that ion correlations induce a non-monotonic
change of the brush height: collapse followed by reexpansion. Our theoretical
predictions are in quantitative agreement with the experimental data reported
in the literature. We perform the scaling analysis to elucidate the essential
role of the electrostatics as the competition between the repulsive osmotic
pressure due to translational entropy and the attraction induced by ion
correlations. Four regimes are identified in the height response: osmotic,
salted, collapsed and neutral. We clarify the absence of causal relationship
between the brush collapse-reexpansion and the charge inversion practically
defined as the reversal of the surface electrostatic potential. Strong ion
correlations also induce microphase separation in PE brushes. Depending on the
parameters, the structural inhomogeneities can be either in the lateral
direction as pinned micelles or in the normal direction as oscillatory layers.
The easy access of our theory to the high dimensional calculation enables us
to explore the rich morphological behaviors of PE brushes and other ion-
containing polymer systems [34, 35, 36, 37, 38, 15, 16]. Furthermore, it is
important to note that non-monotonic behaviors triggered by multivalent ions
have also been reported in the reversal of electro-osmotic flow and the
attractive force between like-charged surfaces [51, 50, 42, 43, 44, 54], which
reveals the universality of the salt concentration dependence of the
electrostatic correlation.

Acknowledgment is made to the donors of the American Chemical Society
Petroleum Research Fund for partial support of this research. This research
used the computational resources provided by the Kenneth S. Pitzer Center for
Theoretical Chemistry.

## References

  * Decher [1997] G. Decher, Fuzzy nanoassemblies: Toward layered polymeric multicomposites, Science 277, 1232 (1997).
  * Stradner  _et al._ [2004] A. Stradner, H. Sedgwick, F. Cardinaux, W. C. K. Poon, S. U. Egelhaaf, and P. Schurtenberger, Equilibrium cluster formation in concentrated protein solutions and colloids, Nature 432, 492 (2004).
  * Kim  _et al._ [2020] H. J. Kim, B. Chen, Z. Suo, and R. C. Hayward, Ionoelastomer junctions between polymer networks of fixed anions and cations, Science 367, 773 (2020).
  * Galvanetto  _et al._ [2023] N. Galvanetto, M. T. Ivanović, A. Chowdhury, A. Sottini, M. F. Nüesch, D. Nettels, R. B. Best, and B. Schuler, Extreme dynamics in a biomolecular condensate, Nature 619, 876 (2023).
  * Xiong  _et al._ [2023] T. Xiong, C. Li, X. He, B. Xie, J. Zong, Y. Jiang, W. Ma, F. Wu, J. Fei, P. Yu, and L. Mao, Neuromorphic functions with a polyelectrolyte-confined fluidic memristor, Science 379, 156 (2023).
  * Xin and Hao [2010] B. Xin and J. Hao, Reversibly switchable wettability, Chem. Soc. Rev. 39, 769 (2010).
  * Stuart  _et al._ [2010] M. A. C. Stuart, W. T. S. Huck, J. Genzer, M. Müller, C. Ober, M. Stamm, G. B. Sukhorukov, I. Szleifer, V. V. Tsukruk, M. Urban, F. Winnik, S. Zauscher, I. Luzinov, and S. Minko, Emerging applications of stimuli-responsive polymer materials, Nat. Mater. 9, 101 (2010).
  * Shabani  _et al._ [2024] S. Shabani, S. Hadjigol, W. Li, Z. Si, D. Pranantyo, M. B. Chan-Park, N. M. O’Brien-Simpson, and G. G. Qiao, Synthetic peptide branched polymers for antibacterial and biomedical applications, Nat. Rev. Bioeng. 2, 343 (2024).
  * Raviv  _et al._ [2003] U. Raviv, S. Giasson, N. Kampf, J.-F. Gohy, R. Jérôme, and J. Klein, Lubrication by charged polymers, Nature 425, 163 (2003).
  * Chen  _et al._ [2009] M. Chen, W. H. Briscoe, S. P. Armes, and J. Klein, Lubrication at physiological pressures by polyzwitterionic brushes, Science 323, 1698 (2009).
  * Rühe  _et al._ [2004] J. Rühe, M. Ballauff, M. Biesalski, P. Dziezok, F. Gröhn, D. Johannsmann, N. Houbenov, N. Hugenberg, R. Konradi, S. Minko, M. Motornov, R. R. Netz, M. Schmidt, C. Seidel, M. Stamm, T. Stephan, D. Usov, and H. Zhang, Polyelectrolyte brushes, in  _Advances in Polymer Science_ (Springer Berlin Heidelberg, 2004) pp. 79–150.
  * Romet-Lemonne  _et al._ [2004] G. Romet-Lemonne, J. Daillant, P. Guenoun, J. Yang, and J. W. Mays, Thickness and density profiles of polyelectrolyte brushes: Dependence on grafting density and salt concentration, Phys. Rev. Lett. 93, 148301 (2004).
  * Messina [2009] R. Messina, Electrostatics in soft matter, J. Phys. Condens. Matter 21, 113102 (2009).
  * Dean  _et al._ [2014] D. Dean, J. Dobnikar, A. Naji, and R. Podgornik,  _Electrostatics of Soft and Disordered Matter_ (Pan Stanford, 2014).
  * Mei  _et al._ [2006] Y. Mei, K. Lauterbach, M. Hoffmann, O. V. Borisov, M. Ballauff, and A. Jusufi, Collapse of spherical polyelectrolyte brushes in the presence of multivalent counterions, Phys. Rev. Lett. 97, 158301 (2006).
  * Kegler  _et al._ [2008] K. Kegler, M. Konieczny, G. Dominguez-Espinosa, C. Gutsche, M. Salomo, F. Kremer, and C. N. Likos, Polyelectrolyte-compression forces between spherical dna brushes, Phys. Rev. Lett. 100, 118302 (2008).
  * Yu  _et al._ [2016] J. Yu, J. Mao, G. Yuan, S. Satija, Z. Jiang, W. Chen, and M. Tirrell, Structure of polyelectrolyte brushes in the presence of multivalent counterions, Macromolecules 49, 5609 (2016).
  * Hsiao and Luijten [2006] P.-Y. Hsiao and E. Luijten, Salt-induced collapse and reexpansion of highly charged flexible polyelectrolytes, Phys. Rev. Lett. 97, 148301 (2006).
  * Yu  _et al._ [2017] J. Yu, N. E. Jackson, X. Xu, B. K. Brettmann, M. Ruths, J. J. de Pablo, and M. Tirrell, Multivalent ions induce lateral structural inhomogeneities in polyelectrolyte brushes, Sci. Adv. 3, eaao1497 (2017).
  * Jackson  _et al._ [2017] N. E. Jackson, B. K. Brettmann, V. Vishwanath, M. Tirrell, and J. J. de Pablo, Comparing solvophobic and multivalent induced collapse in polyelectrolyte brushes, ACS Macro Lett. 6, 155 (2017).
  * Yu  _et al._ [2018] J. Yu, N. E. Jackson, X. Xu, Y. Morgenstern, Y. Kaufman, M. Ruths, J. J. de Pablo, and M. Tirrell, Multivalent counterions diminish the lubricity of polyelectrolyte brushes, Science 360, 1434 (2018).
  * Luo  _et al._ [2006] G. Luo, S. Malkova, J. Yoon, D. G. Schultz, B. Lin, M. Meron, I. Benjamin, P. Vanýsek, and M. L. Schlossman, Ion distributions near a liquid-liquid interface, Science 311, 216 (2006).
  * Bakhshandeh  _et al._ [2011] A. Bakhshandeh, A. P. dos Santos, and Y. Levin, Weak and strong coupling theories for polarizable colloids and nanoparticles, Phys. Rev. Lett. 107, 107801 (2011).
  * Gupta  _et al._ [2020] A. Gupta, A. Govind Rajan, E. A. Carter, and H. A. Stone, Ionic layering and overcharging in electrical double layers in a poisson-boltzmann model, Phys. Rev. Lett. 125, 188004 (2020).
  * Muthukumar [2004] M. Muthukumar, Theory of counter-ion condensation on flexible polyelectrolytes: Adsorption mechanism, J. Chem. Phys. 120, 9343 (2004).
  * Kundagrami and Muthukumar [2008] A. Kundagrami and M. Muthukumar, Theory of competitive counterion adsorption on flexible polyelectrolytes: Divalent salts, J. Chem. Phys. 128, 244901 (2008).
  * Li  _et al._ [2022] M. Li, B. Zhuang, and J. Yu, Effects of ion valency on polyelectrolyte brushes: A unified theory, Macromolecules 55, 10450 (2022).
  * Li  _et al._ [2024] M. Li, X. Xu, S. Yasar, H.-A. Klok, B. Zhuang, and J. Yu, Ion adsorption enhances apparent nonelectrostatic attraction between monomers in polyelectrolyte brushes, Macromolecules 57, 3026 (2024).
  * Shen and Wang [2017] K. Shen and Z.-G. Wang, Electrostatic correlations and the polyelectrolyte self energy, J. Chem. Phys. 146, 084901 (2017).
  * Shen and Wang [2018] K. Shen and Z.-G. Wang, Polyelectrolyte chain structure and solution phase behavior, Macromolecules 51, 1706 (2018).
  * Sing and Qin [2023] C. E. Sing and J. Qin, Bridging field theory and ion pairing in the modeling of polyelectrolytes and complex coacervation, Macromolecules 56, 5941 (2023).
  * Dobrynin and Rubinstein [2005] A. V. Dobrynin and M. Rubinstein, Theory of polyelectrolytes in solutions and at surfaces, Prog. Polym. Sci. 30, 1049 (2005).
  * Duan  _et al._ [2020] C. Duan, W. Li, and R. Wang, Conformation of a single polyelectrolyte in poor solvents, J. Chem. Phys. 153, 064901 (2020).
  * Sing  _et al._ [2013] C. E. Sing, J. W. Zwanikken, and M. O. de la Cruz, Interfacial behavior in polyelectrolyte blends: Hybrid liquid-state integral equation and self-consistent field theory study, Phys. Rev. Lett. 111, 168303 (2013).
  * Sing  _et al._ [2014] C. E. Sing, J. W. Zwanikken, and M. O. de la Cruz, Electrostatic control of block copolymer morphology, Nat. Mater. 13, 694 (2014).
  * Gurovitch and Sens [1999] E. Gurovitch and P. Sens, Adsorption of polyelectrolyte onto a colloid of opposite charge, Phys. Rev. Lett. 82, 339 (1999).
  * Dobrynin  _et al._ [2000] A. V. Dobrynin, A. Deshkovski, and M. Rubinstein, Adsorption of polyelectrolytes at an oppositely charged surface, Phys. Rev. Lett. 84, 3101 (2000).
  * Winkler and Cherstvy [2006] R. G. Winkler and A. G. Cherstvy, Critical adsorption of polyelectrolytes onto charged spherical colloids, Phys. Rev. Lett. 96, 066103 (2006).
  * Li and Wu [2006] Z. Li and J. Wu, Density functional theory for polyelectrolytes near oppositely charged surfaces, Phys. Rev. Lett. 96, 048302 (2006).
  * Wang [2010] Z.-G. Wang, Fluctuation in electrolyte solutions: The self energy, Phys. Rev. E 81, 021501 (2010).
  * Agrawal and Wang [2022a] N. R. Agrawal and R. Wang, Self-consistent description of vapor-liquid interface in ionic fluids, Phys. Rev. Lett. 129, 228001 (2022a).
  * Agrawal and Wang [2022b] N. R. Agrawal and R. Wang, Electrostatic correlation induced ion condensation and charge inversion in multivalent electrolytes, J. Chem. Theory Comput. 18, 6271 (2022b).
  * Agrawal  _et al._ [2023] N. R. Agrawal, R. Kaur, C. Carraro, and R. Wang, Ion correlation-driven like-charge attraction in multivalent salt solutions, J. Chem. Phys. 159, 244905 (2023).
  * Agrawal  _et al._ [2024] N. R. Agrawal, C. Duan, and R. Wang, Nature of overcharging and charge inversion in electrical double layers, J. Phys. Chem. B 128, 303 (2024).
  * Fredrickson [2006] G. H. Fredrickson,  _The Equilibrium Theory of Inhomogeneous Polymers_, International Series of Monographs on Physics (OUP Oxford, 2006).
  * Duan and Wang [2023] C. Duan and R. Wang, Electrostatics-induced nucleated conformational transition of protein aggregation, Phys. Rev. Lett. 130, 158401 (2023).
  * Frigo [1999] M. Frigo,  _A Fast Fourier Transform Compiler_ , Proceedings of the 1999 ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI ’99) (Atlanta, Georgia, 1999).
  * Hoffman and Frankel [2018] J. D. Hoffman and S. Frankel,  _Numerical Methods for Engineers and Scientists_ (CRC Press, 2018).
  * Zhulina  _et al._ [1995] E. B. Zhulina, T. M. Birshtein, and O. V. Borisov, Theory of ionizable polymer brushes, Macromolecules 28, 1491 (1995).
  * Besteman  _et al._ [2004] K. Besteman, M. A. G. Zevenbergen, H. A. Heering, and S. G. Lemay, Direct observation of charge inversion by multivalent ions as a universal electrostatic phenomenon, Phys. Rev. Lett. 93, 170802 (2004).
  * van der Heyden  _et al._ [2006] F. H. J. van der Heyden, D. Stein, K. Besteman, S. G. Lemay, and C. Dekker, Charge inversion at high ionic strength studied by streaming currents, Phys. Rev. Lett. 96, 224502 (2006).
  * Brettmann  _et al._ [2017] B. Brettmann, P. Pincus, and M. Tirrell, Lateral structure formation in polyelectrolyte brushes induced by multivalent ions, Macromolecules 50, 1225 (2017).
  * Tan  _et al._ [2023] H.-G. Tan, L.-D. Guo, H.-Y. Qian, L.-X. Liu, Q.-H. Hao, and B. Miao, Size effect of multivalent counterions on polyelectrolyte brushes in different polar solvents, Macromolecules 56, 9312 (2023).
  * Angelini  _et al._ [2003] T. E. Angelini, H. Liang, W. Wriggers, and G. C. L. Wong, Like-charge attraction between polyelectrolytes induced by counterion charge density waves, Proc. Natl. Acad. Sci. U.S.A. 100, 8634 (2003).