# Theoretical and computational analysis of the electrophoretic polymer
mobility inversion induced by charge correlations

Xiang Yang  Department of Applied Physics, Aalto University, P.O. Box 11000,
FI-00076 Aalto, Finland  Sahin Buyukdagli  Department of Physics, Bilkent
University, Ankara 06800, Turkey  Alberto Scacchi  Department of Applied
Physics, Aalto University, P.O. Box 11000, FI-00076 Aalto, Finland.  Academy
of Finland Center of Excellence in Life-Inspired Hybrid Materials (LIBER),
Aalto University, P.O. Box 16100, FI-00076 Aalto, Finland  Maria Sammalkorpi
Academy of Finland Center of Excellence in Life-Inspired Hybrid Materials
(LIBER), Aalto University, P.O. Box 16100, FI-00076 Aalto, Finland  Department
of Chemistry and Materials Science, Aalto University, P.O. Box 16100, FI-00076
Aalto, Finland  Department of Bioproducts and Biosystems, Aalto University,
P.O. Box 16100, FI-00076 Aalto, Finland  Tapio Ala-Nissila  Department of
Applied Physics, Aalto University, P.O. Box 11000, FI-00076 Aalto, Finland.
Quantum Technology Finland Center of Excellence, Department of Applied
Physics, Aalto University, P.O. Box 11000, FI-00076 Aalto, Finland
Interdisciplinary Centre for Mathematical Modelling and Department of
Mathematical Sciences, Loughborough University, Loughborough, Leicestershire
LE11 3TU, UK

(November 17, 2022)

###### Abstract

Electrophoretic (EP) mobility reversal is commonly observed for strongly
charged macromolecules in multivalent salt solutions. This curious effect
takes place, e.g., when a charged polymer, such as DNA, adsorbs excess
counterions so that the counterion-dressed surface charge reverses its sign,
leading to the inversion of the polymer drift driven by an external electric
field. In order to characterize this seemingly counterintuitive phenomenon
that cannot be captured by electrostatic mean-field theories, we adapt here a
previously developed strong-coupling-dressed Poisson-Boltzmann approach to the
cylindrical geometry of the polyelectrolyte-salt system. Within the framework
of this formalism, we derive an analytical polymer mobility formula dressed by
charge correlations. In qualitative agreement with polymer transport
experiments, this mobility formula predicts that the increment of the
monovalent salt, the decrease of the multivalent counterion valency, and the
increase of the dielectric permittivity of the background solvent, suppress
charge correlations and increase the multivalent bulk counterion concentration
required for EP mobility reversal. These results are corroborated by coarse-
grained molecular dynamics simulations showing how multivalent counterions
induce mobility inversion at dilute concentrations and suppress the inversion
effect at large concentrations. This re-entrant behavior, previously observed
in the aggregation of like-charged polymer solutions, calls for verification
by polymer transport experiments.

††preprint: APS/123-QED

##  I Introduction

Electrostatic correlation effects are ubiquitous in biological systems
involving strongly charged biomolecules and membranes. Counterintuitive
phenomena, such as like-charge attraction, are typically observed in systems
including macromolecules in contact with multivalent counterions Grosberg  _et
al._ (2002); Levin (2002); Angelini  _et al._ (2003); Besteman  _et al._
(2007, 2004); Butler  _et al._ (2003). Several mechanisms mediating the
attraction between like-charged rods have been considered in literature, such
as, covalence-like binding Ray and Manning (1997), Gaussian-fluctuation
correlations Podgornik and Parsegian (1998), and structural correlations
Kornyshev and Leikin (1997, 1999). In dense polymer systems, like-charge
attraction can induce, e.g., bundle formation of F-actin and toroidal
aggregates of DNA Tang  _et al._ (1996); Bloomfield (1997). Multivalent
inorganic ions and polyamines can also act as condensing agents Teif and
Bohinc (2011).

Charge inversion (CI) is another interesting manifestation of charge
correlations. This phenomenon occurs when the macromolecular surface charge
flips its sign upon the adsorption of a sufficient amount of counterions from
the solution. When combined with other physical factors, CI can induce
additional effects. For example, in a press-driven flow through a negatively
charged slit pore, multivalent cation addition can invert the sign of the
monovalent counterion current, generating a like-charged streaming current of
negative sign Van der Heyden  _et al._ (2006). One should also mention the
important role played by CI in biology. In the cell medium of eukaryotic
organisms, the negatively charged DNA and the positively charged histone
proteins can assemble into nucleosome. At high salt concentrations, the stable
structure corresponds to a rope-like DNA wrapped around a bead-like histone
having an inverted net charge Rippe  _et al._ (2008); Grosberg  _et al._
(2002).

Electrophoresis is the motion of dispersed charged particles relative to a
fluid exposed to an electric field. The resulting driven transport is
quantified in terms of the electrophoretic (EP) mobility, defined as the ratio
of the drift velocity and the electric field strength. Based on the EP
mobility strength, one can separate macromolecules in terms of their size and
surface charge Kaper  _et al._ (2003); Danger  _et al._ (2007). Furthermore,
EP transport experiments can be efficiently used to probe the interfacial
charge structure of polyelectrolytes Kabanov and Kabanov (1995); Kabanov  _et
al._ (1996). Additionally, from viral infection to nanopore-based polymer
sensing, polymer transport through confined pores plays a critical role in
biological processes and various nanoscale applications Holm  _et al._ (2001);
Schoch  _et al._ (2008); Jain  _et al._ (2016); Thomas  _et al._ (2016). The
underlying electrohydrodynamic mechanisms have been extensively explored in
recent theoretical works Buyukdagli (2020, 2022).

The surface CI of polyelectrolytes has been previously characterized by EP
transport experiments Wang  _et al._ (2016, 2018). Indeed, these experiments
have shown that the addition of multivalent cations into a polyelectrolyte
solution induces the EP motion of the negatively charged polyelectrolytes
along the external electric field. In this article, we carry out a theoretical
investigation of this seemingly counterintuitive mobility reversal effect. As
the occurrence of this phenomenon requires the presence of multivalent
counterions strongly coupled to the polyelectrolyte charges, the
characterization of the underlying mechanism necessitates the use of a
correlation-corrected electrostatic framework. Thus, considering the linear
response limit of a previously developed strong coupling (SC)-dressed Poisson-
Boltzmann (PB) approach Kanduč  _et al._ (2010); Kanduc  _et al._ (2011);
Buyukdagli (2020, 2022), we derive an analytical polymer transport formula
accounting for the multivalent ion-induced charge correlations.

We show that our correlation-corrected EP transport formula can qualitatively
reproduce various experimentally observed correlation effects on polymer
transport, such as the onset of the mobility reversal by added multivalent
counterions, and its weakening by the increment of the monovalent salt
component, the increase of the solvent permittivity, and the reduction of the
multivalent counterion valency. The analytical structure of our formalism
enables a clear interpretation of the correlation mechanism driving these
effects Wang  _et al._ (2016, 2018).

As a complementary approach, we also perform particle-based coarse-grained
molecular dynamics (MD) simulations. A benchmark between MD simulations and
the dressed-ion theory is crucial, since the latter has been previously tested
for solutions in contact with charged planes Kanduč  _et al._ (2010); Kanduc
_et al._ (2011); Buyukdagli (2020, 2022), but not in the case of charged
cylinders. In addition to reproducing consistently the aforementioned
theoretical predictions, in the dense multivalent counterion regime, our MD
simulations capture a re-entrant phase not covered by the theoretical mobility
formula, which is valid only at dilute counterion densities. Finally, within
our computational framework, we also characterize the effect of the counterion
charge distribution beyond the point-ion approximation. The limitations of our
theory and potential improvements are discussed in the Summary and Conclusions
section.

##  II Molecular Dynamics simulation approach

###  II.1 Simulation system

We start by discussing the details of the MD simulations in this work. The
simulations were run in a box of volume $L_{x}\times L_{y}\times L_{z}$, where
the longitudinal dimension was set to $L_{z}=20$ nm. Depending on the added
salt concentration, the transverse dimensions $L_{x}=L_{y}$ were altered
between 24 nm and 240 nm. As explained in the Supplementary Material (SM), we
chose $L_{z}$ to obviate finite-size effects from the data.

As depicted in Fig. 1, the system includes a coarse-grained (CG) DNA chain
centered at the coordinates $(x,y)=(0,0)$ spanning the cuboid simulation box
along the $z$ axis. The DNA molecule was modeled by a series of spherical
beads distributed uniformly along the polymer axis and separated by $b=0.17$
nm from each other. The value of $b$ is small enough to provide a smooth
potential surface in the longitudinal direction. Moreover, each bead carries
an elementary charge $-e$, resulting in a linear charge density of
$\lambda_{0}=-e/b=-5.9\;e$/nm. The radius of the polymer, or equivalently the
radius of the beads, was set to $a=1.2$ nm, which corresponds to the
characteristic thickness of double-stranded DNA molecules Forrey and
Muthukumar (2006); Mandelkern  _et al._ (1981).

Here the simulations were in the canonical ensemble, i.e. with fixed number of
particles. Therefore, electroneutrality in our simulations was achieved by
neutralizing the DNA charges with $L_{z}/b$ added monovalent counterions (e.g.
Na+). This lead to a counterion concentration of $(L_{x}L_{y}b)^{-1}$.
Additional monovalent salt, such as NaCl, and multivalent salt species of
general chemical structure XCl${}_{q_{\rm c}}$, were added into the solution.
The different cations (X) of equal size and valency $q_{\rm c}$ were divalent,
trivalent (T3+), quadrivalent or octavalent charges. Within the theory
considerations, the bulk concentration of the added monovalent and multivalent
ions will be denoted by $n_{\rm sb}$ and $n_{\rm cb}$, respectively.

Figure 1: Left panel: A schematic of the CG DNA located at $(x,y)=(0,0)$ along
the $z$ axis of the simulation box with total volume $L_{x}\times L_{y}\times
L_{z}$. Right panel: Ions distributed within a cylindrical volume of radius
$r$ centered at $(x,y)=(0,0)$. The red, green and orange spheres represent
trivalent cations, monovalent anions and monovalent cations, respectively. In
this snapshot, $L_{x}=L_{y}=24$ nm, and $L_{z}=20$ nm. The concentration of
multivalent salt (TCl3) is $n_{\rm cb}=200$ mM and that of monovalent
counterions (e.g. Na+) neutralizing the CG DNA charge 17 mM.

Overlapping of the mobile charges was avoided by placing the ions on a regular
grid at the initialization stage. The pairwise interactions between the
polymer beads and the ions were modelled via the standard Weeks-Chandlers-
Andersen potential Weeks  _et al._ (1971)

