# Impact of magnetic field on giant dipole resonance of 40Ca using the EQMD
model

Y. T. Cao(曹雅婷)  Key Laboratory of Nuclear Physics and Ion-beam Application
(MOE), Institute of Modern Physics, Fudan University, Shanghai 200433, China
X. G. Deng(邓先概) ID 111<EMAIL_ADDRESS>Key Laboratory of Nuclear
Physics and Ion-beam Application (MOE), Institute of Modern Physics, Fudan
University, Shanghai 200433, China  Shanghai Research Center for Theoretical
Nuclear Physics， NSFC and Fudan University, Shanghai 200438, China  Y. G.
Ma(马余刚) ID<EMAIL_ADDRESS>Key Laboratory of Nuclear Physics and Ion-
beam Application (MOE), Institute of Modern Physics, Fudan University,
Shanghai 200433, China  Shanghai Research Center for Theoretical Nuclear
Physics， NSFC and Fudan University, Shanghai 200438, China

###### Abstract

By taking into account the magnetic field in the extended quantum molecular
dynamics model (EQMD), we analyzed its effects on giant dipole resonance (GDR)
by studying the responses and strengths of the dipole moments. The selected
system is the 40Ca nucleus which is excited through the Coulomb interaction by
16O. Particle acceleration term in Liénard-Wiechert potential is discussed
which, however, has small impact on magnetic field. The peak energy, strength
and width of GDR, temperature, and angular momentum of 40Ca as a function of
beam energy are investigated. It is found that the magnetic field enhances the
peak energy, strength and width of GDR which is not only due to the
temperature effects but also due to the enhancement of the angular momentum of
nucleus. At beam energy E $>$ 200 MeV/nucleon, magnetic field maintains a
constant value for the strength of GDR. The work sheds light on examining
important roles of the magnetic field on nuclear structure in low-intermediate
energy heavy-ion collisions.

###### pacs:

25.70.-z, 24.10.Lx, 21.30.Fe

##  I Introduction

The electric and magnetic fields generated in heavy-ion collisions as well as
their effects on observable are of research interests for some times [1]. The
strongest magnetic field in the universe could be produced in ultra-
relativistic heavy-ion collisions by theoretical prediction [2, 3]. Some
observational effects and its properties for this strong field have been
discussed in relativistic heavy-ion collisions [4, 5, 6, 7, 8], but less
discussed in low-intermediate energy domain [9, 10, 11, 12, 13]. Even though
the magnetic field is weaker in this energy region in contrast to ultra-
relativistic collisions, but its lifetime is longer, its roles on nuclear
structure or reaction should be further surveyed.

Giant resonances (GR) are resonance states built in the nucleus [14], which
are of particular importance because they provide reliable information about
the bulk behavior of nuclear many-body systems [15, 16]. One of the best-known
modes of giant resonances is isovector giant dipole resonance (IVGDR) in which
neutrons and protons move collectively relative to each other in a nucleus.
IVGDR is an excellent tool for systematical investigation on collective
properties of nuclei because it can give crucial clues to understand nuclear
structure and collective dynamics. Furthermore, IVGDR can also give rise to a
dynamic electric-dipole (E1) moment [17], which reflects asymmetry information
in nuclear equation of state (EoS) [18, 19, 20, 21, 22]. Experimentally, the
IVGDR can be induced by various ways such as photoabsorption [23, 24, 25, 26],
inelastic scattering [27, 28] as well as ultra-peripheral Coulomb excitation
[29] etc., and the data can be fitted using the Lorentzian curve from which
the key parameters consisting of the peak energy, the width, and the peak
strength can be extracted. The peak energy is a good probe of the symmetry
energy $E_{sym}$ around and below the saturation density [19], and it can also
provide direct information about nuclear size and the EoS [14, 20] as well as
the shear viscosity [30, 31, 32, 33]. Meanwhile, it is well established that
the width and peak splitting can be taken as a direct experimental tool to
measure the nuclear deformation [34, 35, 36, 37].

Most studies of GRs are built on the ground state of nuclei. In Ref. [38],
however, it is conjectured that GRs could be built on other nuclear states and
their characteristics would not depend on the detailed microscopic structure
which gives rise to the enthusiasm to study the properties of GRs built on
excited states. In addition, the IVGDR is considered to be one of the most
useful tools to study properties of rotating nuclei at finite temperature
[39]. The resonance properties of IVGDR, such as peak energy, width and peak
strength, are investigated as a function of temperature, contributing on
understanding of the symmetry energy term in the nuclear equation of state
[40, 41].

Several microscopic approaches have been employed to study giant resonances
such as separable random-phase-approximation (SRPA) [42, 43], time-dependent
Skyrme-Hartree-Fock (TSHF) method [44, 45], relativistic quasi-particle random
phase approximation (RQRPA) [48, 46, 47] and some transport models, for
instance, Boltzmann-Uhlenbeck-Uehling (BUU) model [49, 50, 51] and quantum
molecular dynamics (QMD) model [36]. Among them, the extended quantum
molecular dynamics (EQMD) model [52] which was improved based upon the spirit
of QMD has been successfully applied to investigate IVGDR [53, 54].

A variety of experimental and theoretical studies have been carried out by
far, but there are a few discussions about the Coulomb excitation, especially
in the excited states. Besides, the impact of magnetic field is seldom
investigated at low and intermediate energy heavy-ion reaction which is far
from pion production energies [10]. So the aim of this paper is to study the
effect of electromagnetic field on the IVGDR of the excited 40Ca through
ultra-peripheral 40Ca + 16O collision at energy ranging from 50 to 500
MeV/nucleon.

The organization of the paper is as follows: In Sect. II, we give a brief
introduction of our simulation model and method including electromagnetic
field and the algorithm of IVGDR. Results of magnetic field effects are
discussed in Sect. III. Finally, conclusion is given in Sect. IV.

##  II Model and methodology

###  II.1 EQMD model

Quantum molecular dynamics (QMD) model has been widely used to study
intermediate energy nuclear reactions and nuclear fragmentation [55]. It has
the microscopic basis and high flexibility of this type of models [56] that
make them successful application in studies of giant resonances such as giant
dipole resonance (GDR), pygmy dipole resonance (PDR) as well as giant monopole
resonance (GMR) [57, 58, 59]. Nevertheless, standard QMD model shows
insufficient stability, for which the phase space obtained from the samples of
Monte Carlo is usually not in the lowest point of energy [52]. To compensate
for this shortcoming, two features are of great importance, one is the
capability to describe nuclear ground state, another is the stability of
nucleus in the model description [52]. The above both features have been well
solved in the EQMD model. As we know, in order to cancel the zero-point energy
caused by the wave packet broadening in the standard QMD, the friction cooling
process can be used to keep the mathematical ground state, but the Pauli
principle is broken [52]. However, the EQMD model embodies the Pauli potential
into effective interactions, making saturation property and $\alpha$-cluster
structures can be obtained after energy cooling [53]. Besides, the EQMD model
takes the kinetic-energy term of the momentum variance of wave packets in the
Hamiltonian into account, which is ignored as the spurious constant term in
the standard QMD [60, 61]. Furthermore, in the EQMD model, the Hamiltonian
introduces the wave packet width as a complex variable, and treats it as an
independent dynamic variable, and in contrast, a uniform and static width is
applied for all nucleons in the standard QMD model.

In the EQMD model, the wave packet of nucleon is taken as the form of
Gaussian-like, of which the total wave function of the system is treated as a
direct product of all nucleons [52]

$\Psi=\prod_{i}\varphi(\bm{r}_{i})\,,$ (1)

$\varphi(\bm{r}_{i})=(\frac{\nu_{i}+\nu^{\ast}_{i}}{2\pi})^{3/4}{\rm
exp}[-\frac{\nu_{i}}{2}(\bm{r}_{i}-\bm{R}_{i})^{2}+\frac{i}{\hbar}\bm{P}_{i}\cdot\bm{r}_{i}]\,,$
(2)

where $\bm{R}_{i}$ and $\bm{P}_{i}$ are the centers of position and momentum
of the $i$-th wave packet, respectively. The Gaussian width $\nu_{i}$ is
introduced as

$\nu_{i}\equiv\frac{1}{\lambda_{i}}+i\delta_{i}\,,$ (3)

where $\lambda_{i}$ and $\delta_{i}$ are dynamical variables in the process of
initialization.

The energy-minimum state is considered as the ground state of initial nucleus.
At the first step, a random configuration is given to each nucleus.
Furthermore, under the time-dependent variation principle (TDVP) [62],
propagation of each nucleon can be described as [52]

$\begin{split}\bm{\dot{R}}_{i}=\frac{\partial{H}}{\partial{\bm{P}_{i}}}+\mu_{\bm{R}}\frac{\partial{H}}{\partial{\bm{R}_{i}}},\bm{\dot{P}}_{i}=-\frac{\partial{H}}{\partial{\bm{R}_{i}}}+\mu_{\bm{P}}\frac{\partial{H}}{\partial{\bm{P}_{i}}},\\\
\frac{3\hbar}{4}\dot{\lambda}_{i}=-\frac{\partial{H}}{\partial{\delta_{i}}}+\mu_{\lambda}\frac{\partial{H}}{\partial{\lambda_{i}}},\frac{3\hbar}{4}\dot{\delta}_{i}=\frac{\partial{H}}{\partial{\lambda_{i}}}+\mu_{\delta}\frac{\partial{H}}{\partial{\delta_{i}}},\end{split}$
(4)