$V_{ij}(r)=\left\\{4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r}\right)^{12}-\left(\frac{\sigma_{ij}}{r}\right)^{6}\right]+\epsilon_{ij}\right\\}\theta(r_{ij}^{\rm
cut}-r),$ (1)

where the indices $i$ and $j$ correspond to the multivalent, sodium, chloride
ions, and DNA beads; $\theta(x)$ is the Heaviside step function. In the
implementation of the potential (1), we used the Lorentz-Berthelot mixing
rules $\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$ and
$\epsilon_{ij}=\sqrt{\epsilon_{i}\epsilon_{j}}$, and set the cut-off radii of
the pairwise interactions $r_{ij}^{\rm cut}=2^{1/6}\sigma_{ij}$. The
parameters $\epsilon_{i}$ and $\sigma_{i}$ are reported in Table 1. The
approach follows our previous work, see Ref. Vahid  _et al._ (2022) for more
details.

The pairwise electrostatic interactions between charges $q_{i}$ and $q_{j}$,
with respective position vectors $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$, were
taken into account with the Coulomb potential

$V^{\rm
c}(\mathbf{r}_{i}-\mathbf{r}_{j})=\frac{q_{i}q_{j}e^{2}}{4\pi\varepsilon_{\rm
w}||\mathbf{r}_{i}-\mathbf{r}_{j}||}.$ (2)

These potentials treat the solution as a continuous dielectric medium with
permittivity $\varepsilon_{\rm w}$. Unless stated otherwise, the relative
dielectric constant of the Coulomb potential was set to the relative
permittivity of water $\varepsilon_{\rm w}/\varepsilon_{0}=78$, where
$\varepsilon_{0}$ stands for the vacuum permittivity. The MD simulations were
performed using the LAMMPS Jan2020 package Plimpton (1995); Thompson  _et al._
(2022).

The long-range electrostatic interactions were calculated with the particle-
particle particle-mesh (PPPM) method Hockney and Eastwood (2021). Up to the
characteristic split distance $d^{\rm s}$, the pairwise Coulomb interactions
in Eq. (2) were evaluated in real space, whereas beyond this distance in
reciprocal space. The PPPM accuracy parameter (relative error) was set to
$10^{-5}$ and the stencil size parameter to 5. The setting of the distance
$d^{\rm s}$ is discussed in the SM. All simulations were performed in the
$NVT$ ensemble. During the simulations, the temperature was controlled by the
Nose-Hoover thermostat with $300$ K as the reference temperature Nosé (1984);
Hoover (1985). After an initial placement of the ions into the simulation box,
the system energy was minimized by the conjugate gradient method. This was
followed by a 24 ns $NVT$ simulation run for data analysis. The first 4 ns of
this run were disregarded. A 2 fs time step was used.

Variable | DNA bead | Na+ | Cl- | multivalent ion
---|---|---|---|---
$\epsilon$ (kcal/mol) | 0.1 | 0.13 | 0.124 | 0.1
$\sigma$ (nm) | 2.4 | 0.234 | 0.378 | 0.5
Table 1: The parameters of the Weeks-Chandlers-Andersen potential used in the
CG model.

###  II.2 DNA Charge Inversion

In our DNA-liquid system, the mobile charge configuration is set by the
collective effect of the electrostatic and steric interactions between the
ions and the DNA beads, and thermal fluctuations suppress the electrostatic
coupling of the charged entities. It has been previously shown that in the
presence of a sufficient amount of multivalent counterions in the liquid, the
electrostatic many-body effects take over the thermal fluctuations, leading to
the overcompensation of the bare DNA charges by the counterions bound to the
DNA molecule Netz (2001); Wang  _et al._ (2016, 2018); Hsiao (2008); Kanduc
_et al._ (2011); Buyukdagli (2022); Aksimentiev (2010).

Figure 2: Radial functions related to charge inversion from the MD
simulations. (a) Normalized density profiles of the T3+ cations (solid curves)
and Cl- anions (dashed curves). (b) Cumulative charge density in Eq. (3). (c)
Electrostatic potential in Eq. (4) at various trivalent salt concentrations
$n_{\rm cb}$ indicated in the legend of (b). The inset in (a) displays the
dimensionless Cl- densities on a smaller linear scale. The dashed vertical
line in (c) marks the no-slip boundary located at $r=a+\delta$. The gray
regions indicate the radius of the polymer $a$. The relative liquid
permittivity is $\varepsilon_{\rm w}/\varepsilon_{0}=78$.

With the aim of illustrating the corresponding CI effect, and of relating the
latter to the EP mobility inversion of the molecule, we focus first on the
configuration of trivalent salt (TCl3). Figure 2(a) displays the dimensionless
density profiles $n_{i}(r)/n_{i\rm{b}}$ of the multivalent cations (solid
curves) and monovalent anions (dashed curves) at various bulk TCl3
concentrations from the MD simulations. The density profiles $n_{i}(r)$ are
calculated from counting the ions in a cylindrical volume around the DNA (cf.
Fig. 1). One can see that the T3+ adsorption peaks, located at $r\approx 1.4$
nm, are followed by the Cl- density peaks at $r\approx 1.8$ nm. Moreover, the
inset shows that these peaks correspond to the interfacial excess of Cl- ions.
Thus, the Cl- attraction by the adsorbed T3+ cations leads to an apparent
like-charge binding of the Cl- anions onto the anionic DNA molecule.

The DNA zeta potential setting the EP mobility of the polymer can be obtained
from the cumulative charge density $\lambda_{t}(r)$ of the molecule. The
latter corresponds to the total charge enclosed by a cylindrical volume of
radius $r$ divided by the length of the cylinder, $L_{z}$ (see Fig. 1), i.e.

$\lambda_{t}(r)=\lambda_{0}+2\pi
e\int_{0}^{r}dr^{\prime}r^{\prime}[n_{\rm+}(r^{\prime})-n_{\rm-}(r^{\prime})+q_{\rm
c}n_{\rm c}(r^{\prime})],$ (3)

where $n_{+}(r)$, $n_{-}(r)$ and $n_{c}(r)$ denote the number density of the
monovalent Na+ and Cl- ions, and the multivalent X${}^{q_{\rm c}+}$ cations,
respectively. Figure 2(b), displaying the cumulative charge density in Eq.
(3), shows that CI occurs in the bulk T3+ concentration regime $n_{\rm
cb}\gtrsim 20$ mM.

Figure 3: EP mobility $\mu_{\rm p}$ from MD simulations as a function of the
multivalent counterion concentration $n_{\rm cb}$ for various values of (a)
the cation valency $q_{\rm c}$, (b) the dielectric constant, and (c) the
monovalent salt concentration $n_{\rm sb}$. The dashed curves correspond to a
quadratic fitting function (see main text), and the crosses represent the
reversal concentrations $n_{\rm cb}^{*}$, which are reported in Table 2.

By integrating the radial Poisson equation, one can relate the radial
component of the electric field $E(r)$ to the cumulative charge density as
$E(r)=\lambda_{t}(r)/2\pi r\varepsilon_{\rm w}$. Upon integration of the
latter equality, the electrostatic potential $\phi(r)$ takes the form

$\phi(r)=\int_{r}^{\infty}\frac{\lambda_{t}(r^{\prime})}{2\pi\varepsilon_{\rm
w}r^{\prime}}dr^{\prime}.$ (4)

The resulting potential profiles displayed in Fig. 2(c) show that CI manifests
itself as the emergence of the electrostatic potential peaks, where the slope
of $\phi$ switches from positive to negative.

As described in Sec. III, the coupled solutions of the electrostatic Poisson
equation and the hydrodynamic Stokes relation yield the Helmholtz-Smoluchowski
identity, which relates the EP polymer mobility to the DNA zeta potential
$\zeta$ via

$\mu_{\rm p}=\frac{\varepsilon_{\rm w}\zeta}{\eta},$ (5)

where $\eta=8.91\times 10^{-4}$ Pa s is the dynamic viscosity of water. In Eq.
(5), the zeta potential corresponds to the electrostatic potential value at
the no-slip surface separating the mobile ions from the ones bound to the
charged polymer Bhattacharjee (2016). Given this definition, the zeta
potential can be obtained from our potential profiles as
$\zeta=\phi(a+\delta)$. Based on the results of previous experiments Yamaguchi
and Kobayashi (2016); Galla  _et al._ (2014), in our calculations the
thickness of the no-slip region was set to $\delta=0.5$ nm.

###  II.3 Electrophoretic mobility inversion

We next focus on the EP mobility under a variety of conditions from the MD
simulations. Figure 3(a) illustrates the dependence of the EP mobility on the
multivalent cation at concentration $n_{\rm cb}$ for various valencies. First,
the plot shows that for $q_{\rm c}=2-4$ the increment of the concentration
$n_{\rm cb}$ increases the negative polymer mobility and switches it to
positive. The characteristic concentration $n_{\rm cb}^{*}$ at which the
mobility is inverted can be identified by a quadratic fit of the curve
$\mu_{\rm p}(n_{\rm cb})=0$ in the vicinity of the inversion point. In Fig. 3
these fitting functions are shown by dashed lines. The reversal concentrations
provided by the fits are $n_{\rm cb}^{*}=280$, $5$, and $0.05$ mM, for
divalent, trivalent, and quadrivalent salt, respectively. This is in agreement
with prior findings on asymmetricity of the salt and charge reversal Antila
_et al._ (2017). Additionally, in agreement with polymer transport experiments
Wang  _et al._ (2016, 2018), the increase of the counterion valency lowers the
critical salt concentration $n_{\rm cb}^{*}$ required for the occurrence of
the mobility reversal, i.e. $q_{\rm c}\uparrow n_{\rm cb}^{*}\downarrow$.

The intensification of the DNA CI and the resulting mobility reversal by ion
valency stems from the amplification of charge correlations. Indeed, in Sec.
III.3, we show that the weight of the charge correlations responsible for CI
is proportional to the electrostatic coupling parameter defined as $\Xi=q_{\rm
c}^{2}\ell_{\rm B}/\mu_{\rm GC}$, where $\ell_{\rm
B}=e^{2}/(4\pi\varepsilon_{\rm w}k_{\rm B}T)$ stands for the Bjerrum length,
and $\mu_{\rm GC}=1/(2\pi q_{\rm c}\ell_{\rm B}\sigma_{\rm p})$ is the Gouy-
Chapman (GC) length. Here $\sigma_{\rm p}$ is the surface-charge number
density of the polymer, and $\sigma_{\rm p}=1/2\pi ab$ in the simulations.
These identities imply that the coupling parameter is a cubic function of the
ion valency ($\Xi\propto q_{\rm c}^{3}$), explaining the sharp emergence of
mobility reversal upon the rise of the counterion valency at constant
concentration $n_{\rm cb}$, as can be seen in Fig. 3(a).