where $H$ is the expected value of the Hamiltonian, and $\mu_{\bm{R}}$,
$\mu_{\bm{P}}$, $\mu_{\lambda}$ and $\mu_{\delta}$ are various friction
coefficients. During the friction cooling process, these coefficients are
negative, making the system goes to its local energy minimum point [36]. In
contrast, in the subsequent nuclear reaction simulation stage, these
coefficients are zero to maintain the energy conservation of the system [63].

As a consequence, the expected value of Hamiltonian can be expressed as

$\begin{split}H&=\left\langle\Psi\left|\sum_{i}-\frac{\hbar^{2}}{2m}\bigtriangledown^{2}_{i}-\hat{T}_{zero}+\hat{H}_{int}\right|\Psi\right\rangle\\\
&=\sum_{i}\frac{\bm{P}^{2}_{i}}{2m}+\frac{3\hbar^{2}(1+\lambda^{2}_{i}\delta^{2}_{i})}{4m\lambda_{i}}-T_{zero}+H_{int},\end{split}$
(5)

where the first, second and third term are the center momentum of the wave
packet, the contribution of the dynamic wave packet, and the zero point
center-of-mass kinetic energy $-T_{zero}$, respectively. The first term can be
expressed as $\left\langle\hat{\bm{p}}_{i}\right\rangle^{2}/2m$, the second
term can be treated as
$\left(\left\langle\hat{\bm{p}}_{i}^{2}\right\rangle-\left\langle\hat{\bm{p}}_{i}\right\rangle^{2}\right)/2m$,
and the form of the third term can be found in details in Ref. [52].

For the effective interaction $H_{int}$, it consists of the Skyrme potential,
the Coulomb potential, the symmetry energy, and the Pauli potential as follows

$H_{int}=H_{Skyrme}+H_{Coulomb}+H_{Symmetry}+H_{Pauli}.$ (6)

The form of Skyrme interaction is the simplest, written as

$H_{Skyrme}=\frac{\alpha}{2\rho_{0}}\int\rho^{2}(\bm{r})d^{3}r+\frac{\beta}{(\gamma+1)\rho_{0}^{\gamma}}\int\rho^{\gamma+1}(\bm{r})d^{3}r,$
(7)

where $\alpha=-124.3$ MeV, $\beta=70.5$ MeV, and $\gamma=2$, which can be
obtained by the fits to the ground state properties of finite nuclei.

A phenomenological repulsive Pauli potential is introduced to prevent nucleons
with the same spin-$S$ and isospin-$I$ to come close to each other in the
phase space, which can be presented as

$H_{Pauli}=\frac{c_{P}}{2}\sum_{i}(f_{i}-f_{0})^{\mu}\theta(f_{i}-f_{0}),$ (8)
$f_{i}\equiv\sum_{j}\delta(S_{i},S_{j})\delta(I_{i},I_{j})\left|\left\langle\phi_{i}|\phi_{j}\right\rangle\right|^{2},$
(9)

where $f_{i}$ is the overlap of the $i$-th nucleon with other nucleons having
the same spin and isospin, $\theta$ is the unit step function, and $c_{P}=15$
MeV is a coefficient denoting strength of the Pauli potential. For the other
two parameters, we take $f_{0}=1.0$ and $\mu=1.3$. And the symmetry potential
can be written as

$H_{Symmetry}=\frac{C_{S}}{2\rho_{0}}\sum_{i,j\neq{i}}\int[2\delta(I_{i},I_{j})-1]\rho_{i}(\bm{r})\rho_{j}(\bm{r})d^{3}r,$
(10)

where $C_{S}$ is the symmetry energy coefficient which is 25 MeV in this work.

###  II.2 Liénard-Wiechert potential

Liénard and Wiechert proposed that the electromagnetic potential of a moving
particle can be obtained from a scalar potential field $\phi$ and a vector
potential field $\bm{A}$ by Lorentz gauge respectively as follows [64]

$\phi(\bm{r},t)=\frac{1}{4\pi\epsilon_{0}}\left(\frac{q}{(1-\bm{n}\cdot\bm{v}_{s}/c)|\bm{r}-\bm{r}_{s}|}\right)_{t_{s}}\,,$
(11)
$\bm{A}(\bm{r},t)=\frac{\mu_{0}}{4\pi}\left(\frac{q\bm{v}_{s}}{(1-\bm{n}\cdot\bm{v}_{s}/c)|\bm{r}-\bm{r}_{s}|}\right)_{t_{s}}\,,$
(12)

where $\epsilon_{0}$ and $\mu_{0}$ are the dielectric constant and
permeability in vacuum, the speed of light is $1/\sqrt{\epsilon_{0}\mu_{0}}$,
and $q$ is the charge of the particle. $\bm{n}$ is the unit vector represented
by $(\bm{r}-\bm{r}_{s})/|\bm{r}-\bm{r}_{s}|$, and $\bm{v}_{s}$ and
$\bm{r}_{s}$ are the velocity and position of the particle at time $t_{s}$,
respectively. The potential field $\phi(\bm{r},t)$ and $\bm{A}(\bm{r},t)$ at
$t=t_{s}+|\bm{r}-\bm{r}_{s}|/c$ which is excited by particle $q$ at $t_{s}$ is
called the Liénard-Wiechert potential and often referred as the retarded
potential.

The electromagnetic field generated by moving charged particles can be
expressed as

$\displaystyle e\bm{E}(\bm{r},t)$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{4\pi\epsilon_{0}}\sum_{n}\frac{Z_{n}}{(cR_{n}-\bm{R}_{n}\cdot\bm{v}_{n})^{3}}\\{(c^{2}-v_{n}^{2})(\bm{R}_{n}-R_{n}\bm{v}_{n})+(\bm{R}_{n}\cdot\bm{a}_{n})(c^{2}\bm{R}_{n}-cR_{n}\bm{v}_{n})$
(13) $\displaystyle-$
$\displaystyle[c^{2}R_{n}^{2}-cR_{n}(\bm{R}_{n}\cdot\bm{v_{n}})]\bm{a}_{n}\\}\,,$
$\displaystyle e\bm{B}(\bm{r},t)$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{4\pi\epsilon_{0}c}\sum_{n}\frac{Z_{n}}{R_{n}(cR_{n}-\bm{R}_{n}\cdot\bm{v}_{n})^{3}}\\{R_{n}(c^{2}-v_{n}^{2})(\bm{v}_{n}\times\bm{R}_{n})+cR_{n}(\bm{R}_{n}\cdot\bm{a}_{n})(\bm{v}_{n}\times\bm{R}_{n})$
(14) $\displaystyle+$
$\displaystyle\bm{R}_{n}(c\bm{R}_{n}-cR_{n}\bm{v}_{n})(\bm{a}_{n}\times\bm{R}_{n})\\}\,,$

where $Z_{n}$ and $\bm{a}_{n}$ is the number of charges and acceleration on
the particle $n$, respectively; and $\bm{R}_{n}=\bm{r}-\bm{r}_{sn}$, where
$\bm{r}_{sn}$ and $\bm{v}_{n}$ are the position and velocity of the particle
$n$ at retarded time $t_{sn}=t-|\bm{r}-\bm{r}_{sn}(t_{sn})|/c$, respectively.
In particular, the charge constant $q$ is multiplied on the Eq. (15) and the
Eq. (16), so that the right sides of them contain the fine structure constant
$\alpha=e^{2}/(4\pi)=1/137$ $(\epsilon_{0}=\hbar=c=1)$. If we ignore the
acceleration of charged particles, the electromagnetic field can be simplified
as [9, 3]

$e\bm{E}(\bm{r},t)=\frac{e^{2}}{4\pi\epsilon_{0}}\sum_{n}Z_{n}\frac{c^{2}-v_{n}^{2}}{(cR_{n}-\bm{R}_{n}\cdot\bm{v}_{n})^{3}}(c\bm{R}_{n}-R_{n}\bm{v}_{n})\,,$
(15)
$e\bm{B}(\bm{r},t)=\frac{e^{2}}{4\pi\epsilon_{0}c}\sum_{n}Z_{n}\frac{c^{2}-v_{n}^{2}}{(cR_{n}-\bm{R}_{n}\cdot\bm{v}_{n})^{3}}\bm{v}_{n}\times\bm{R}_{n}\,,$
(16)

In the non-relativistic approximation, that is, $v\ll{c}$, the Eq. (15) and
the Eq. (16) can be further simplified as

$e\bm{E}(\bm{r},t)=\frac{e^{2}}{4\pi\epsilon_{0}}\sum_{n}Z_{n}\frac{\bm{e}_{n}}{R_{n}^{3}}\,,$
(17)
$e\bm{B}(\bm{r},t)=\frac{e^{2}}{4\pi\epsilon_{0}c}\sum_{n}Z_{n}\frac{\bm{v}_{n}\times\bm{R}_{n}}{R_{n}^{3}}\,.$
(18)

For self-consistency, we use the original Coulomb potential of the EQMD model
for the electric field part, expressed as

$H_{Coulomb}=\frac{e^{2}}{2}\sum_{i}\sum_{i\neq{j}}Z_{i}Z_{j}\frac{1}{r_{ij}}{\rm
erf}(\frac{r_{ij}}{\sqrt{4L}})\,,$ (19)

where

$r_{ij}=|\bm{r}_{i}-\bm{r}_{j}|\,,$ (20) ${\rm
erf}(x)=\frac{2}{\sqrt{\pi}}\int^{x}_{0}e^{-u^{2}}du\,,$ (21)

and the magnetic field can be obtained from the Eq. (14) or Eq. (16) for
making a comparison.

With the Eq. (16), we displayed time evolution of the magnetic field at
different beam energies. For an example, Fig. 1(a) shows time evolution of the
magnetic field at the central point (0,0,0) produced in 16O + 40Ca collision
process at 100 MeV/nucleon and impact parameter of 7.13 fm. One can see that
the maximum strength of the magnetic field reaches $2.75\times 10^{15}$ Gauss
(in Lorentz units, 1 MeV${}^{2}\approx$ 5.11 $\times$ 1013 Gauss) at $\sim$ 60
$\rm fm/c$ for 100 MeV/nucleon collision and the lifetime of the magnetic
field can extend to $\sim$ 120 $\rm fm/c$. By integrating over the time
evolution of the magnetic field, we can get the magnetic field integral as a
function of beam energy as shown in Fig. 1(b), from which we observed an
almost linear correlation.

Figure 1: Time evolution of the magnetic field at the central point (0,0,0) of
16O + 40Ca at beam energy of 100 MeV/nucleon and impact parameter of 7.13 fm
(a). The dependence between magnetic field integral over time and beam energy
for 16O + 40Ca at 7.13 fm (b).

###  II.3 GDR algorithm

Macroscopic description of IVGDR by the Goldhaber-Teller model [65] is applied
to calculate GDR [66, 67] from the nuclear phase space obtained by the EQMD
model. There are two ways to calculate IVGDR spectrum, the one is to give the
nucleus a boost to obtain the dipole oscillation, the another is to simulate a
Coulomb excitation with a nucleus [54]. In this work, we chose the latter one.

We get the initial state wave function $\Psi(0)$ of nucleus by the EQMD model
and let 16O interact with the target nucleus 40Ca at some energies and an
impact parameter $7.13$ fm, which is the sum of radii of 40Ca and 16O.
Besides, the process to obtain the evolution of excited wave function to the
final state is described in Ref. [52] in details.

The dipole moments of the system in the coordinate space $D_{G}(t)$ and
momentum space $K_{G}(t)$ can be written as [53, 54, 36, 58, 68]

$D_{G}(t)=\frac{NZ}{A}[R_{Z}(t)-R_{N}(t)]\,,$ (22)
$K_{G}(t)=\frac{NZ}{A\hbar}\left[\frac{P_{Z}(t)}{Z}-\frac{P_{N}(t)}{N}\right]\,,$
(23)

where $R_{Z}(t)[P_{Z}(t)]$ and $R_{N}(t)[P_{N}(t)]$ are the center-of-mass of
the protons and neutrons in the coordinate (momentum) space, respectively, $A$
is the sum of mass numbers of the target and projectile nuclei, $N$ and $Z$
are the neutron and proton numbers of the compound system, respectively.
$K_{G}(t)$ is the canonically conjugate momentum of $D_{G}(t)$ [53].

Figure 2: Time evolution of the dipole moments (a) and $\gamma$-emission
probability (b) for 40Ca excited by 16O at incident energy $100$ MeV/nucleon.

The $\gamma$-emission probability of nuclear system at photon energy
$E_{\gamma}=\hbar\omega$ can be derived from Eq. (22) as

$\frac{dP}{dE_{\gamma}}=\frac{2e^{2}}{3\pi\hbar
c^{3}E_{\gamma}}|D^{\prime\prime}(\omega)|^{2}\,,$ (24)

where $D^{\prime\prime}(\omega)$ is obtained from the Fourier transform of the
second derivative of $D_{G}(t)$ with respect to time, i.e.,

$D^{\prime\prime}(\omega)=\int^{t_{2}}_{t_{1}}D^{\prime\prime}_{G}(t)e^{i\omega
t}dt\,.$ (25)

Considering the lifetime of GDR excitation, we take the final state at
$t_{2}=200\,\text{fm}/c$, which corresponds with the lifetime of GDR
excitation.

##  III Results and discussion

Figure 3: The dependence of peak energy (a), strength (b) and width of GDR
spectrum (c) on incident energy for 40Ca excited by 16O in the cases whether
or not to consider acceleration. The symbol denoted “a” represents the results
with particle acceleration effect, while “neglect a” represents the ones
without. The lines are fits to the dots. The hatched area represents the
experimental result. Figure 4: The dependence of peak energy (a), strength (b)
and width of GDR spectrum (c) on incident energy for 40Ca excited by 16O in
the cases without or with magnetic field.

###  III.1 Impacts on GDR of 40Ca

We discuss ultra-peripheral collision case (here it means b = $7.13$ fm) of
40Ca and 16O, and focus on the GDR analysis for excited 40Ca. Many previous
studies [53, 54] have shown the existence of $\alpha$-cluster structure in
16O. For simplicity, we only consider 16O whose nucleons are distributed by
the Woods-Saxon distribution without any $\alpha$-clusters. The
electromagnetic fields come from the 40Ca and 16O interactions. In our
simulations, more than 90% of 40Ca remain complete, and the following results
have excluded the events that 40Ca undergoes fragmentation. First of all, we
check the dipole moment $D_{G}(t)$ and $\gamma$-emission probability for 40Ca
with definitions in Eqs. (22) and (24) which are shown in Fig. 2. Seeing from
Fig. 2 (a), dipole moment $D_{G}(t)$ seems not to be affected by the magnetic
field. However, one can see that $\gamma$-emission probability is slightly
affected by the magnetic field and the peak value increases a little bit.

Figure 5:  Temperature (a) and the $y$-component of angular momentum (b) as a
function of incident energy for 40Ca excited by 16O in cases of without and
with magnetic field.

By Gaussian fitting to the GDR spectrum, the peak energy $E_{\gamma}^{c}$, the
strength $(dP/dE_{\gamma})^{c}$, and the width $\Gamma_{\gamma}^{c}$ can be
quantitatively extracted. In Eq. (14), particle acceleration is considered.
One could notice that accelerated charged particles emit EM waves and lose
their energies. As shown in Fig. 3, one can see that when particle
acceleration is taken into account, energy position and strength of the GDR
decrease and the width of the GR becomes wider. However, such effects are
small. So we neglect particle acceleration effect in the following
calculations.

Fig. 4 (a) shows the comparison between our calculations and the experimental
data from Ref. [69, 70]. The value of $E_{\gamma}^{c}$ is fitted to the data
around 100 MeV/nucleon. As incident energy increases, $E_{\gamma}^{c}$
increases and is higher than the data when E $>$ 200 MeV/nucleon. And one can
notice that for all of peak energy, strength and width of GDR spectrum, they
are enhanced when the magnetic field is switched on. In addition, one can see
that the width decreases with the increasing of incident energy when the
incident energy is less than 200 MeV/nucleon but keeps stable as incident
energy further increases. Also the gap is different between the cases for more
and for less than 250 MeV/nucleon. Such a different trend could be related to
different effects of the magnetic field generated by the projectile on the
target.

###  III.2 Impacts on temperature and angular momentum of 40Ca

In Fig. 4, we find that the magnetic field enhances the peak energy, strength
and width of GDR spectrum. As we know, here the GDR mode is based on the
excited state of the nucleus. Therefore here we check the excitation energy of
target nucleus 40Ca. Fig. 5(a) displays the incident energy dependence of
temperature which is obtained by $E^{\ast}$ = $\it{k}T^{2}$, where $E^{\ast}$
= $(E_{total}-E_{ground})/A$ is the calculated excitation energy per nucleon
[71]. The level density parameter $k$ is taken as $1/8$ and $A$ represents the
mass number of 40Ca [72]. Since we are just focusing on ultra-peripheral
collision, the higher the incident energy is, the shorter the interaction time
between 40Ca and 16O is, and then the lower the excitation energy of 40Ca is.
Thus, the $T$ of 40Ca decreases with the increasing of incident energy.
However, the influence of magnetic field on the temperature differs from the
results of GDR width as shown in Fig. 4(c). It indicates that something
affected by the magnetic field other than temperature, such as rotation of
nucleus, caused the different trend of the width [39]. Thus we check the
angular momentum of nucleus, which can be calculated by $\bm{J}$ =
$\sum_{i}\bm{r_{i}}\times\bm{p_{i}}$ where $i$ represents each nucleon [73,
74] under the center of mass system of 40Ca, and obtain the results as shown
in Fig. 5(b). One can see that the dependence of $y$-component of angular
momentum $J_{y}$ on incident energy with magnetic field is higher than that
one without magnetic field. As in Refs. [39, 75], the width of GDR increases
with angular momentum of nucleus. So here the enhancement of $J_{y}$ by
magnetic field could be related to the results in Fig. 4(c). Also, combining
Fig. 5(a) and Fig. 5(b), we found that the magnetic field reduces $T$ and
increases $J_{y}$ below 250 MeV/nucleon. While above 250 MeV/nucleon, the
magnetic field increases both of $T$ and $J_{y}$ that could explain why the
enhancement of width of GDR becomes larger above 250 MeV/nucleon as shown in
Fig. 4. One can also notice that for the low beam energy, the angular momentum
is large even without magnetic field. So in low energy region, it may be not
only vibrational mode for GDR but also mixed modes due to rotation that
determine the width of the GDR. In this case, it is hard to extract a pure
vibrational mode for GDR. However, as beam energy increases, the angular
momentum becomes small which would reduce such a mixed mode effect. Thus it
could be the reason for the fact that at higher beam energy region, the energy
position of the “GDR” is close to the experimental result. Note that the
experimental data corresponds to the case with only the electric field.