Figure 3(a) also shows that for $q_{\rm c}\geq 2$, the EP mobility exhibits a
non-monotonic dependence on $n_{\rm cb}$. Namely, as the multivalent cation
concentration is increased beyond the value $n_{\rm cb}^{*}$, the reversed
mobility rises, reaches a peak, and drops monotonically. This peculiarity is
qualitatively similar to the re-entrance phenomenon observed in like-charge
polymer interactions; according to prior experimental studies and theoretical
analysis of these systems, multivalent cations triggering DNA condensation at
low concentrations reverse the effect at large concentrations, resulting in
the segregation of the condensates in the charged liquid Hsiao (2008);
Buyukdagli (2017); Truzzolillo  _et al._ (2018).

The permittivity of the solvent is an additional control parameter previously
investigated by transport experiments Wang  _et al._ (2016). The effect of the
liquid permittivity on the EP polymer velocity is displayed in Fig. 3(b). One
can see that at fixed multivalent ion concentration, the reduction of the
solvent permittivity rises the negative DNA mobility and switches the latter
to positive, i.e. $\varepsilon_{\rm w}\downarrow\mu_{\rm p}\uparrow$. This
implies that the mixing of the water solvent with a lower permittivity liquid
can solely trigger DNA mobility inversion. This result is in agreement with
polymer transport experiments in Ref. 26, where the negative DNA velocity in
water-ethanol mixtures was observed to rise and reverse with increasing the
volume fraction of ethanol in the liquid.

The occurrence of mobility reversal upon the reduction of the solvent
permittivity is a consequence of stronger charge correlations in lower
permittivity liquids. Indeed, according to its definition above, the
electrostatic coupling parameter scales quadratically with the inverse
dielectric permittivity, i.e. $\Xi\propto\varepsilon_{\rm w}^{-2}$. Hence, the
reduction of the liquid permittivity enhances the weight of the electrostatic
many-body interactions responsible for CI and mobility reversal. Owing to this
mechanism, Fig. 3(b) shows that the multivalent ion concentration at the
reversal drops with the liquid permittivity, i.e. $\varepsilon_{\rm
w}\downarrow n_{\rm cb}^{*}\downarrow$.

Finally, in Fig. 3(c), we investigate the effect of added monovalent salt,
such as NaCl, at concentration $n_{\rm sb}$, on DNA mobility. The plot
indicates that monovalent salt ions counteract the multivalent cations and
monotonically suppress the mobility inversion. Consequently, in accordance
with polymer transport experiments Wang  _et al._ (2016, 2018), the minimum
counterion concentration for mobility reversal increases with the amount of
monovalent salt, i.e. $n_{\rm sb}\uparrow n^{*}_{\rm cb}\uparrow$. In Sec.
III, within the framework of our correlation-augmented EP transport theory, we
show that this feature originates from the attenuation of charge correlations
by monovalent ions.

Parameters |  $n^{*}_{\rm cb}$ [mM]
---|---
Fig. 3 (a) |  |  |  |
$q_{\rm c}=2,3,4,8$ |  n/a |  0.055 |  5 |  280
Fig. 3 (b) |  |  |  |
$\varepsilon_{\rm w}/\varepsilon_{0}=100,78,60,40$ |  n/a |  0.3 |  5 |  35
Fig. 3 (c) |  |  |  |
$n_{\rm sb}=0,6,20,50$ mM |  5 |  7.2 |  10 |  15
Fig. 4 |  |  |  |
rod A, B, C, D |  70 |  80 |  140 |  260
Table 2: The reversal concentrations $n^{*}_{\rm cb}$ from MD simulations for
different parameters in Figs. 3 and 4. n/a indicates that $n^{*}_{\rm cb}$ is
too small to be reliably estimated from the simulations.

###  II.4 Influence of charge distribution in multivalent counterions

In order to extend our understanding of the mechanism behind mobility
inversion beyond the counterion valency, we scrutinize the role played by the
spatial charge distribution in the multivalent counterions on polymer
mobility. To this end, we ran MD simulations by replacing the spherical
trivalent counterions by rod-like charges with the same valency. Each rod-like
ion of total length $L=1.6$ nm consists of seven spheres linearly distributed.
The diameter of each sphere is $0.5$ nm, same as previous multivalent
spherical ions (see Table I), and the distance between two adjacent spheres is
$0.18$ nm.

First, we analyze the purely steric effect originating from the finite length
of the trivalent ions. To this aim, we place three elementary charges on the
central sphere of the rod (rod A in Fig. 4), and compare the resulting polymer
mobility (cyan curve) to the case with spherical trivalent ions (red curve).
The plot indicates that as the rotational penalty experienced by the rod-like
ions close to the polymer surface reduces their density and their degree of
condensation, the finite counterion size decreases the polymer mobility whilst
increasing the critical counterion concentration required for mobility
inversion from $n_{\rm cb}^{*}=7$ mM to $n_{\rm cb}^{*}=70$ mM.

Second, we investigate the electrostatic effect associated with the surface
charge density of the multivalent counterions. To this purpose, we split the
three elementary charges on the ion by moving two unit charges from the center
to the end in a symmetric fashion (see the corresponding configurations in the
legend of Fig. 4). One can see that the resulting reduction of the ionic
surface charge density weakens the DNA screening by the counterions, lowering
the EP polymer mobility monotonically at all concentrations. As a result, the
critical ion concentration for mobility reversal rises from $n_{\rm
cb}^{*}=70$ mM for rod A-like counterions to $n_{\rm cb}^{*}=80$ mM, $140$ mM,
and $260$ mM for rod-like counterions of type B, C, and D, respectively.

Figure 4: EP mobility from MD simulations in the presence of trivalent
counterions of different intramolecular structures. The reversal
concentrations are found at $n_{\rm cb}^{*}=7,70,80,140$ and 260 mM for
sphere, rod A, B, C and D, respectively. The reversal concentrations $n_{\rm
cb}^{*}$ are shown in Table 2. Here the concentration of monovalent salt is
$n_{\rm sb}=0$. In the inset the red circle denotes a spherical trivalent
counterion, whereas the yellow circles are monovalent spherical counterions
located at different sites in the seven-bead rodlike molecule.

##  III SC-dressed EP transport theory

In order to physically understand the various aspects of the EP mobility
reversal investigated in the previous section by MD simulations, we develop
here a SC-corrected analytical theory of EP polymer transport. The theory is
based on the dressed-ion limit Kanduc  _et al._ (2011) of the SCPB formalism
Buyukdagli (2020, 2022), treating the monovalent salt ions at the weak-
coupling level while taking into account the SC correlations mediated by the
multivalent charge species. We emphasize that due to its grand-canonical
nature, the present theory can consistently account for the bulk charge
reservoir present in the real system. This feature leads into quantitative
differences between the theory and the canonical MD simulations here. In
contrast to the MD simulations, where the finite number of charges in the
setup required to force the electroneutrality condition via the addition of
extra monovalent cations, in our grand-canonical formalism the inclusion of
the chemical equilibrium between the interfacial and bulk charges allows to
satisfy automatically both the DNA and the bulk electroneutrality conditions.

###  III.1 Derivation of the SC-dressed polymer mobility

The EP mobility of a cylindrical polymer with radius $a$ translocating through
a nanopore of radius $d$ has been previously calculated by the coupled
solution of the Navier-Stokes and Poisson equations in Ref. Buyukdagli (2022).
Taking the limit $d\to\infty$ to remove the membrane interface from the model,
the EP polymer mobility in the bulk liquid follows from this previous result
as

$\mu_{\rm p}=\mu_{\rm ep}\left[\phi(a^{*})-\phi(d\to\infty)\right],$ (6)

where $\phi(\mathbf{r})=eV(\mathbf{r})/(k_{\rm B}T)$ is the dimensionless
electrostatic potential. Moreover, we introduced the EP mobility coefficient
defined as $\mu_{\rm ep}=\varepsilon_{\rm w}\varepsilon_{0}k_{\rm
B}T/(e\eta)$, with the dielectric permittivity of vacuum $\varepsilon_{0}$ and
water $\varepsilon_{\rm w}$, the thermal energy $k_{\rm B}T$, the electron
charge $e$, and the water viscosity $\eta$. In Eq. (6), the effective polymer
radius $a^{*}$ is related to the physical polymer radius $a$ by
$a^{*}=a+\delta$, where $\delta$ is the hydrodynamic no-slip length. We also
note that the relation (6) is naturally equivalent to Eq. (5).

The evaluation of the mobility (6) requires the calculation of the average
potential. In the dressed-ion limit Kanduc  _et al._ (2011) of the SCPB
formalism Buyukdagli (2020) describing the electrostatics of multivalent
electrolyte mixtures, the SC-dressed average electrostatic potential is given
by

$\phi(\mathbf{r})=\phi_{\rm s}(\mathbf{r})+\phi_{\rm c}(\mathbf{r}),$ (7)

where the average potential component associated with the monovalent salt and
satisfying the linear PB equation reads

$\phi_{\rm
s}(\mathbf{r})=\int\mathrm{d}^{3}\mathbf{r}^{\prime}G(\mathbf{r},\mathbf{r}^{\prime})\sigma(\mathbf{r}^{\prime}),$
(8)

and the SC potential induced by the multivalent-ion component of valency
$q_{\rm c}$ and bulk concentration $n_{\rm cb}$ is

$\phi_{\rm c}(\mathbf{r})=q_{\rm c}n_{\rm
cb}\int\mathrm{d}^{3}\mathbf{r}^{\prime}G(\mathbf{r},\mathbf{r}^{\prime})k_{\rm
c}(\mathbf{r}^{\prime}).$ (9)

In Eqs. (8)-(9), we introduced the fixed surface charge of the polymer
$\sigma(\mathbf{r})=-\sigma_{\rm p}\delta(r-a)$, and the electrostatic Green
function satisfying the kernel equation

$\left[\nabla\cdot\varepsilon(\mathbf{r})\nabla-\kappa^{2}\varepsilon(\mathbf{r})\theta(r-a)\right]G(\mathbf{r},\mathbf{r}^{\prime})=-\frac{e^{2}}{k_{\rm
B}T}\delta(\mathbf{r}-\mathbf{r}^{\prime}).$ (10)

In Eq. (10), we defined the radial dielectric permittivity profile

$\varepsilon(\mathbf{r})=\varepsilon_{\rm p}\theta(a-r)+\varepsilon_{\rm
w}\theta(r-a),$ (11)