Figure 6: The dependences of peak energy (a), strength (b) and width of GDR
spectrum (c) on temperature for 40Ca in the cases without or with magnetic
field.

###  III.3 Correlation between the GDR and temperature

Considering that the GDR with finite temperature is an interesting issue [75],
we plot temperature dependence of GDR variables in Fig. 6. It is clear to see
that for $E_{\gamma}^{c}$, it decreases as the temperature increases.
Actually, in analogy with a classical oscillator with a friction force, the
higher the temperature, the greater the probability of two-body dissipation,
which means that the frequency of collective motion is reduced [76, 31], thus
the $E_{\gamma}^{c}$ is decreased. Certainly, for $\Gamma_{\gamma}^{c}$, it
increases with temperature, which is consistent with the results in Refs. [77,
78]. Additionally, the enhancement of two-body dissipation reduces the
available phase space of the $\Gamma_{\gamma}^{c}$ resulting in the decrease
of the $(dP/dE_{\gamma})^{c}$ at higher $T$ [33]. For the magnetic field
dependence at different temperatures, we can refer to Fig. 5(a). Moreover, the
magnetic field does not only affect temperature but also angular momentum of
nucleus as shown in Fig. 5(b). Thus the magnetic effect on GDR indicates a
mechanism of mixture of temperature and angular momentum effects.

##  IV Conclusion

In summary, by taking into account the magnetic field in the extended quantum
molecular dynamics model (EQMD), we analyzed its effects on GDR of 40Ca
through the Coulomb excitation process of 16O + 40Ca collision. We found that
the magnetic field enhances the peak energy, strength and width of GDR
spectrum. By checking the effects on temperature and angular momentum, we
found that magnetic field increases the angular momentum of 40Ca. Thus how the
magnetic field acts on angular momentum and nucleus shape evolution is an
interesting issue and should be explored in future study. As for the
correlation between the GDR and temperature for excited 40Ca under electric
and magnetic fields, the peak energy and strength decreases with the
increasing of temperature, but the width enhances, which is in overall
agreement with other theoretical work. The current work sheds light on
resurveying the magnetic field effects on nuclear collective excitation in
low-intermediate energy heavy-ion collisions.

This work was supported in part by the National Natural Science Foundation of
China under contract Nos. 11890710, 11890714, and 12147101, and the Guangdong
Major Project of Basic and Applied Basic Research No. 2020B0301030008.