with the water permittivity $\varepsilon_{\rm w}$ and the polymer permittivity
$\varepsilon_{\rm p}$. Moreover, the salt screening parameter is defined as
$\kappa^{2}=4\pi\ell_{\rm B}\left(n_{+{\rm b}}+n_{-{\rm b}}\right)$, where
$n_{\pm{\rm b}}$ stands for the bulk concentration of the monovalent salt
ions. As discussed at the beginning of Sec. III, due to the grand-canonical
nature of our formalism, the ion concentrations satisfy automatically the bulk
electroneutrality condition, i.e. $n_{+{\rm b}}-n_{-{\rm b}}+q_{\rm c}n_{\rm
cb}=0$.

Figure 5: (Color online) (a) Polymer mobility $\mu_{\rm p}$ in Eq. (21) versus
the multivalent counterion concentration $n_{\rm cb}$ for various values of
the (a) monovalent salt concentration $n_{+{\rm b}}$, (b) the electrolyte
permittivity $\varepsilon_{\rm w}$, (c) the multivalent charge valency $q_{\rm
c}$, and (d) the polymer permittivity $\varepsilon_{\rm p}$. The model
parameters are $a=1.2$ nm, $\delta=0.1$ nm, and $\sigma_{\rm p}=0.783$ $e/{\rm
nm}^{2}$. The charge valency in (a)-(b) and (d) is $q_{\rm c}=3$, and the salt
concentration in (b)-(d) is $n_{+{\rm b}}=300$ mM.

In Eq. (9), we introduced the partition function of the multivalent ions

$k_{\rm c}(\mathbf{r})=\theta(r-a)e^{-q_{\rm c}\phi_{\rm s}(\mathbf{r})-q_{\rm
c}^{2}\delta G(\mathbf{r})/2}$ (12)

related to the counterion density as $\rho_{\rm c}(r)=n_{\rm cb}k_{\rm c}(r)$,
with the ionic self-energy

$\delta
G(\mathbf{r})=\lim_{\mathbf{r}^{\prime}\to\mathbf{r}}\left[G(\mathbf{r},\mathbf{r}^{\prime})-G_{\rm
b}(\mathbf{r}-\mathbf{r}^{\prime})\right],$ (13)

where $G_{\rm b}(\mathbf{r}-\mathbf{r}^{\prime})=\ell_{\rm
B}/|\mathbf{r}-\mathbf{r}^{\prime}|$ stands for the electrostatic Green’s
function in the bulk region. The self-energy (13) takes into account two
separate electrostatic effects. The first one is the repulsive ionic solvation
force on the multivalent charges originating from the screening deficiency in
the salt-free polymer volume. The second effect is the strongly repulsive
image-charge forces associated with the dielectric contrast between the
polymer and the solvent.

Figure 6: (Color online) (a) Dimensionless multivalent ion density $k_{\rm
c}(r)=n_{\rm c}(r)/n_{\rm cb}$ against the separation distance from the
polymer axis at the model parameters in Figs. 5(a)-(d).

The evaluation of the average potential components (8)-(9) requires the
knowledge of the Green’s function satisfying Eq. (10). The latter can be
solved by exploiting the cylindrical symmetry of the system. The details of
this solution in Fourier space can be found in Refs. 48; 25. The Fourier
expansion of the Green’s function and its bulk limit read