## References

  * [1] J. Rafelski and B. Müller, Phys. Rev. Lett. 36, 517 (1976).
  * [2] V. Skokov, A. Y. Illarionov, V. Toneev, Int. J. Mod. Phys. A 24, 5925 (2009).
  * [3] V. Voronyuk, V. D. Toneev, W. Cassing et al. Phys. Rev. C 83, 054911 (2011).
  * [4] K. Fukushima, D. E. Kharzeev, H. J. Warringa, Phys. Rev. D 78, 074033 (2008).
  * [5] D. E. Kharzeev, L. D. McLerran, H. J. Warringa, Nucl. Phys. A 803, 227 (2008).
  * [6] M. Asakawa, A. Majumder, B. Müller, Phys. Rev. C 81, 064912 (2010).
  * [7] X. L. Zhao, G. L. Ma, Y. G. Ma, Phys. Rev. C 99, 034903 (2019).
  * [8] K. Hattori and X. G. Huang, Nucl. Sci. Tech. 28, 17 (2017).
  * [9] L. Ou and B. A. Li, Phys. Rev. C 84, 064605 (2011).
  * [10] X. G. Deng and Y. G. Ma, Eur. Phys. J. A 54, 204 (2018).
  * [11] G. F. Wei, B. A. Li, G. C. Yong, L. Ou, X. W. Cao, X. Y. Liu, Phys. Rev. C 97, 034620 (2018).
  * [12] Y. L. Sun, Y. J. Wang, Q. F. Li, F. Q. Wang, Phys. Rev. C 99, 064607 (2019).
  * [13] X. G. Deng and Y. G. Ma, Phys. Rev. C 101, 054610 (2020).
  * [14] M. N. Harakeh and A. Van der Woude, Giant Resonances: Fundamental High-Frequency Modes of Nuclear Excitation (Oxford University Press, 2001).
  * [15] N. Paar, D. Vretenar, E. Khan, G. Coló, Rep. Prog. Phys. 70, R02 (2007).
  * [16] Ait Ben Mennana and M. Oulne, Eur. Phys. J. Plus 136, 85 (2021).
  * [17] J. H. Bai, Z. M. Niu, B. Y. Sun, Y. F. Niu, Phys. Lett. B 815, 136147 (2021).
  * [18] A. Bracco, E. G. Lanza, A. Tamii, Prog. Part. Nucl. Phys. 106, 360-433 (2019).
  * [19] L. Trippa, G. Coló, and E. Vigezzi, Phys. Rev. C77, 061304(R) (2008).
  * [20] J. Piekarewicz, B. K. Agrawa, G. Coló, W. Nazarewicz, N. Paar, P.-G. Reinhard, X. Roca-Maza, and D. Vretenar, Phys. Rev. C 85, 041302(R) (2012).
  * [21] A. Tamii, I. Poltoratska, P. von Neumann-Cosel et al., Phys. Rev. Lett. 107, 062502 (2011).
  * [22] Jun Xu, Chin. Phys. Lett. 38, 042101 (2021).
  * [23] P. Carlos, H. Beil, R. Bergere, A. Leprétre, A. Veyssiére, Nucl. Phys. A 172, 437 (1971).
  * [24] P. Carlos, H. Beil, R.Bergére, A.Leprétre, A.De Miniac, A. Veyssiére, Nucl. Phys. A 225, 171 (1974).
  * [25] B. L. Berman, S. C. Fultz, Rev. Mod. Phys. 47, 713 (1975).
  * [26] V. M. Masur and L. M. Mel’nikova, Phys. Part. Nucl. 37, 923 (2006).
  * [27] L. Donaldson, C. Bertulani, J. Carter et al., Phys. Lett. B 776, 133 (2018).
  * [28] E. Ramakrishnan, T. Baumann, A. Azhari et al., Phys. Rev. Lett. 76, 2025 (1996).
  * [29] C. A. Bertulani, S. R. Klein, J. Nystrand, Annu. Rev. Nucl. Part. Sci. 55, 271 (2005).
  * [30] N. Auerbach and S. Shlomo, Phys. Rev. Lett. 103, 172501 (2009).
  * [31] C. Q. Guo, Y. G. Ma, W. B. He et al., Phys. Rev. C 95, 054622 (2017).
  * [32] D. Mondal, D. Pandit, S. Mukhopadhyay et al., Phys. Rev. Lett. 118, 192501 (2017).
  * [33] S. Bhattacharya D. Pandit, B. Dey et al., Phys. Rev. C 103, 014305 (2021).
  * [34] N. Loiudice, F. Palumbo, Phys. Rev. Lett. 41, 1532 (1978).
  * [35] D. Pandit, B. Dey, D. Mondal et al., Phys. Rev. C 87, 044325 (2013).
  * [36] S. S. Wang, Y. G. Ma, X. G. Cao, W. B. He, H. Y. Kong, C. W. Ma, Phys. Rev. C 95, 054615 (2017).
  * [37] B. Dey, S. S. Wang, D. Pandit, S. Bhattacharya, X. G. Cao, W. B. He, Y. G. Ma, N. Quang Hung, and N. Dinh Dang, Phys. Rev. C 102, 031301(R) (2020).
  * [38] P. Axel, Phys. Rev. C 126, 671 (1962).
  * [39] Y. Alhassid, Nucl. Phys. A 649, 107-114 (1999).
  * [40] B. A. Li, L. W. Chen, C. M. Ko, Phys. Rep. 464, 113 (2008).
  * [41] J. Liu, C. Gao, N. Wan, C. Xu, Nucl. Sci. Tech. 32, 117 (2021).
  * [42] W. Kleinig, V. O. Nesterenko, J. Kvasil, P.-G. Reinhard, and P. Vesely, Phys. Rev. C 78, 044313 (2008).
  * [43] V. Nesterenko, W. Kleinig, J. Kvasil, P.-G. Reinhard, and P. Vesely, Int. J. Mod. Phys. E 16, 624 (2007).
  * [44] J. A. Maruhn, P. G. Reinhard, P. D. Stevenson, J. Rikovska Stone, and M. R. Strayer, Phys. Rev. C 71, 064328 (2005).
  * [45] S. Fracasso, E. B. Suckling, P. Stevenson, Phys. Rev. C 86, 044303 (2012).
  * [46] D. Sarchi, P. F. Bortignon, G. Coló, Phys. Lett. B 601, 27 (2004).
  * [47] D. Vretenar, A.V. Afanasjev, G.A. Lalazissis, P. Ring, Phys. Rep. 409, 101 (2005)..
  * [48] D. Pena Arteaga, E. Khan, P. Ring, Phys. Rev. C 79, 034311 (2009).
  * [49] Rui Wang, Zhen Zhang, Lie-Wen Chen, Che Ming Ko, Yu-Gang Ma, Phys. Lett. B 807, 135532 (2020).
  * [50] Rui Wang, Zhen Zhang, Lie-Wen Chen, Yu-Gang Ma, Frontiers in Physics. 8, 330 (2020).
  * [51] Yi-Dan Song, Rui Wang, Zhen Zhang, Yu-Gang Ma, Phys. Rev. C 104, 044603 (2021).
  * [52] T. Maruyama, K. Niita, A. Iwamoto, Phys. Rev. C 53, 297 (1996).
  * [53] W. B. He, Y. G. Ma, X. G. Cao, X. Z. Cai, G. Q. Zhang, Phys. Rev. Lett. 113, 032506 (2014).
  * [54] W. B. He, Y. G. Ma, X. G. Cao, X. Z. Cai, and G. Q. Zhang, Phys. Rev. C 94, 014301 (2016).
  * [55] J. Aichelin, Phys. Rep. 202, 233 (1991).
  * [56] T. Furuta, K. H. O. Hasnaoui, F. Gulminelli, C. Leclercq, and A. Ono, Phys. Rev. C 82, 034307 (2010).
  * [57] Y. Kanada-En’yo and M. Kimura, Phys. Rev. C 72, 064301 (2005).
  * [58] H. L. Wu, W. D. Tian, Y. G. Ma et al. Phys. Rev. C 81, 047602 (2010).
  * [59] B. S. Huang and Y. G. Ma, Phys. Rev. C 103, 054318 (2021).
  * [60] J. Aichelin, H. Stöcker, Phys. Lett. B 176, 14 (1986).
  * [61] C. Hartnack, R. K. Puri, J. Aichelin, J. Konopka, S. A. Bass, H. Stöcker, W. Greiner Eur. Phys. J. A 1, 151 (1998).
  * [62] A. Kerman and S. Koonin, Ann. Phys. 100, 332 (1976).
  * [63] C. Z. Shi and Y. G. Ma, Nucl. Sci. Tech. 32, 66 (2021).
  * [64] David J. Griffiths, Introduction to electrodynamics (Prentice Hall, 3rd edition, January 1999).
  * [65] M. Goldhaber and E. Teller, Phys. Rev. 74, 1046 (1948).
  * [66] S. Bacca, N. Barnea, G. Hagen, G. Orlandini, and T. Papenbrock, Phys. Rev. Lett. 111, 122502 (2013).
  * [67] S. Bacca, N. Barnea, G. Hagen, M. Miorelli, G. Orlandini, and T. Papenbrock, Phys. Rev. C 90, 064619 (2014).
  * [68] V. Baran, M. Cabibbo, M. Colonna, M. Di Toro, N. Tsonev, Nucl. Phys. A679, 373 (2001).
  * [69] J. Ahrens, H. Borchert, K. H. Czocktt et al., Nucl. Phys. A 251, 479 (1975).
  * [70] V. A. Erokhova, M. A. Yolkin et al. Bull. Russ. Acad. Sci. Phys. 67, 1636 (2003).
  * [71] Y. G. Ma, A. Siwek and J. Péter et al., Phys. Lett. B 390, 41 (1997.
  * [72] W. Reisdorf, Z. Phys. A 300, 227 (1981).
  * [73] Y. G. Ma, W. Q. Shen, Z. Y. Zhu, Phys. Rev. C 51, 1029 (1995).
  * [74] Zhi-Wan Xu, Song Zhang, Yu-Gang Ma et al. Nucl. Sci. Tech. 29, 186 (2018).
  * [75] D. Santonocito and Y. Blumenfeld, Eur. Phys. J. A 56:279 (2020)
  * [76] M. Di Toro, V. M. Kolomietz, and A. B. Larionov, Phys. Rev. C 59, 3099 (1999).
  * [77] N. D. Dang, J. Phys. G: Nucl. Part. Phys. 31, S1697 (2005).
  * [78] D. Mondal, D. Pandit, S. Mukhopadhyay et al., Phys. Lett. B 784, 423 (2018).


# A Double-Spring Model for Nanoparticle Diffusion in a Polymer Network

Yu Lu  School of Mechanics and Engineering Science, Shanghai Institute of
Applied Mathematics and Mechanics, Shanghai Key Laboratory of Mechanics in
Energy Engineering, Shanghai University, Shanghai, China  Xin-Yue Liu  School
of Mechanics and Engineering Science, Shanghai Institute of Applied
Mathematics and Mechanics, Shanghai Key Laboratory of Mechanics in Energy
Engineering, Shanghai University, Shanghai, China  Guo-Hui Hu
<EMAIL_ADDRESS>School of Mechanics and Engineering Science, Shanghai
Institute of Applied Mathematics and Mechanics, Shanghai Key Laboratory of
Mechanics in Energy Engineering, Shanghai University, Shanghai, China

###### Abstract

The transport of nanoparticles (NPs) in polymer networks, as a typical
simplified model describing various structures in living systems, is
profoundly important in biomedical engineering and nanotechnology. Predicting
the effective diffusivity of NP confined in an ordered network has been an
intriguing focus in this frontier field. In the present study, the diffusion
of NPs in an unentangled polymer network for different NP radii and network
stiffness is numerically investigated by single particle dissipative particle
dynamics (DPD). It is found that, the deformation due to the junction
deviation contributes significantly to the the potential barrier $U$ for the
NP to overcome during hopping, and it is dominated over the strain energy
induced by loop stretching for larger NPs and lower network rigidity. Analyses
based on the theory of continuum mechanics reveal that the relation between
this deformed energy and the junction deviation can be described by a non-
linear spring. Taking into account both effects of the loop stretching and
junction deviation, a double-spring model is proposed to characterize the
diffusivity of the NPs in the ordered network. The theoretical prediction is
in good agreement with our numerical simulations, and qualitatively consistent
with the investigations available. This model is helpful to improve our
understanding on the dynamic behavior of nanoparticle in complex biological
environment, and provide theoretical guidance in designing biomedical
applications.

††preprint: APS/123-QED

The diffusion of nanoparticles (NP) in heterogeneous medium, such as polymer
liquids [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], porous medium [11, 12, 13, 14, 15,
16, 17, 18, 19] is of great significance in targeting delivery of
nanomedicines, various fundamental biological processes and plenty of physical
chemical applications [20, 21, 22, 23, 24, 25, 26, 27]. These heterogeneous
media are usually described as local confinements in theoretical models[10, 9,
19, 28], leading to the diffusion properties distinct from the classic
Brownian motion. Previous investigations [19, 29] indicated that at short
timescales the particles show a transition from ballistic to Brownian motion,
while at intermediate time the NPs’ motion is restricted by the surrounding
structures, manifesting a sublinear mean square displacement (MSD). At long
timescales, the NP in complicated medium behaviors like Brownian motion with
an effective diffusivity $D^{L}$.

Figure 1:  (a) A schematic depiction of the ordered network constructed by
polymer strands and NP (blue). (b) The average distance between two junctions
(cyan), which is connected by strand (yellow), is $a_{x}$. A double spring
model (e) is proposed to simultaneously characterize the deformation of the
polymer network and the potential barriers overcome by the nanoparticle during
hopping. The first spring measures the loop stretching (c) in the mesh of NP
occupied, which can be simplified by a linear spring with stiffness
$k_{l}\simeq\chi\beta^{-1}$. The junction deviation (d) when NP tries to
squeeze in the polymer loop will induce another kind of deformation in polymer
network. By using the continuum mechanics, the deformation energy of junction
deviation can be represented by a sublinear spring with stiffness
$k_{n}\simeq\beta^{-2}$.

The motion of a NP in a heterogeneous medium with strong local confinement,
such as permanently cross-linked polymer network, has been described by the
nonlinear Langevin equation (NLE), and the local constraints can be considered
as a potential barrier $U$ [9, 10, 28, 19]. Activated hopping of nanoparticles
can be observed [25, 28, 9, 1, 10, 19] when it diffuses by overcoming the free
energy barrier. The mechanism of NP hopping in polymer networks is worth
investigating in depth, which might be closely related to the physical origin
of the potential barrier. However, the theories available[28, 9] presented
different scaling law of the potential barrier dependent on the NP radius. Cai
et al.[28] suggested the mechanism underlying is the loop stretching in
polymer networks, and proposed a scaling laws based on the “elementary”
network theory under the assumption that the restoring force of the loops is
linearly elastic with the stiffness proportional to the NP size, and the
deformed energy can be written in a parabolic form. Dell and Schweizer[9]
developed a theory of hopping based on a nonlinear Langevin equation and
polymer reference interaction site model (PRISM) [30], showing that the
potential barrier is linearly increase with NP size, following by a cubical
increase. Utilizing molecular dynamics simulations, Sorichetti et al. [31]
reported that both theories present satisfying predictions, and it is hard to
evaluate their reliability. The numerical study conducted by Chen et al.[29]
found that the diffusion behaviors of NP could be qualitatively different as
the length of polymer chains in the network varies. Using the strand
fluctuation distance $d_{f}$ as a normalized length scale, they found that the
polymer networks with shorter polymer strands (lower $d_{f}$) perform higher
local confinement on NPs, which corresponds to larger potential barrier.

Although extensive investigations have been conducted to explore the NP
diffusion in polymer networks, a significant question remaining open is to
predict the diffusive behavior theoretically. An effective way to solve this
problem might be to formulate an accurate mathematical expression for the
potential barrier in nonlinear Langevin equation, which requires an in-depth
understanding of the deformation of the polymer network when the nanoparticle
traverses its loop. As indicated by previous literatures [28, 19], the
potential barrier needed to be overcome by the NPs is related to the
deformation energy of the polymer network as the NPs hopping. Thus two crucial
questions have to be examined. The first one is to quantitatively understand
the physics behind the polymer network deformation during NP squeezing the
loop. This goal can be achieved by analyzing the numerical information based
on mesoscopic simulations, whereas would be difficult for experimental
measurements. Secondly, how the properties of the polymer network, especially
mechanical properties, influence the local confinement on NP? The theories
developed by James & Guth[32] and Flory[33] are helpful to connect these two
questions. In their theory, the positions of junctions in polymer network is
in the form of Gaussian probability density function, and their standard
deviations are defined as the fluctuation length $d_{f}$. Then the elasticity
modulus $\nu$ of polymer network has magnitude of $\nu\sim d_{f}^{-2}$[34]. In
our recent numerical study[19], we carried out a detailed analysis on the
trajectories of NP and the junctions in polymer network. It is found that the
junctions in polymer network fluctuate around their origin locations with a
non-zero mean value when NP tries to hop in the polymer loop, suggesting that
other than the loop stretching, the deformation energy due to the junction
deviation during hopping contributes significantly to the potential barrier.

The single-particle Dissipative Particle Dynamics (DPD)[35, 36, 37] is used to
simulate the diffusion of NP in an unentangled polymer network in the present
study. As shown in Fig. 1(a)(b), the network with cubic topology is
constructed by the polymer strands modelled by the classical bead-spring
model. A single strand in the polymer network is consisted of $N$ DPD
particles, which represent the Kuhn monomer with size $b$. The polymer strands
in the network are permanently connected by the junctions. The average
distance between two neighbouring junctions is the mesh size $a_{x}$, which is
given by $a_{x}\simeq bN^{\mu}$, in which $\mu$ is the Flory exponent
depending on the solvent properties. Usually $\mu=0.588$ for athermal and good
solvent. The system is filled by the solvent molecules represented by single
DPD particles, and their number density $n$ is kept to be a constant. The
volume fraction of polymer network can be estimated by
$\varphi=(3N-5)/na_{x}^{3}$. The correlation length $\xi$ can be calculated by
$\xi\simeq b\varphi^{-\mu/(3\mu-1)}$[38], which yields the tube diameter of
$a_{e}\approx 5\xi$[38]. In present work, $N$ and $\varphi$ are selected to
keep the mesh size $a_{x}$ smaller than tube diameter $a_{e}$, thus the
polymer networks are unentangled. Full details of the numerical techniques are
described in the Supplement Information S1.

Since the potential barrier in nanoparticle diffusion is closely correlated
with the polymer network deformation, a normalized fluctuation length
$\beta=d_{f}/a_{x}$ is introduced aiming at considering the influences of the
stiffness on the deformation. It should be mentioned that there is a slight
difference in definition of fluctuation length in previous literature. Flory
et al.[33] proposed that the local constraints on junctions limit the
fluctuation of strands, which affects the elasticity of polymer network. They
considered $d_{f}$ as the mean square displacement(MSD) of the junctions,
while Chen et al.[29] defined it as the MSD of the central monomer of strands.
It is shown in Supporting Information S2 that the results obtained from these
two definitions are similar for the unentangled polymer network, and the
formulation of Flory et al.[33] is used in the present study.

As shown in Fig. 2, the mean square displacement of NP varying with the
normalized NP size $\chi=d_{NP}/a_{x}$ is depicted against time in log-log
coordinates. It is found that, at the short-time stage, the dynamics of NP is
dominated by inertia, the ballistic motion can be observed by $MSD\simeq
t^{2}$ when the NP has not yet collided with polymer strands. At the
intermediate stage, the local confinement constructed by polymer network
limits the diffusion of NPs, and the MSD is sublinear with elapsed time. In
the long-time stage, the NP walks randomly cross the polymer network and MSD
becomes linear, which gradually returns a long-time diffusivity $D^{L}$.

\begin{overpic}[height=113.81102pt]{MSD_vs_axae_cp1.1.eps}
\put(84.0,42.0){(a)} \end{overpic}
\begin{overpic}[height=113.81102pt]{MSD_vs_axae_cp1.3.eps}
\put(85.0,42.0){(b)} \end{overpic}
\begin{overpic}[height=113.81102pt]{U_DL_vs_cp.eps} \put(85.0,40.0){(c)}
\end{overpic}

Figure 2:  MSD for smaller NP size $\chi=1.1$(a) and larger NP size
$\chi=1.3$(b) varied by normalized fluctuation length $\beta$ is plotted
against time in log-log coordinate. (c) The potential barrier (right, dash
line) for NP to overcome and the corresponding long time diffusivity (left,
solid line) varied by normalized fluctuation length $\beta$ are depicted
against the normalized NP size $\chi$.

The NP diffusion in polymer network exhibits distinct dynamic behavior as the
fluctuation length $\beta$ differs. For smaller NP with $\chi=1.1$ (Fig
.2(a)), although the discrepancy of nanoparticle MSD between different $\beta$
is negligible for both intermediate- and long-time stage[29], it is slightly
greater for larger fluctuation length $\beta$. For the larger NP with
$\chi=1.3$ (Fig .2(b)), smaller $\beta$, corresponding to larger stiffness of
the polymer network according to Flory’s theory[33], leads to stronger local
confinement on NP, which can be observed by a plateau in the MSD for
intermediate-time stage. This implies the occurrence of hopping phenomenon,
and results in the rapid reduction of the long time diffusivity $D^{L}$ if the
nanoparticle radii $\chi$ increases, as shown in Fig. 2. The intermediate- and
long-time behaviors of NP are qualitatively different as $\beta$ varies.
Similar results have been reported by Chen et al. recently[29].

Utilizing the Nonlinear Langevin equation (NLE) and Kramer’s theory, it will
yield the effective diffusivity once the potential barrier $U$ is obtained.
For strong local confinement ($U>3$), it has a simplified form given by
$D^{L}\sim e^{-U}$[9, 10, 19]. As demonstrated by previous study[19], the
deformations due to both loop stretching and junction deviation contribute to
the potential barrier $U$ for NP to overcome, as shown in Fig. 1, and the
prediction of the effective diffusivity based on the potential barrier $U$
calculated according to this idea is in good agreement with the simulation
results of DPD. For the influences of the network stiffness on NP diffusion,
it is shown in Fig. 2 that the potential barrier is generally larger for
smaller $\beta$, and consequently the effective diffusivity will be smaller.
Furthermore, as the fluctuation length becoming smaller, $U$ grows more
rapidly with $\chi$. Therefore, both normalized NP radius $\chi$, which is
widely used as the measurement of the confinement by previous study[28, 9,
13], and normalized fluctuation length $\beta$, determine the local
confinement constructed by polymer network on NPs.

\begin{overpic}[height=142.26378pt]{displacement_field.eps}
\put(24.0,50.0){(a)} \end{overpic}
\begin{overpic}[height=142.26378pt]{F32.eps} \put(26.0,50.0){(b)}
\end{overpic} \begin{overpic}[height=113.81102pt]{Wdnu_and_dv.eps}
\put(80.0,40.0){(c)} \end{overpic}
\begin{overpic}[height=113.81102pt]{sublinear_spring.eps} \put(14.0,40.0){(d)}
\end{overpic}

Figure 3: (a) Displacement and (b) corresponding shearing strain of the slice
in network during NP hopping. (c) The variation of the the deformation energy
$W/\nu$ (left, solid line) and the normalized deviation length $d_{v}/a_{x}$
(right, dash line) for different $\chi$ and $\beta$. (d) Deformation energy
$W/\nu$ vs. the normalized deviation length $(\frac{d_{v}}{a_{x}})^{2}$
follows a power law, and the slope is $0.75$ in log-log coordinates.

The potential barrier for NP to overcome during translocation between
neighbouring meshes is physically originated from the deformation energy. To
measure the deformation in the framework of continuum mechanics[32], the
coordinate $\bf{X}$ of a junction in equilibrium locations is used as the
reference states, which can be obtained by the ensemble average of the
junction’s location in the undeformed polymer network. It fluctuates from the
equilibrium location in Gaussian distribution with the variance $d_{f}^{2}$.
The NP hopping will induce the deformation of the polymer network, meanwhile
the junction $\bf{X}$ moves to a place with coordinate $\bf{x(X)}$. The
deformation gradient is calculated by

$F_{iK}({\bf X})=\frac{\partial x_{i}({\bf X})}{\partial X_{K}},$ (1)

in which subscripts $i$ and $K$ are $1,2,3$, representing for three
dimensions. Then the deformation energy normalized by $k_{B}T$ can be
calculated by the well known phantom network model proposed by James and
Guth[32],

$W=\frac{\nu}{2}\int(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}-3)dV$ (2)

in which $dV$ is the volume differential. $\lambda_{i}$ is the eigenvalues of
the deformation gradient tensor, which stands for the stretches of $dV$ in
principal directions. The network stiffness is estimated by
$\nu\sim\beta^{-2}$ for Gaussian chains[34, 39].

Focusing on the slice in the network where the NP passes, its deformation and
the component $F_{32}$ of the deformation gradient are shown in the Fig.
3(a)(b) for $\beta=0.2991$ and $\chi=1.1$. Without loss of generality, we
assume the hopping direction of NP is in $X_{3}$ direction and the NP locates
at the center of the network. Then the junction in the slice is deformed in
the hopping direction with a normalized displacement $(x_{3}-X_{3})/a_{x}$. It
is found that the displacement and its gradient field are in axis-symmetry.
The maximum displacement is denoted as $d_{v}/a_{x}$ for the center junctions
${\bf{X}}_{c}$, which is the deviation length found in our previous study[19].
An important fact is, in ($X_{1}$ \- $X_{2}$) plane perpendicular to the
nanoparticle traversal, the displacement is almost zero even inside the
element occupied by NP. This suggests that the deformed energy due to the
junction deviation is independent of the mechanism of the loop stretching
proposed by Cai et al. [28] It is also found that, the components $F_{23}$ and
$F_{32}$ of deformation gradient tensor $F_{iK}$ are the main contributions in
this deformation energy.