$\displaystyle G(\mathbf{r},\mathbf{r}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{\ell_{\rm
B}}{\pi}\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\;e^{in(\phi-\phi^{\prime})}e^{ik(z-z^{\prime})}$
$\displaystyle\hskip
45.52458pt\times\left[I_{n}(pr_{<})+\Delta_{n}K_{n}(pr_{<})\right]K_{n}(pr_{>});$
$\displaystyle G_{\rm b}(\mathbf{r},\mathbf{r}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{\ell_{\rm
B}}{\pi}\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\;e^{in(\phi-\phi^{\prime})}e^{ik(z-z^{\prime})}$
$\displaystyle\hskip 56.9055pt\times I_{n}(pr_{<})K_{n}(pr_{>}),$

where we used the modified Bessel functions of the first kind $I_{n}(x)$ and
of the second kind $K_{n}(x)$, the auxiliary screening parameter
$p=\sqrt{\kappa^{2}+k^{2}}$, the shortcut notations $r_{<}={\rm
min}(r,r^{\prime})$ and $r_{>}={\rm max}(r,r^{\prime})$ for the radial
coordinates, and the dielectric jump function

$\Delta_{n}=\frac{pI_{n}(ka)I^{\prime}_{n}(pa)-\gamma
kI_{n}(pa)I^{\prime}_{n}(ka)}{-pI_{n}(ka)K^{\prime}_{n}(pa)+\gamma
kK_{n}(pa)I^{\prime}_{n}(ka)}$ (16)

including the dielectric coefficient $\gamma=\varepsilon_{\rm
p}/\varepsilon_{\rm w}$. Injecting Eqs. (III.1)-(III.1) into Eq. (13), the
self-energy in Eq. (12) follows as

$\delta G(r)=\frac{\ell_{\rm
B}}{\pi}\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\;\Delta_{n}K^{2}_{n}(pr).$
(17)

Substituting now the Green’s function (III.1) into Eq. (8), the mean-field
(MF) potential follows as the standard solution of the linear PB equations
around a charged cylinder,

$\phi_{\rm s}(r)=-\frac{2}{q_{\rm c}\kappa\mu_{\rm GC}}\frac{K_{0}(\kappa
r)}{K_{1}(\kappa a)},$ (18)

as expected. Moreover, the correlation component (9) of opposite sign
responsible for CI becomes

$\phi_{\rm c}(r)=4\pi\ell_{\rm B}q_{\rm c}n_{\rm cb}J(r),$ (19)

where we introduced the integral function

$\displaystyle J(r)$ $\displaystyle=$ $\displaystyle K_{0}(\kappa
r)\int_{a}^{r}\mathrm{d}r^{\prime}r^{\prime}k_{\rm c}(r^{\prime})$
$\displaystyle\hskip 45.52458pt\times\left\\{I_{0}(\kappa
r^{\prime})+\frac{I_{1}(\kappa a)}{K_{1}(\kappa a)}K_{0}(\kappa
r^{\prime})\right\\}$ $\displaystyle+\left\\{I_{0}(\kappa
r)+\frac{I_{1}(\kappa a)}{K_{1}(\kappa a)}K_{0}(\kappa r)\right\\}$
$\displaystyle\hskip
8.53581pt\times\int_{r}^{\infty}\mathrm{d}r^{\prime}r^{\prime}k_{\rm
c}(r^{\prime})K_{0}(\kappa r^{\prime}).$

Plugging Eqs. (18)-(19) into Eqs. (7)-(9), the polymer mobility in Eq. (6)
finally takes the form

$\mu_{\rm p}=-\frac{e\sigma_{\rm p}}{\kappa\eta}\frac{K_{0}(\kappa
a^{*})}{K_{1}(\kappa a)}+\frac{eq_{\rm c}n_{\rm
cb}}{\kappa^{2}\eta}\left[\kappa^{2}J(a^{*})-1\right].$ (21)

Equation (21) is the main result of the present work. We emphasize that to our
knowledge this identity is the first analytical beyond-MF mobility formula
accounting for the SC correlations induced by the multivalent charges. The
first term of Eq. (21) corresponds to the MF-level EP mobility, and the second
term of opposite sign is the multivalent counterion contribution responsible
for mobility reversal. From Eq. (21), the critical multivalent cation
concentration at the mobility reversal follows as

$n^{*}_{\rm cb}=\frac{\kappa\sigma_{\rm s}}{q_{\rm c}}\frac{K_{0}(\kappa
a^{*})}{K_{1}(\kappa a)}\frac{1}{\kappa^{2}J(a^{*})-1}.$ (22)

In the simplest case, where the no-slip length vanishes, i.e. $a^{*}=a$, Eqs.
(21)-(22) simplify to

$\displaystyle\mu_{\rm p}$ $\displaystyle=$ $\displaystyle-\frac{e\sigma_{\rm
p}}{\kappa\eta}\frac{K_{0}(\kappa a)}{K_{1}(\kappa a)}$
$\displaystyle+\frac{eq_{\rm c}n_{\rm
cb}}{\kappa^{2}\eta}\left\\{\frac{\kappa}{aK_{1}(\kappa
a)}\int_{a}^{\infty}\mathrm{d}r^{\prime}r^{\prime}k_{\rm
c}(r^{\prime})K_{0}(\kappa r^{\prime})-1\right\\};$ $\displaystyle n_{\rm
cb}^{*}$ $\displaystyle=$ $\displaystyle\frac{\kappa\sigma_{\rm s}}{q_{\rm
c}}\frac{K_{0}(\kappa a)}{K_{1}(\kappa a)}$ (25)
$\displaystyle\times\left\\{\frac{\kappa}{aK_{1}(\kappa
a)}\int_{a}^{\infty}\mathrm{d}r^{\prime}r^{\prime}k_{\rm
c}(r^{\prime})K_{0}(\kappa r^{\prime})-1\right\\}^{-1}.$

Figure 7: (Color online) (a) Polymer mobility (30) and (b)-(c) critical
concentration (32) versus the electrostatic coupling parameter $\Xi$. The
polymer radius is $\bar{a}=10.0$ and the no-slip length is $\bar{\delta}=1.0$.
The remaining parameters are given in the legends.

###  III.2 Characterization of the EP mobility inversion

Figure 5 displays the polymer mobility in Eq. (21) versus the multivalent
counterion concentration, and Fig. 6 illustrates the multivalent ion density
profile in Eq. (12) at various model parameters. These mobility and counterion
density plots should be interpreted together.

The inspection of Figs. 5(a)-(b) indicates that the theoretical mobility
formula (21) agrees qualitatively with the MD simulation results in Figs.
3(b)-(c). Namely, Figs. 6(a) and (b) respectively show that the increment of
the monovalent salt or the solvent permittivity suppressing the interfacial
potential attenuates the counterion adsorption, i.e. $n_{+{\rm b}}\uparrow
k_{\rm c}(r)\downarrow$ and $\varepsilon_{\rm w}\uparrow k_{\rm
c}(r)\downarrow$. In Figs. 5(a)-(b), one sees that due to the resulting
reduction of charge correlations, the larger the monovalent salt concentration
or the liquid permittivity, the larger the multivalent counterion
concentration required for the mobility reversal, i.e. $n_{+{\rm b}}\uparrow
n_{\rm cb}^{*}\uparrow$ and $\varepsilon_{\rm w}\uparrow n_{\rm
cb}^{*}\uparrow$.

Figures 5(c) and 6(c) indicate that charge valency brings an opposite effect
to polymer mobility. Indeed, one sees that multivalent cations of larger
valency exhibit a stronger adsorption and higher interfacial density, i.e.
$q_{\rm c}\uparrow k_{\rm c}(r)\uparrow$. Due to the resulting intensification
of the surface charge correlations, the higher the charge valency, the lower
the counterion concentration at the mobility reversal, i.e. $q_{\rm c}\uparrow
n_{\rm cb}^{*}\downarrow$. We note that this trend is equally in agreement
with the computational result of Fig. 3(a).

Finally, we investigate the surface polarization effects associated with the
low dielectric permittivity of the polymer on its mobility. Figure 6(d) shows
that for $\varepsilon_{\rm p}<\varepsilon_{\rm w}$, the repulsive image-charge
forces embodied by the self-energy term $\delta G(r)$ of Eq. (12) exclude the
multivalent counterions from the surface of the polymer, i.e.
$\varepsilon_{\rm p}\downarrow k_{\rm c}(r)\downarrow$. As a result, Fig. 5(d)
indicates that the lower the polymer permittivity, the higher the required
bulk counterion concentration for the mobility reversal, i.e.
$\varepsilon_{\rm p}\downarrow n^{*}_{\rm cb}\uparrow$. We finally note that
the present theory does not capture the re-entrance regime of the simulation
results in Figs. 3(a)-(b). This limitation of the SC formalism is a
consequence of the virial treatment of the multivalent charges, restricting
the validity of the theory to dilute multivalent counterion concentrations.

###  III.3 Scaling of the mobility and the reversal concentration

In our MD simulation results presented in Sec. II, the effect of the
experimentally controllable physical parameters on the mobility inversion was
explained qualitatively in terms of the electrostatic coupling parameter
$\Xi$. In order to provide a solid mathematical background supporting this
analysis, we investigate here the dependence of the analytically derived EP
mobility formula (21) on this coupling parameter. To this aim, we first scale
all lengths by the GC length as $\bar{r}=r/\mu_{\rm GC}$, and define the
dimensionless salt screening parameter $s=\kappa\mu_{\rm GC}$ and the
electrostatic coupling parameter $\Xi=q_{\rm c}^{2}\ell_{\rm B}/\mu_{\rm GC}$.
The counterion partition function (12) becomes

$k_{\rm c}(\bar{r})=\theta(\bar{r}-\bar{a})e^{-\psi(\bar{r})-\Xi U_{\rm
s}(\bar{r})},$ (26)

with the scaled average potential and self-energy

$\displaystyle\psi(\bar{r})$ $\displaystyle=$
$\displaystyle-\frac{2}{s}\frac{K_{0}(s\bar{r})}{K_{1}(s\bar{a})};$ (27)
$\displaystyle U_{\rm s}(\bar{r})$ $\displaystyle=$
$\displaystyle\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{\mathrm{d}\bar{k}}{2\pi}\;\Delta_{n}K^{2}_{n}(\bar{p}\bar{r}),$
(28)

where we introduced the additional dimensionless parameters
$\bar{a}=a/\mu_{\rm GC}$, $\bar{k}=k\mu_{\rm GC}$,
$\bar{p}=\sqrt{s^{2}+\bar{k}^{2}}$, and

$\Delta_{n}=\frac{\bar{p}I_{n}(\bar{k}\bar{a})I^{\prime}_{n}(\bar{p}\bar{a})-\gamma\bar{k}I_{n}(\bar{p}\bar{a})I^{\prime}_{n}(\bar{k}\bar{a})}{-\bar{p}I_{n}(\bar{k}\bar{a})K^{\prime}_{n}(\bar{p}\bar{a})+\gamma\bar{k}K_{n}(\bar{p}\bar{a})I^{\prime}_{n}(\bar{k}\bar{a})}.$
(29)

In terms of these scaled parameters, the polymer mobility (21) takes the form

$\mu_{\rm p}=\frac{\mu_{\rm ep}}{q_{\rm
c}}\left\\{-\frac{2}{s}\frac{K_{0}(s\bar{a}^{*})}{K_{1}(s\bar{a})}+4\pi\left[I(a^{*})-s^{-2}\right]\Xi\;\bar{n}_{\rm
cb}\right\\},$ (30)

with the dimensionless counterion concentration $\bar{n}_{\rm cb}=\mu_{\rm
GC}^{3}n_{\rm cb}$, and the auxiliary integral

$\displaystyle I(\bar{r})$ $\displaystyle=$ $\displaystyle
K_{0}(s\bar{r})\int_{\bar{a}}^{\bar{r}}\mathrm{d}\bar{r}^{\prime}\bar{r}^{\prime}k_{\rm
c}(\bar{r}^{\prime})$ $\displaystyle\hskip
45.52458pt\times\left\\{I_{0}(s\bar{r}^{\prime})+\frac{I_{1}(s\bar{a})}{K_{1}(s\bar{a})}K_{0}(s\bar{r}^{\prime})\right\\}$
$\displaystyle+\left\\{I_{0}(s\bar{r})+\frac{I_{1}(s\bar{a})}{K_{1}(s\bar{a})}K_{0}(s\bar{r})\right\\}$
$\displaystyle\hskip
8.53581pt\times\int_{\bar{r}}^{\infty}\mathrm{d}\bar{r}^{\prime}\bar{r}^{\prime}k_{\rm
c}(\bar{r}^{\prime})K_{0}(s\bar{r}^{\prime}).$

We first consider the case $\varepsilon_{\rm p}=\varepsilon_{\rm w}$, where
ionic image-charge forces are absent. In Eq. (30), one sees that the
electrostatic coupling parameter $\Xi$ is involved in the positive correlation
component corresponding to the second term in the bracket at two different
levels. First, the coupling parameter $\Xi$ comes into play as a linear
prefactor scaling the multivalent cation concentration $\bar{n}_{\rm cb}$. The
corresponding dependence on the parameter $\Xi$ originates from the
correlations between the polymer charges and the multivalent counterions.
These correlations driving the counterion adsorption onto the polymer surface
give rise to the polymer CI and the mobility reversal. Indeed, one notes that
in the limit of vanishing coupling parameter ($\Xi\to 0$), the correlation
component in Eq. (30) disappears, and the EP mobility tends to its MF
counterpart corresponding to the first term in the bracket.

Secondly, Eq. (III.3) shows that $\Xi$ equally contributes to the mobility
(30) non-linearly as the self-energy magnitude of the ionic partition function
in Eq. (26). This non-linear contribution takes into account the correlations
of the multivalent counterions with their monovalent salt cloud and the
polarization charges. These correlations generate the repulsive solvation and
image-charge forces suppressing the multivalent counterion adsorption.

Figure 7(a) shows that this competition gives rise to the non-uniform trend of
the polymer mobility (30) with respect to the electrostatic coupling
parameter. Namely, one sees that at large enough multivalent cation
concentrations (blue and red), the negative mobility initially governed by MF
electrophoresis rises linearly with $\Xi$, switches from negative to positive
at the CI point, and reaches a peak located at $\Xi\approx 100$. In the regime
of larger coupling strengths, where the repulsive ionic solvation energy in
Eq. (26) significantly excludes multivalent counterions from the polymer
surface, the resulting weakening of the polymer CI suppresses the reversed
mobility and turns the latter from positive back to negative.

From Eq. (30), the critical concentration for the mobility reversal follows as

$\bar{n}_{\rm cb}^{*}=\frac{1}{2\pi
s\Xi}\frac{K_{0}(s\bar{a}^{*})}{K_{1}(s\bar{a})}\frac{1}{I(\bar{a}^{*})-s^{-2}}.$
(32)

Figure 7(b) shows that due to the amplification of the interfacial counterion
adsorption responsible for CI, the rise of $\Xi$ leads to the linear drop of
the characteristic counterion concentration (32) up to $\Xi=\Xi_{\rm c}\approx
100$. At larger $\Xi$ values, where the multivalent counterion adsorption is
significantly weakened by the repulsive ionic solvation forces, the
characteristic concentration $\bar{n}_{\rm cb}^{*}$ reverses its slope and
rises quickly with the coupling parameter. Figure 7(c) indicates that if one
also takes into account the low dielectric permittivity of the polymer, the
addition of the strongly repulsive image-charge interactions to the ionic
solvation forces drops substantially the boundary of the self-energy-dominated
coupling parameter regime, i.e. $\varepsilon_{\rm p}\downarrow\Xi_{\rm
c}\downarrow$.

##  IV Summary and Conclusions

In this work, we have combined correlation-corrected transport theory and
particle-based numerical simulations to characterize the electrostatic
mechanisms behind various experimentally observed features of DNA mobility
inversion. Our MD simulations were run in the $NVT$ ensemble. On the other
hand, our polymer transport theory was derived in the dressed-ion limit Kanduc
_et al._ (2011) of the SCPB formalism Buyukdagli (2020, 2022) incorporating
mutually the weak and SC interactions induced by the mono- and multivalent
ions, respectively. As this theory is based on the ${\rm\mu}VT$ ensemble, the
underlying grand-canonical picture includes the bulk charge reservoir present
in the physical system. The main result of this theoretical approach is the
identity of Eq. (21), providing the first-known-to-us analytical EP mobility
formula taking into account the contribution from SC correlations responsible
for mobility inversion. We summarize below the main predictions of this
analytical formula, equally supported by our numerical simulations, providing
qualitatively equivalent results.

The EP mobility formula and its scaled form in Eq. (30) indicate – as
anticipated – that EP polymer mobility reversal is caused by the excess of
multivalent counterions adsorbed onto the DNA surface. The strength of this
effect originating from the strong coupling of the multivalent counterions and
the DNA charges is set by the magnitude of the coupling parameter $\Xi$ in the
second term of Eq. (30). Our MD simulation results show that as the adsorption
excess is accompanied with the CI of the macromolecule, the resulting like-
charge coion attraction by the anionic polyelectrolyte comes into play as an
additional macroscopic signature of mobility reversal (see Fig. 2).

We have also investigated the effect of the experimentally controllable system
components on the polymer mobility. We found that the increment of the
monovalent salt component weakens charge correlations and the multivalent
charge adsorption onto DNA, resulting in the suppression of the EP mobility
reversal. Consequently, the amount of multivalent counterions $n_{\rm cb}^{*}$
required for the occurrence of mobility inversion increases with the
monovalent salt concentration $n_{+{\rm b}}$, i.e. $n_{+{\rm b}}\uparrow
n_{\rm cb}^{*}\uparrow$. However, as the multivalent counterion-DNA coupling
is enhanced by the increase in the ion valency or the decrease of the
dielectric permittivity, i.e. $q_{\rm c}\uparrow\Xi\uparrow$ and
$\varepsilon_{\rm w}\downarrow\Xi\uparrow$, the increase of the counterion
charge or the decrease of the solvent permittivity lowers the minimum amount
of multivalent charges required for mobility inversion, i.e. $q_{\rm
c}\uparrow n_{\rm cb}^{*}\downarrow$ and $\varepsilon_{\rm w}\downarrow n_{\rm
cb}^{*}\downarrow$. It is noteworthy that the aforementioned predictions are
in qualitative agreement with the observations of the EP polymer transport
experiments Wang  _et al._ (2016, 2018).

In the presence of curved interfaces separating water and solvent-free
dielectric cavities, such as the surface of DNA with permittivity
$\varepsilon_{\rm p}\approx 2-5$ in contact with the high permittivity
electrolyte, the computational cost associated with the presence of an
infinite number of image-charge interactions does not allow the explicit
inclusion of the surface polarization effects into the MD simulation
framework. However, due to the underlying continuum field representation, our
SC-dressed transport theory naturally incorporates the image-charge effects on
the polymer mobility. Based on this generality of the theory, we showed that
the strongly repulsive polarization forces originating from the dielectric
contrast at the polymer-solvent interface repel the multivalent counterions
and reduce their coupling to the DNA charges. Consequently, the dielectric
cavity created by the low permittivity polymer rises the minimum amount of
trivalent counterions required for mobility inversion by almost two orders of
magnitude.

Finally, we have considered two additional peculiarities exclusively
accessible by our simulation approach. In Fig. 3, we showed that the
multivalent counterions triggering mobility inversion at dilute concentrations
suppress the inversion at large concentrations. The corresponding re-entrance
phenomenon previously observed in experiments on DNA condensation Hsiao
(2008); Buyukdagli (2017); Truzzolillo  _et al._ (2018) calls for verification
by EP transport experiments. Additionally, in order to characterize the effect
of the typically extended charge structure of multivalent cations, such as
spermidine, we investigated the impact of the counterion density on the EP
mobility. Our simulations revealed that as the extended structure of the
multivalent cations reduces their surface charge density, also resulting in
their rotational penalty close to the DNA surface, the finite counterion size
weakens the DNA-counterion coupling, therefore suppressing the mobility
reversal. Consequently, the characteristic multivalent counterion
concentration for mobility reversal increases with the ion charge magnitude.
It should be noted that the modeling here simplifies the molecular level
dependencies in nucleic acid – ion interactions, see, e.g., Lipfert  _et al._
(2014); Antila and Sammalkorpi (2014); Antila  _et al._ (2015).

In the present study, the MD framework chosen for our numerical computations
required the simulation of the charged system in the $NVT$ ensemble,
characterized by the constraint of conserved particle number. Future Monte-
Carlo simulations using the particle insertion method would enable the
simulation of the model in the grand-canonical ensemble, thereby allowing
direct quantitative comparison of the theoretical and simulation results.
Moreover, our theoretical and computational models neglect the conformational
polymer fluctuations. It has been experimentally shown that the stiff polymer
approximation holds for polymer lengths extending up to the persistence length
of $l_{\rm p}\approx 30-55$ nm Brunet  _et al._ (2015). Therefore, the
consideration of long polymers necessitates the incorporation of the polymer
conformations into the model. The latter is a highly challenging task, and is
beyond the scope of the present article.

## SUPPLEMENTARY MATERIAL

The Supplementary Material consist of assessment of finite-size effects in the
MD simulations and additional details regarding convergence of the PPPM
electrostatics calculation method in the simulations.

## Acknowledgements

This work was supported by the Academy of Finland through its Centres of
Excellence Programme (2022-2029, LIBER) under project no. 346111 (M.S.). The
work was also supported by Technology Industries of Finland Centennial
Foundation TT2020 grant (T.A-N. and X.Y.). We are grateful for the support by
FinnCERES Materials Bioeconomy Ecosystem. Computational resources by CSC IT
Centre for Finland and RAMI – RawMatters Finland Infrastructure are also
gratefully acknowledged.

## Author Declarations

### Conflict of Interest

The authors have no conflicts of interest to disclose.

## Data Availability

Link to plotted simulations data is provided at https://research.aalto.fi/. If
using the open data, we request acknowledging the authors by a citation to the
original source (this publication).

## Supplementary Material

### Finite-size effects in molecular dynamics simulations

In this section, we explore the influence of the size of the simulation box on
the ion density and the average potential profiles. A major constraint for the
convergence of the system size $V$ is the minimum amount of counterions
required to screen the polymer charges. Therefore, according to the identity
$N_{\rm c}=\rho_{\rm c}V$, low charge concentrations $\rho_{\rm c}$ require a
large system size to keep the counterion number $N_{\rm c}$ above some
treshold value needed to compensate for the DNA charges.

The aforementioned effect can be illustrated by comparing the ion number
density at the boundary of the simulation box with its bulk value. According
to this criterion, finite-size effects can be considered to be negligible if
the normalized multivalent ion density defined as $k_{\rm c}(r)=n_{\rm
c}(r)/n_{\rm cb}$ tends to unity at the system boundary located at $r^{*}$. As
an example, we observed that at high concentrations located in the regime
$n_{\rm cb}\gtrsim 100$ mM, a simulation box of size $L_{x}=24$ nm was able to
provide a sufficient amount of counterions to compensate the PE charges.
Consequently, at this box size, the ion density at the system boundary was
sufficiently close to unity ($k_{\rm c}(r^{*})\gtrsim 0.98$) for finite-size
effects to be negligible.

Figure 8: Normalized trivalent ion density at the bulk concentrations $n_{\rm
cb}=50$ mM (left) and $n_{\rm cb}=5$ mM (right), and for different dimensions
of the simulation box (see the legends) displaying reduction in the finite-
size effects at $L_{x}\gtrsim 60$ nm (left) and $L_{x}\gtrsim 112.8$ nm
(right). The relative permittivity and the longitudinal system size are
$\varepsilon_{\rm w}/\varepsilon_{0}=78$ and $L_{z}=20$ nm.  Figure 9:
Cumulative charge density with trivalent counterions of bulk concentrations
$n_{\rm cb}=1$ mM (left) and $n_{\rm cb}=50$ mM (right), and different
longitudinal system sizes (see the legend) indicating the numerical
convergence at $L_{z}\gtrsim 10$ nm. The relative permittivity is
$\varepsilon_{\rm w}/\varepsilon_{0}=78$, and the transverse system size is
$L_{x}=112.8$ nm (left) and $L_{x}=60$ nm (right).

The left plot of Fig. 8 shows that for the same box size $L_{x}=24$ nm but at
the lower counterion concentration $n_{\rm cb}=50$ mM (orange curve), due to
the decrease of the counterion number below the critical amount required to
screen the DNA charges, the normalized ion density at the system boundary
drops below unity. One sees that by increasing the box size to $L_{x}=60$ nm
(green curve), the density at the boundary increases to unity, indicating the
elimination of the visible finite-size effects from the simulations. In the
right panel of Fig. 8, one sees that decreasing further the counterion
concentration down to $n_{\rm cb}=5$ mM, the box size $L_{x}=60$ nm becomes
insufficient to satisfy the thermodynamic limit in the transverse direction
(green curve). In this case, the increase of the system size up to
$L_{x}=112.8$ nm allowed again to restore the thermodynamic limit by rising
the boundary value of the normalized ion density significantly close to unity
(blue curve). In Table 3, we report the characteristic system sizes enabling
the substantial suppression of the finite-size effects by keeping the boundary
value of the normalized ion density above the value of $0.95$.

An additional finite-size effect is the finite boundary value of the
electrostatic potential. In the physical system, the bulk electroneutrality
condition leads to a vanishing potential in the bulk reservoir. However, in
the simulation box of finite transverse dimensions, the electrostatic
potential always has a finite value at the boundaries of the system. In order
to minimize this finite-size effect, the transverse size of the box should be
chosen significantly larger than the characteristic screening lengths of the
potential. Thus, the elimination of this effect required us to choose the
system size such that the latter satisfies the condition $L_{x}\gg\mu_{\rm
GC}$ for the salt-free liquid and $L_{x}\gg\kappa^{-1}$ in the presence of
salt. The corresponding values are displayed in Table 3.

$n_{\rm cb}$ [mM] |  $L_{x}$ [nm] | max $\kappa^{-1}$ [nm] |  $\mu_{\rm GC}$ [nm] |  $d^{\rm s}$ [nm]
---|---|---|---|---
$0.1-0.5$ | 240 | 25 | 0.1 | 3.0
$1-5$ | 112.8 | 7.9 | 0.1 | 3.0
$10-50$ | 60 | 2.5 | 0.1 | 1.5
$>100$ | 24 | 0.8 | 0.1 | 1.5
Table 3: Characteristic lengths and parameters: the horizontal box size
$L_{x}$, the screening length $\kappa^{-1}$, the GC length $\mu_{\rm GC}$, and
the PPPM distance $d^{\rm s}$ for various trivalent salt concentrations and
the relative permittivity $\varepsilon_{\rm w}/\varepsilon_{0}=78$.

Finally, in order to eliminate finite-size effects associated with the length
of the polymer, we had to identify the characteristic polymer length above
which the results remain unaffected by an increase in the molecule length. In
Fig. 9, we show that the increase of the polymer length above $L_{z}\approx
10$ nm does not lead to a significant change of the cumulative charge density.
Thus, in order to ensure that our ion condensation results are free of finite
length effects, we always set the polymer length to $L_{z}=20$ nm.

### Additional details regarding convergence of the PPPM electrostatics
calculation method

We provide here details regarding convergence considerations of the
particle–particle particle-mesh (PPPM) method used to evaluate the pairwise
electrostatic interactions between the charges in the polymer-ion complex. The
PPPM approach splits the electrostatic Coulomb interactions into a short-range
and a long-range component. The short-range component defined as the
interactions between each charge and its neighbors located within a distance
$d^{\rm s}$ is evaluated directly in the method. The long-range interactions
of the central charge with the other charges located at a larger distance
$r>d^{\rm s}$ are evaluated by mapping the interacting particles into a grid
and carrying-out a fast Fourier transform Brown  _et al._ (2012). The Key
parameter for the convergence of the calculation is $d^{\rm s}$. Grid spacing
is generated automatically from $d^{\rm s}$ by LAMMPS. For computational
efficiency, the parameter $d^{\rm s}$ is often in coarse-grained modelling set
to a value slightly larger than the particle size. However, for elevated local
charge density such as in this study, the choice can lead to loss of
convergence when $d^{\rm s}>4$ nm. We carefully checked the ion distribution
convergence in our work. Based on the outcome of the convergence checking, in
our simulations, we use $d^{\rm s}=1.5$ nm for system sizes in the range
$L_{x}\leq 60$ nm. In larger systems where the simulations are significantly
more time consuming if extensive weight is on the mesh part, we use a larger
value of $d^{\rm s}=3$ nm. This allows fast numerical computation by keeping
practically the same accuracy. The ion distribution remains the same as
$d^{\rm s}=1.5$ nm. The corresponding values of the parameter $d^{\rm s}$ are
reported in Table 3. The parameter of relative error was set to $10^{-5}$ and
the stencil size was set to 5. No other PPPM calculation parameters were
modified from LAMMPS algorithm implementation defaults.

## References

  * Grosberg  _et al._ (2002) A. Y. Grosberg, T. Nguyen, and B. Shklovskii, Reviews of modern physics 74, 329 (2002).
  * Levin (2002) Y. Levin, Reports on progress in physics 65, 1577 (2002).
  * Angelini  _et al._ (2003) T. E. Angelini, H. Liang, W. Wriggers, and G. C. Wong, Proceedings of the National Academy of Sciences 100, 8634 (2003).
  * Besteman  _et al._ (2007) K. Besteman, K. Van Eijk, and S. Lemay, Nature Physics 3, 641 (2007).
  * Besteman  _et al._ (2004) K. Besteman, M. A. Zevenbergen, H. A. Heering, and S. G. Lemay, Physical review letters 93, 170802 (2004).
  * Butler  _et al._ (2003) J. C. Butler, T. Angelini, J. X. Tang, and G. C. Wong, Physical review letters 91, 028301 (2003).
  * Ray and Manning (1997) J. Ray and G. S. Manning, Macromolecules 30, 5739 (1997).
  * Podgornik and Parsegian (1998) R. Podgornik and V. Parsegian, Physical review letters 80, 1560 (1998).
  * Kornyshev and Leikin (1997) A. Kornyshev and S. Leikin, The Journal of chemical physics 107, 3656 (1997).
  * Kornyshev and Leikin (1999) A. Kornyshev and S. Leikin, Physical review letters 82, 4138 (1999).
  * Tang  _et al._ (1996) J. X. Tang, S. Wong, P. T. Tran, and P. A. Janmey, Berichte der Bunsengesellschaft für physikalische Chemie 100, 796 (1996).
  * Bloomfield (1997) V. A. Bloomfield, Biopolymers 44, 269 (1997).
  * Teif and Bohinc (2011) V. B. Teif and K. Bohinc, Progress in biophysics and molecular biology 105, 208 (2011).
  * Van der Heyden  _et al._ (2006) F. H. Van der Heyden, D. Stein, K. Besteman, S. G. Lemay, and C. Dekker, Physical review letters 96, 224502 (2006).
  * Rippe  _et al._ (2008) K. Rippe, J. Mazurkiewicz, and N. Kepper, DNA Interactions with Polymers and Surfactants. John Wiley & Sons, Isnc , 135 (2008).
  * Kaper  _et al._ (2003) H. J. Kaper, H. J. Busscher, and W. Norde, Journal of Biomaterials Science, Polymer Edition 14, 313 (2003).
  * Danger  _et al._ (2007) G. Danger, M. Ramonda, and H. Cottet, Electrophoresis 28, 925 (2007).
  * Kabanov and Kabanov (1995) A. Kabanov and V. Kabanov, Bioconjugate chemistry 6, 7 (1995).
  * Kabanov  _et al._ (1996) V. Kabanov, A. Yaroslavov, and S. Sukhishvili, Journal of controlled release 39, 173 (1996).
  * Holm  _et al._ (2001) C. Holm, P. Kékicheff, and R. Podgornik,  _Electrostatic effects in soft matter and biophysics_ , Vol. 46 (Springer Science & Business Media, 2001).
  * Schoch  _et al._ (2008) R. B. Schoch, J. Han, and P. Renaud, Reviews of modern physics 80, 839 (2008).
  * Jain  _et al._ (2016) M. Jain, H. E. Olsen, B. Paten, and M. Akeson, Genome biology 17, 1 (2016).
  * Thomas  _et al._ (2016) T. Thomas, H. Tajmir-Riahi, and T. Thomas, Amino Acids 48, 2423 (2016).
  * Buyukdagli (2020) S. Buyukdagli, The Journal of Physical Chemistry B 124, 11299 (2020).
  * Buyukdagli (2022) S. Buyukdagli, Langmuir 38, 022302 (2022).
  * Wang  _et al._ (2016) Y. Wang, R. Wang, B. Cao, Z. Guo, and G. Yang, Sci. Rep. 6, 1 (2016).
  * Wang  _et al._ (2018) Y. Wang, R. Wang, T. Gao, and G. Yang, Polymers 10, 244 (2018).
  * Kanduč  _et al._ (2010) M. Kanduč, A. Naji, J. Forsman, and R. Podgornik, The Journal of chemical physics 132, 124701 (2010).
  * Kanduc  _et al._ (2011) M. Kanduc, A. Naji, J. Forsman, and R. Podgornik, Physical Review E 84, 011502 (2011).
  * Forrey and Muthukumar (2006) C. Forrey and M. Muthukumar, Biophy. J. 91, 25 (2006).
  * Mandelkern  _et al._ (1981) M. Mandelkern, J. G. Elias, D. Eden, and D. M. Crothers, J. Mol. Biol. 152, 153 (1981).
  * Weeks  _et al._ (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
  * Vahid  _et al._ (2022) H. Vahid, A. Scacchi, X. Yang, T. Ala-Nissila, and M. Sammalkorpi, The Journal of Chemical Physics 156, 214906 (2022).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, _et al._ , Comput. Phys. Commun. 271, 108171 (2022).
  * Hockney and Eastwood (2021) R. W. Hockney and J. W. Eastwood,  _Computer simulation using particles_ (crc Press, 2021).
  * Nosé (1984) S. Nosé, Mol. Phys. 52, 255 (1984).
  * Hoover (1985) W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
  * Netz (2001) R. R. Netz, Eur. Phys. J. E 5, 557 (2001).
  * Hsiao (2008) P.-Y. Hsiao, J. oPhys. Chem. B 112, 7347 (2008).
  * Aksimentiev (2010) A. Aksimentiev, Nanoscale 2, 468 (2010).
  * Bhattacharjee (2016) S. Bhattacharjee, J. Control Release 235, 337 (2016).
  * Yamaguchi and Kobayashi (2016) A. Yamaguchi and M. Kobayashi, Colloid Polym. Sci. 294, 1019 (2016).
  * Galla  _et al._ (2014) L. Galla, A. J. Meyer, A. Spiering, A. Sischka, M. Mayer, A. R. Hall, P. Reimann, and D. Anselmetti, Nano Lett. 14, 4176 (2014).
  * Antila  _et al._ (2017) H. S. Antila, P. R. Van Tassel, and M. Sammalkorpi, J. Chem. Phys. 147, 124901 (2017).
  * Buyukdagli (2017) S. Buyukdagli, Physical Review E 95, 022502 (2017).
  * Truzzolillo  _et al._ (2018) D. Truzzolillo, S. Sennato, S. Sarti, S. Casciardi, C. Bazzoni, and F. Bordi, Soft Matter 14, 4110 (2018).
  * Buyukdagli  _et al._ (2011) S. Buyukdagli, M. Manghi, and J. Palmeri, The Journal of chemical physics 134, 074706 (2011).
  * Lipfert  _et al._ (2014) J. Lipfert, S. Doniach, R. Das, and D. Herschlag, Annual Review of Biochemistry 83, 813 (2014).
  * Antila and Sammalkorpi (2014) H. S. Antila and M. Sammalkorpi, J. Phys. Chem. B 118, 3226 (2014).
  * Antila  _et al._ (2015) H. S. Antila, M. Härkönen, and M. Sammalkorpi, Phys. Chem. Chem. Phys. 17, 5279 (2015).
  * Brunet  _et al._ (2015) A. Brunet, C. Tardin, L. Salomé, P. Rousseau, N. Destainville, and M. Manghi, Macromolecules 48, 3641 (2015).
  * Brown  _et al._ (2012) W. M. Brown, A. Kohlmeyer, S. J. Plimpton, and A. N. Tharrington, Computer Physics Communications 183, 449 (2012).


Eliminating black holes singularity in Brane world Cosmology.
Poula Tadros1 and Mohamed Assaad Abdel-Raouf2

1 Department of Physics and Astronomy, University of Turku, 20500,
Turku,Finland.

2 Physics Department, Faculty of Science, Ain Shame University, Cairo, Egypt.

Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In recent works, a construction was proposed resulting emergent universes
inside black holes. This result can be obtained from a 4D black hole embedded
in a 5D spacetime with the fifth dimension compactified on a circle $[0,2]$ (0
and 2 are identified) such that the two branes are at 0 and 1. In the present
work, we study this setup by deriving particles’ equations of motion in the
new universes, based on redefining energy and angular momentum. This leads to
disappearance of the singularity in centers of black holes in classical
General Relativity.

Keywords: brane world cosmology, string cosmology, Black holes.

##  1 Introduction

Black holes were introduced as solutions for Einstien’s field equations by K.
Schwarzchild (For English translation see Ref. $[1]$). Schwarzchild solution
represents a non rotating, non charged spherically symmetric black hole.
However, this solution suffers from a singularity at the center of the black
hole. In the context of string theory it is argued that black holes emerge
naturally (Reviewed in Ref. $[2,3]$).
Recently, It was suggested that inside the horizon of a black hole there is an
emergent cosmology i.e. a new universe $[4]$ using a topological soliton
called S-branes $[5]$ to violate the null energy condition so eliminating the
singularity.

Cosmological perturbations were studied in the emergent universe $[6]$, and a
mechanism of reheating $[7]$ was proposed using the decay of the S-brane as an
initiator for the big bang and the radiation dominated era. This construction
solves the black hole information paradox as the information falling into the
black hole emerges in the new universe without any loss.

In the present work, we study a particular metric of a 4D black hole embedded
in a 5D spacetime within the context of Brane world cosmology. Introducing an
additional compact dimension as the fifth dimension we prove that particles
falling into the black hole from the visible or the hidden brane do not
experience a singularity with the assumption that the cosmological constant is
negative or zero inside the horizon. Thus, an emergent cosmology can be
defined inside the horizon between the two branes.

##  2 Mathematical Analysis

In this setup a 4D black hole in our universe is viewed as embedded in a 5D
world. The metric for this black hole is defined as in Ref. [8] to be

$ds^{2}=\frac{\Lambda
y^{2}}{3}(f(r)dt^{2}-\frac{dr^{2}}{f(r)}-r^{2}(d\theta^{2}+\sin^{2}(\theta)d\phi^{2}))-dy^{2}$
(1)

where $(t,r,\theta,\phi)$ are the usual spherical coordinates, y is the
coordinate through the fifth dimension(the compact dimension), and
$f(r)=1-\frac{2m}{r}-\frac{\Lambda r^{2}}{3}$ where $\Lambda$ is the
cosmological constant.

By spherical symmetry of the system on each y=constant plane, we can set
$\theta=\frac{\pi}{2}$.

The first geodesic equation (the equation for t) is given by

$\frac{d^{2}t}{d\tau^{2}}+\frac{f^{\prime}(r)}{2f(r)}\frac{dt}{d\tau}\frac{dr}{d\tau}+\frac{1}{y}\frac{dt}{d\tau}\frac{dy}{d\tau}=0,$
(2)

where f’(r) is the derivative of f with respect to r.
This equation can be integrated to give

$\frac{dt}{d\tau}=\frac{E}{y\sqrt{|f|}}$ (3)

where E is a constant which represents the total energy of the system in our
model.

The second geodesic equation(the equation for $\phi$)is given by

$\frac{d^{2}\phi}{d\tau^{2}}+\frac{1}{r}\frac{d\phi}{d\tau}\frac{dr}{d\tau}+\frac{1}{y}\frac{d\phi}{d\tau}\frac{dy}{d\tau}=0,$
(4)

Which can be integrated to

$\frac{d\phi}{d\tau}=\frac{L}{ry},$ (5)

where L is a constant which represents the total angular momentum of the
system in our model.
The equation for r is given by

$\frac{d^{2}r}{d\tau^{2}}+\frac{f(r)f^{\prime}(r)}{2}(\frac{dt}{d\tau})^{2}-\frac{f^{\prime}(r)}{2f(r)}(\frac{dr}{d\tau})^{2}+rf(r)(\frac{d\phi}{d\tau})^{2}+\frac{1}{y}\frac{dr}{d\tau}\frac{dy}{d\tau}=0.$
(6)

After using equations (3),(5) and dividing by $E^{2}$ (recalling that the
impact parameter b is given by $b=\frac{L}{E}$), we get

$\frac{1}{E^{2}}\frac{d^{2}r}{d\tau^{2}}+\frac{f^{\prime}(r)}{2y^{2}}-\frac{1}{2E^{2}}\frac{dr}{d\tau}\frac{d}{d\tau}(\ln|f(r)|)+\frac{b^{2}f(r)}{ry^{2}}+\frac{1}{E^{2}}\frac{dr}{d\tau}\frac{d}{d\tau}(\ln|y|)=0.$
(7)

The equation for y is given by

$\frac{d^{2}y}{d\tau^{2}}+\frac{\Lambda
yf(r)}{3}(\frac{dt}{d\tau})^{2}-\frac{\Lambda
y}{3f(r)}(\frac{dr}{d\tau})^{2}-\frac{\Lambda
r^{2}y}{3}(\frac{d\phi}{d\tau})^{2}=0.$ (8)

With similar calculations as the previous equation we get

$\frac{d^{2}y}{d\tau^{2}}-\frac{\Lambda
y}{3f(r)}(\frac{dr}{d\tau})^{2}=\frac{\Lambda E^{2}}{3y}(b^{2}-1).$ (9)

In the Ekpyrotic scenario where there exist two nearly parallel D-branes one
containing our universe(the visible brane), and another brane(the hidden
brane). We can redefine y such that the distance between the two branes is 1.

Since particles (except gravitons) are restricted to move on one of the
branes, we have two possible values for the impact parameter(b) assuming that
the black hole is in the visible brane. If the particle comes from the Hidden
brane to the position ”above” the black hole, then $b=1$. if the light comes
from the visible brane, then $b=0$. Let us examine both cases.

Case 1: $b=1$.
Equations (7) and (9) will be

$\frac{1}{E^{2}}\frac{d}{d\tau}(\ln(|\frac{dr}{d\tau})|)-\frac{1}{2E^{2}}\frac{d}{d\tau}(\ln(|f(r)|))+\frac{1}{E^{2}}\frac{d}{d\tau}(\ln(|y|))+\frac{f}{ry^{2}}+\frac{f^{\prime}(r)}{2y^{2}}=0.$
(10)

and

$\frac{d^{2}y}{d\tau^{2}}=\frac{\Lambda y}{3f(r)}(\frac{dr}{d\tau})^{2}.$ (11)

Using coordinate time instead of the parameter $\tau$ we get

$\frac{d}{dr}(\frac{1}{\sqrt{|f(r)|}y}\frac{dy}{dt})=\frac{\Lambda}{3f(r)\sqrt{|f(r)|}}(\frac{dr}{dt})^{2}.$
(12)

$\frac{d}{dt}(\ln(|\frac{E}{f(r)}\frac{dr}{dt}|)=\frac{\sqrt{|f(r)|}}{Ey^{2}}(\frac{f(r)}{r}+\frac{f^{\prime}(r)}{2}).$
(13)

These are integrated to

$\frac{dr}{dt}=\frac{f(r)}{E}\exp(\frac{1}{E}\int\frac{\sqrt{|f(r)|}}{y^{2}}(\frac{f(r)}{r}+\frac{f^{\prime}(r)}{2})dr)$
(14)

$y=\exp(\int dt(\sqrt{|f(r)|}\int
dt^{\prime}(\frac{\Lambda}{3f\sqrt{|f(r)|}}(\frac{dr}{dt})^{2}))$ (15)

We examine the limits of y as r tends to zero and infinity.
Far before the interaction the velocity of the particle is constant i.e.
$\frac{dr}{dt}\rightarrow 0$ and $f\rightarrow-\infty$ so

$\lim_{r\to\infty}(y)=1$ (16)

which is exactly what is expected as the particle is restricted to be on the
brane with y=1 before the interaction.

The limit as r tends to zero can be deduced from the fact that the velocity of
particles tend to zero on approaching the singularity ${[9]}$ i.e.
$\frac{dr}{dt}\rightarrow 0$ and $f\rightarrow-\infty$.

$\lim_{r\to 0}(y)=1,$ (17)

i.e. a particle in the hidden brane will remain in the hidden brane in the
discussed limits.

Case 2: $b=0$.
Equations (7) and (9) give

$\frac{1}{E^{2}}\frac{d^{2}r}{d\tau^{2}}+\frac{f^{\prime}(r)}{2y^{2}}-\frac{1}{2E^{2}}\frac{dr}{d\tau}\frac{d}{d\tau}(\ln|f(r)|)+\frac{1}{E^{2}}\frac{dr}{d\tau}\frac{d}{d\tau}(\ln(y))=0$
(18)

and

$\frac{d^{2}y}{d\tau^{2}}-\frac{\Lambda
y}{3f(r)}(\frac{dr}{d\tau})^{2}=-\frac{\Lambda E^{2}}{3y}$ (19)

Using coordinate time and with similar calculations as the first case we get

$y=\exp(\int dt\sqrt{|f(r)|}[\int
dt\frac{\Lambda\sqrt{|f(r)|}}{3}(\frac{1}{E^{2}}\exp(-\int
dr\frac{f(r)f^{\prime}(r)}{(\frac{dr}{dt})^{2}})-1)]).$ (20)

Examining the same limits we get

$\lim_{r\to\infty}(y)=0,$ (21)

because as r approaches infinity $f\rightarrow-\infty$, $f`\rightarrow-\infty$
and $\frac{dr}{dt}\rightarrow 0$, as expected as particles are restricted to
be on the visible brane.

$\lim_{r\to 0}(y)=\left\\{\begin{array}[]{cc}1&\quad\Lambda=0\\\
0&\quad\Lambda<0\end{array}\right.$ (22)

This is because $r\rightarrow 0$ implies $f\rightarrow-\infty$,
$f`\rightarrow\infty$. This converges only if $\Lambda\leq 0$ inside the
horizon.

##  3 Discussion

We have shown that in the case of a non rotating black hole in 4D in 5D
spacetime with one compact dimension, the setup of a model where there are two
branes on which all particles (except gravitons) are restricted to move, a
particle falling radially into the black hole ends up in the bulk between the
two branes.

This suggests that a new emergent cosmology can be defined in the bulk inside
the horizon. The shape of the emergent universe is determined by the solution
of eq.(20). The most important feature is that this universe is free of
physical singularities. Moreover, particles in the hidden brane never reach
the visible brane, that is why we do not observe white holes in our universe
in a finite time.

In the new universe the particles velocity decreases as they approach a
certain point ($r=0$), this can be interpreted as an expanding universe with
the center of expansion at $r=0$. From the viewpoint of the hidden brane the
black hole looks like a potential well and the emergent universe is trapped
inside it. while the visible brane sees the black hole as a bridge to the
hidden brane through the bulk. Thus the two viewpoints see a universe in the
bulk. Adopting the setup in which the two branes have universe-antiuniverse
pair (universe on the visible brane and antiuniverse on the hidden brane),
then on any of the two branes one sees a universe, or an antiuniverse. This
setup has useful features such as that the matter-antimatter asymmetry is
avoided, and the system is CPT symmetric [10].

These results agree with Ref. [4] that an emergent universe can be defined
inside the horizon from the point of view of a three dimensional observer.
However, in this paper we argue that it is in fact in an additional dimension.
Thus, we can identify the emergent universe discussed here and the emergent
universe in Ref. [4], this can be done by adopting the following timeline: At
the beginning of the formation of the black hole the high density deforms the
two branes in a certain region until they collide. An S-brane was formed in
the interface of the collision, then the S-brane decays inducing reheating as
in Ref. [5].

##  4 References.

$[1]$ K. Schwarzschild, Sitzungsberichte der Königlich Preuss. Akad. der
Wissenschaften $(1999)$ 189\.
$[2]$ S. De Haro, J. van Dongen, Visser, M. and J. Butterfield, Stud. Hist.
Philos. Sci. Part B - Stud. Hist. Philos. Mod. Phys. $69(2020)$ 82\.
$[3]$ J. van Dongen, S. De Haro, M. Visser, and J. Butterfield, Stud. Hist.
Philos. Sci. Part B - Stud. Hist. Philos. Mod. Phys. $69(2020)112$.
$[4]$ R. Brandenberger, L. Heisenberg, and J. Robnik, J. High Energy Phys.
2021 (2021) 90.
$[5]$ M. Gutperle, and A. Strominger, J. High Energy Phys. 6 (2002) 395.
$[6]$ R. Brandenberger, and Z. Wang, Phys. Rev. D 102 (2020) 23516.
$[7]$ R. Brandenberger, K. Dasgupta, and Z. Wang, Phys. Rev. D 102 (2020)
063514.
$[8]$ P. Wesson, H. Mashhoon, and H. Liu, Modern Physics Letters A, 12 (1997)
2309.
$[9]$ A. T. Augousti, P. Gusin, B. Kuśmierz, J. Masajada and A. Radosz, Gen.
Relativ. Gravit. (2018) 50.
$[10]$ L. Boyle,K. Finn, and N. Turok, Phys. Rev. Lett. 121 (2018) 1.