As shown in Fig. 3(c), the variation of both the normalized deviation length
$d_{v}/a_{x}$ and the deformation energy shows similar tendency as $\chi$ and
$\beta$ varying, i.e., they increase for larger nanoparticle and softer
network, or larger $\chi$ and $\beta$. Moreover, the deformation energy $W$
can be related to the deviation length by a power law with an exponent
$\alpha$ independent of the fluctuation length $\beta$,
$W/\nu\sim({d_{v}}/{a_{x}})^{2\alpha}$, as shown in Fig. 3(d). Therefore, the
relation between $W$ and the elongation $d_{v}/a_{x}$ can be simplified as a
non-linear spring:

$W(x_{c})=\frac{k_{n}}{2}(\frac{d_{v}}{a_{x}})^{2\alpha},$ (3)

in which the stiffness of the spring $k_{n}$ is proportional to
$\nu\sim\beta^{-2}$, and the nonlinear parameter $\alpha\approx 0.75$ is
obtained by fitting.

\begin{overpic}[height=113.81102pt]{dv_vs_alkn.eps} \put(80.0,42.0){(a)}
\end{overpic} \begin{overpic}[height=113.81102pt]{U_vs_alkn.eps}
\put(80.0,42.0){(b)} \end{overpic}

Figure 4: The variations of the deviation length $d_{v}/a_{x}$ (a) and the
potential barrier $U/k_{n}$ (b) with the ratio between the stiffness of two
springs, which ensures the theoretical model (dash lines) will obtain
diffusivity quantitatively agreed with the numerical simulations.

These results indicate, there are two kinds of deformations contributing to
the potential barrier during the NP diffusion, which are the loop stretching
proposed by Cai et. al. [28], and the junction deviation due to NP squeezing,
as shown in the Fig. 1(d). In the first kind of deformation, the polymer loop
in the mesh the NP occupies elongates in the $(X_{1}-X_{2})$ plane, and its
length scale changes from $a_{x}$ to $d_{NP}$. The second kind of deformation
describes that, the particle locating at $x_{n0}$ moves to $x_{n}$ in $X_{3}$
direction, resulting in a displacement of the polymer strand by $d_{v}$,
simultaneously the junctions originally locating at $X_{c}$ move to
$x_{c}=X_{c}+d_{v}$. Therefore, the potential barrier during hopping process
can be characterized by a double-springs system (Fig. 1(e)), which is given by

$U(x_{n})=U_{l}+W.$ (4)

The first term $U_{l}$ is the deformation energy of loop stretching normalized
by $k_{B}T$, which has been modelled as a linear spring, and can be written as
(see Appendix C in Ref. [28]),

$U_{l}(x_{n},x_{c})=\frac{k_{l}}{2}(\frac{d_{nc}-d_{nc0}}{a_{x}})^{2},$ (5)

in which $d_{nc}=x_{n}-x_{c}$ is the elongation of this spring is the distance
between NP and junctions. The equilibrium distance $d_{nc0}=x_{n0}-X_{c}$ is
the half of mesh size $d_{nc0}=0.5a_{x}$. $k_{l}$ is the stiffness of the
linear spring, which is considered to be proportional to the confinement
degree of polymer network on the nanoparticle. The ratio $\chi$ between the NP
diameter and the mesh size is a usual choice for confinement parameter.
Alternatively, Chen et al.[29] suggested that the actual local confinement
that NPs can feel is correlated to the fluctuation length of polymer strand
$d_{f}$, rather than mesh size $a_{x}$. So we assume that
$k_{l}\sim\frac{d_{NP}}{d_{f}}=\chi\beta^{-1}$. The detailed derivation of
equation (5) is provided in the Supplement Information (S3). The second term
$W$ in equation (4) is the strain energy induced by the junction deviation,
which is described as the sublinear spring in equation (3). Substituting
equation (3) and (5) into equation (4) yields

$\begin{split}U(x_{n},x_{c})&=U_{l}(x_{n})+W(x_{c})\\\
&=\frac{k_{l}}{2}(\frac{d_{nc}-d_{nc0}}{a_{x}})^{2}+\frac{k_{n}}{2}(\frac{x_{c}-X_{c}}{a_{x}})^{2\alpha}\\\
\end{split}$ (6)

The equilibrium state of the system requires the free energy $U(x_{n},x_{c})$
reaches a minimum, i.e.,

$\frac{\partial U(x_{n},x_{c})}{\partial x_{c}}=0.$ (7)

If we denote the arrival of the nanoparticle at the junction location $x_{c}$
as a successful hopping, i.e., $d_{nc}=x_{n}-x_{c}=0$, and consider the non-
linear parameter $\alpha=0.75$, as shown in Fig. 3(d), the elastic
coefficients of the two springs can be related by the deviation length
$d_{v}/a_{x}$ using equation (7),

$d_{v}/a_{x}=\frac{4}{9}(\frac{k_{l}}{k_{n}})^{2}$ (8)

Then the normalized potential barrier $U/k_{n}$ is

$\frac{U}{k_{n}}=(U_{l}+W)/k_{n}=\frac{1}{8}\frac{k_{l}}{k_{n}}+\frac{2}{9}(\frac{k_{l}}{k_{n}})^{3}$
(9)

in which the parameters $k_{l}=A\beta^{-1}\chi$ and $k_{n}=B\beta^{-2}$ are
determined by the properties of polymer network $\beta=d_{f}/a_{x}$ and
normalized NP size $\chi=d_{NP}/a_{x}$[34, 29]. Based on the potential barrier
we obtained, the coefficients $A$ and $B$ can be estimated by utilizing a
gradient descent algorithm and the least square method, as presented in Fig.
4(b), which yields $A=6.799\pm 0.5637$ and $B=3.9669\pm 0.4367$. To confirm
the effectiveness of the two-spring model, the simulation results of
$d_{v}/a_{x}$ are compared with the theoretical prediction of equation (8) in
Fig. 4(a), which shows well consistence with an error less than $9\%$. This
ensures the theoretical model will obtain diffusivity quantitatively agreed
with the DPD simulation results.

As above mentioned, there are different opinions on the scaling law of the
potential barrier with nanoparticle size in previous studies[28, 9]. Here, we
try to clarify this issue by analyzing the contributions of the two kinds of
deformation energy in our model. In the double-spring model, the deformation
energy of loop stretching $U_{l}\sim\chi\beta^{-1}$, and the strain energy
induced by the junction deviation $W\sim\chi^{3}\beta$. The ratio between the
two parts of potential barrier $r_{U}=\frac{W}{U_{l}}\sim\chi^{2}\beta^{2}$,
which implies that $U_{l}$ will be dominated for smaller NP and more rigid
network. In this case, $U_{l}$ is linear varying with the nanoparticle radius,
which is consistent with the prediction of Cai et al. [28], and inversely
proportional to the stiffness of the polymer network. For larger $r_{U}$, the
junctions deviation becomes more important than the loop stretching on the
deformation energy, thus the potential barrier $U\approx W$ grows cubically
with NP size or linearly with NP volume. This is agreed with the result of
Dell et al.[9], i.e., linear increase of $U$ with NP size followed by a rapid
increase. Moreover, contrary to its effect on $U_{l}$, the rigidity of the
network will enhance the growth of $W$.

In summary, NP mobility in the unentangled polymer network with an ordered
cubic topology is investigated by the coarse-grained molecular dynamics method
for different NP size and network rigidity. The network rigidity is
characterized by the normalized fluctuation length $d_{f}/a_{x}$ of the
polymer strands. Results show that the junction deviation, which is measured
by the averaged displacements of the junctions in the direction of the NP
traversal, have profound contribution to the deformation of the polymer
network, especially for larger nanoparticles and softer networks. This kind of
deformation is analyzed based on the phantom network model[40] in the
conventional framework of continuum mechanics, which indicates that it can be
simplified as a non-linear spring. Therefore, the double-spring model is
presented to predict the effective diffusivity for the NP in polymer network,
combining the mechanism of loop stretching proposed by Cai et al.[28], in
which the related potential barrier is described as a linear spring. The
theoretical prediction has good agreement with our numerical results, and
qualitatively consistent with previous theoretical study by Dell et al.[9].
This model is found to be helpful to clarify the controversy over this issue
in literature available[28, 9]. This theoretical framework is promising to be
deepen our understanding on diffusion of nanoparticle in living systems with
strong heterogeneity.

## References

  * [1] Chundong Xue, Xu Zheng, Kaikai Chen, Yu Tian, and Guoqing Hu.  Probing non-gaussianity in confined diffusion of nanoparticles.  The Journal of Physical Chemistry Letters, 7(3):514–519, 2016.
  * [2] Jagannathan T Kalathi, Umi Yamamoto, Kenneth S Schweizer, Gary S Grest, and Sanat K Kumar.  Nanoparticle diffusion in polymer nanocomposites.  Physical Review Letters, 112(10):108301, 2014.
  * [3] Ting Ge, Gary S. Grest, and Michael Rubinstein.  Nanorheology of entangled polymer melts.  Phys. Rev. Lett., 120:057801, Feb 2018.
  * [4] François Amblard, Anthony C Maggs, Bernard Yurke, Andrew N Pargellis, and Stanislas Leibler.  Subdiffusion and anomalous local viscoelasticity in actin networks.  Physical Review Letters, 77(21):4470, 1996.
  * [5] John van Zanten, Samiul Amin, and Ahmed Abdala.  Brownian motion of colloidal spheres in aqueous peo solutions.  Macromolecules, 37:3874–3880, 07 2004.
  * [6] Yunhong Dong, Xiaoqing Feng, Nanrong Zhao, and Zhonghuai Hou.  Diffusion of nanoparticles in semidilute polymer solutions: A mode-coupling theory study.  The Journal of chemical physics, 143(2):024903, 2015.
  * [7] Umi Yamamoto and Kenneth S Schweizer.  Microscopic theory of the long-time diffusivity and intermediate-time anomalous transport of a nanoparticle in polymer melts.  Macromolecules, 48(1):152–163, 2015.
  * [8] Igor V. Volgin, Sergey V. Larin, Enrique Abad, and Sergey V. Lyulin.  Molecular dynamics simulations of fullerene diffusion in polymer melts.  Macromolecules, 50(5):2207–2218, 2017.
  * [9] Zachary E. Dell and Kenneth S. Schweizer.  Theory of localization and activated hopping of nanoparticles in cross-linked networks and entangled polymer melts.  Macromolecules, 47(1):405–414, 2013.
  * [10] Chundong Xue, Xinghua Shi, Yu Tian, Xu Zheng, and Guoqing Hu.  Diffusion of nanoparticles with activated hopping in crowded polymer solutions.  Nano Letters, 20(5):3895–3904, 2020.
  * [11] Aljaž Godec, Maximilian Bauer, and Ralf Metzler.  Collective dynamics effect transient subdiffusion of inert tracers in flexible gel networks.  New Journal of Physics, 16(9):092002, 2014.
  * [12] Anand Rahalkar and Murugappan Muthukumar.  Diffusion of polyelectrolytes in polyelectrolyte gels.  Macromolecules, 50(20):8158–8168, 2017.
  * [13] Kavindya K Senanayake, Ehsan Akbari Fakhrabadi, Matthew W Liberatore, and Ashis Mukhopadhyay.  Diffusion of nanoparticles in entangled poly (vinyl alcohol) solutions and gels.  Macromolecules, 52(3):787–795, 2019.
  * [14] Andrey G. Cherstvy, Samudrajit Thapa, Caroline E. Wagner, and Ralf Metzler.  Non-gaussian, non-ergodic, and non-fickian diffusion of tracers in mucin hydrogels.  Soft Matter, 15:2526–2551, 2019.
  * [15] Daniel S. Banks, Charmaine Tressler, Robert D. Peters, Felix HöFling, and Cécile Fradin.  Characterizing anomalous diffusion in crowded polymer solutions and gels over five decades in time with variable-lengthscale fluorescence correlation spectroscopy.  Soft Matter, 12(18):4190–4203, 2016.
  * [16] Praveen Kumar, Ligesh Theeyancheri, Subhasish Chaki, and Rajarshi Chakrabarti.  Transport of probe particles in a polymer network: effects of probe size, network rigidity and probe–polymer interaction.  Soft Matter, 15:8992–9002, 2019.
  * [17] Toshihiro Toyota, David A Head, Christoph F Schmidt, and Daisuke Mizuno.  Non-gaussian athermal fluctuations in active gels.  Soft Matter, 7(7):3234–3239, 2011.
  * [18] Johann Hansing, Joseph R. Duke, Emily B. Fryman, Jason E. DeRouchey, and Roland R. Netz.  Particle diffusion in polymeric hydrogels with mixed attractive and repulsive interactions.  Nano Letters, 18(8):5248–5256, 2018.
  * [19] Yu Lu and Guo-Hui Hu.  A potential barrier in the diffusion of nanoparticles in ordered polymer networks.  Soft Matter, pages –, 2021.
  * [20] Francesco Turci and Nigel B. Wilding.  Wetting transition of active brownian particles on a thin membrane.  Phys. Rev. Lett., 127:238002, Dec 2021.
  * [21] Jonathan Garamella, Kathryn Regan, Gina Aguirre, Ryan J. McGorty, and Rae M. Robertson-Anderson.  Anomalous and heterogeneous dna transport in biomimetic cytoskeleton networks.  Soft Matter, 16:6344–6353, 2020.
  * [22] Bo Wang, James Kuo, and Steve Granick.  Bursts of active transport in living cells.  Physical Review Letters, 111(20):208102, 2013.
  * [23] Felix Höfling and Thomas Franosch.  Anomalous transport in the crowded world of biological cells.  Reports on Progress in Physics, 76(4):046602, 2013.
  * [24] Wei He, Hao Song, Yun Su, Ling Geng, Bruce J Ackerson, HB Peng, and Penger Tong.  Dynamic heterogeneity and non-gaussian statistics for acetylcholine receptors on live cell membrane.  Nature Communications, 7:11701, 2016.
  * [25] Bo Wang, Stephen M Anthony, Sung Chul Bae, and Steve Granick.  Anomalous yet brownian.  Proceedings of the National Academy of Sciences, 106(36):15160–15164, 2009.
  * [26] Bo Wang, James Kuo, Sung Chul Bae, and Steve Granick.  When brownian diffusion is not gaussian.  Nature Materials, 11(6):481, 2012.
  * [27] Raffaele Pastore, Antonio Ciarlo, Giuseppe Pesce, Francesco Greco, and Antonio Sasso.  Rapid fickian yet non-gaussian diffusion after subdiffusion.  Phys. Rev. Lett., 126:158003, Apr 2021.
  * [28] Li-Heng Cai, Sergey Panyukov, and Michael Rubinstein.  Hopping diffusion of nanoparticles in polymer matrices.  Macromolecules, 48(3):847–862, 2015.
  * [29] Yulong Chen, Rui Ma, Xin Qian, Ruoyu Zhang, Xifu Huang, Haohao Xu, Mi Zhou, and Jun Liu.  Nanoparticle mobility within permanently cross-linked polymer networks.  Macromolecules, 53(11):4172–4184, 2020.
  * [30] Justin B. Hooper and Kenneth S. Schweizer.  Contact aggregation, bridging, and steric stabilization in dense polymer-particle mixtures.  Macromolecules, 38(21):8858–8869, 2005.
  * [31] Valerio Sorichetti, Virginie Hugouvieux, and Walter Kob.  Dynamics of nanoparticles in polydisperse polymer networks: from free diffusion to hopping.  Macromolecules, 54(18):8575–8589, 2021.
  * [32] Hubert M. James.  Statistical properties of networks of flexible chains.  The Journal of Chemical Physics, 15(9):651–668, 1947.
  * [33] P. J. Flory.  Theory of elasticity of polymer networks. the effect of local constraints on junctions.  The Journal of Chemical Physics, 66(12):5720–5729, 1977.
  * [34] Michael Rubinstein and Sergei Panyukov.  Nonaffine deformation and elasticity of polymer networks.  Macromolecules, 30(25):8036–8044, 1997.
  * [35] PJ Hoogerbrugge and JMVA Koelman.  Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics.  EPL (Europhysics Letters), 19(3):155, 1992.
  * [36] W. Pan, I. V. Pivkin, and G. E. Karniadakis.  Single-particle hydrodynamics in dpd: A new formulation.  EPL (Europhysics Letters), 84(1):548–551, 2008.
  * [37] Wenxiao Pan, Bruce Caswell, and George Em Karniadakis.  Rheology, microstructure and migration in brownian colloidal suspensions.  Langmuir, 26(1):133–142, 2009.
  * [38] Li-Heng Cai, Sergey Panyukov, and Michael Rubinstein.  Mobility of nonsticky nanoparticles in polymer liquids.  Macromolecules, 44(19):7853–7863, 2011.
  * [39] Andrei A. Gusev.  Numerical estimates of the topological effects in the elasticity of gaussian polymer networks and their exact theoretical description.  Macromolecules, 52(9):3244–3251, 2019.
  * [40] P. J. Flory, M. Gordon, P. J. Flory, and N . G. McCrum.  Statistical thermodynamics of random networks.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 351(1666):351–380, 1976.