# Implications of Gamma Ray Burst GRB221009A for Extra Dimensions

Janus Capellan Aban, Chuan-Ren Chen, Yuan-Feng Hsieh, Chrisna Setyo Nugroho
Department of Physics, National Taiwan Normal University, Taipei 116, Taiwan

###### Abstract

Anomalous high energy photons, known as GRB221009A, with 18 TeV and 251 TeV
were observed by LHAASO and Carpet-2 recently. Such observation of high energy
gamma-ray bursts from distant source causes a mystery since high energy
photons suffer severe attenuation before reaching the earth. One possibility
is the existence of axion-like particles (ALP), and high energy photons at the
source can convert to these ALPs which travel intergalactically. In this
paper, we study the effects of extra dimensions on the conversion probability
between photon and ALPs. The conversion probability saturates and may reach
almost $100\%$ for high energy photons. We show that the size of extra
dimension affects the energy at which saturation occurs. The observations of
high-energy photons may support the possibility of smaller extra dimension.

##  I Introduction

Recently ultra high energy photons known as GRB221009A have been detected
S.Dichiara  through a range of instruments located on Earth, such as Fermi GBM
P.Veres , Carpet-2 D.Dzhappuev , and LHAASO Y.Huang ; LHAASO:2023kyg .
Notably, the LHAASO collaboration recorded the observation of high-energy
photons, exhibiting energies extending to around 18 TeV while Carpet-2
detected a 251 TeV photons shower. However, photons of cosmic origin cannot
traverse vast cosmological distances and reach Earth, as they undergo
substantial scattering due to the cosmic radiation background, transforming
into electron and positron pairs Gould:1966pza ; Fazio:1970pr ; A.Nikishov .
Consequently, the observation of energetic photons from distant sources is
prohibited by conventional cosmology, unless there exists non-standard new
physics that permits the occurrence of this phenomenon.

To account for the abundant detection of highly energetic photons from the
gamma-ray burst GRB221009A, a mechanism is required to prevent the attenuation
of photon flux as it traverses huge cosmological distances. One plausible
approach involves the incorporation of axion-like particle (ALP) $a$ that
interacts with photons via the relationship
$\dfrac{-g_{a\gamma}a}{4}F_{\mu\nu}\tilde{F}^{\mu\nu}=g_{a\gamma}a\vec{E}\cdot\vec{B}$.
The ALP represents a collective term that encompasses particles similar to the
Peccei-Quinn axion $a_{PQ}$, which is renowned for its solution of the strong
CP problem Peccei:1977hh . Given that the host galaxy of the GRB221009A source
possesses a galactic magnetic field of approximately $\textit{O}(1)\>\mu$G,
this magnetic field can induce the conversion of photons into ALP as caused by
ALP-photon interaction. Subsequently, this ALP can travel unhindered straight
to the Milky Way galaxy, and their conversion back into photons occurs due to
the influence of our galactic magnetic field, enabling the observation of
these highly energetic photons Maiani:1986md ; Raffelt:1987im ; Zhang:2022zbm
.

Since the ALPs are Standard Model (SM) singlets, they may propagate to the
large extra dimensions just like the gravitons and sterile neutrinos
Davoudiasl:2002fq ; Esmaili:2014esa ; Deffayet:2000pr . The theory of large
extra dimensions (LED) was pioneered in Arkani-Hamed:1998jmv ;
Antoniadis:1998ig  by proposing that the new scale of quantum gravity $M_{F}$
could be much lower compared to conventional Planck scale $M_{P}=10^{19}$ GeV
to solve the gauge hierarchy problem between weak scale and gravity. The scale
$M_{F}$ is related to the Planck scale $M_{P}$ by the fundamental equation
Antoniadis:1998ig  $\dfrac{M_{P}^{2}}{M_{F}^{\delta+2}}\simeq V_{\delta}\,$,
where $V_{\delta}$ is the volume of the $\delta$ compact extra dimensions. In
particular, for torus configuration, the volume of the extra dimensions is
denoted by $V_{\delta}=(2\pi)^{\delta}R_{1}R_{2}...R_{\delta}$, where $R_{i}$
for $i=1,2,...,\delta$ represents the radii of the extra dimensions. In this
case, the axion can propagate into the bulk or $(4+\delta)$ dimensions while
all the SM fields are confined in the brane or 4D flat-Minkowski spacetime.
The higher dimensional bulk axion field has Fourier decomposition that
produces a tower of axion Kaluza-Klein (KK) excitations or KK modes
Dienes:1999gw  as manifested in the brane after the compactification of the
$\delta$ extra dimensions, and the lowest KK mode is identified as the
ordinary ALP that we know in the brane. In this context, all the ALP KK modes
possess a coupling strength $g_{a\gamma}$ to photons and possess the ability
to undergo conversions to photons and vice versa, particularly in the presence
of external magnetic fields.

To see the impact of KK axions on the highly energetic photons observations on
Earth, we assume the same mechanism of the conventional photon-axion mixing.
The energetic photons from cosmic sources are mixed with KK axions in the
presence of the external magnetic field of the host galaxy. In return, some of
the photons are converted to KK axions that continue the journey unimpeded
across cosmological distance. The unconverted photons are severely attenuated
due to background light by annihilating into electron-positron pairs, meaning
such photons cannot travel far enough. The KK axions will arrive at the Milky
Way Galaxy and encounter galactic magnetic fields that will be reconverted
into energetic photons detectable on Earth.

The paper is structured as follows: In Section II, we briefly discuss the
model containing KK axion and photon interaction Deffayet:2000pr ;
Dienes:1999gw ; Chang:1999si  and derive the coupled Klein-Gordon and Maxwell
equations relevant for photon to KK axion conversion. The linearization of
this system of equations provides a compact expression leading to the mixing
matrix $M$. After the diagonalization of $M$ in Section III, we compute the
conversion probabilities of photon to KK axions. The discussions of Section IV
are about the constraints on axion parameters provided by Galanti:2022pbg .
Here, we study the effect of axion zero mode mass $m_{0}$, the extra dimension
radius $R$, photon-axion coupling $g_{a\gamma}$, as well the photon energy $E$
on photon to KK axions conversion. Finally, we present our conclusion in
Section V.

##  II Model

In the framework of LED, where hierarchical radii are present with a dominant
radius denoted as $R$ much larger than the others, the bulk system is
effectively equivalent to a 5D extra-dimensional model with extra-dimensional
radius $R$, as discussed in works such as Barbieri:2000mg . In this paper, we
focus on this type of LED model. We recall the generic bulk axion and photon
dynamics considered in Deffayet:2000pr ; Dienes:1999gw ; Chang:1999si , the
action is given as

$\displaystyle S_{4}=\int
d^{4}x\Bigg{[}\sum_{n}\Bigg{(}-\dfrac{1}{2}\\{\partial^{\mu}a^{(n)}\partial_{\mu}a^{(n)}+m_{n}^{2}a^{(n)^{2}}\\}-\dfrac{g_{a\gamma}a^{(n)}}{4}F_{\mu\nu}\tilde{F}^{\mu\nu}\Bigg{)}-\dfrac{1}{4}F_{\mu\nu}F^{\mu\nu}\Bigg{]}$
(1)

such that $m_{n}^{2}=\dfrac{n^{2}}{R^{2}}$ for $n\geq 1$ are the masses square
of the nonzero mode of Kaluza-Klein (KK) axions $a^{(n)}$, while
$\tilde{F}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}$ is
the dual of the electromagnetic field tensor
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, and
$\epsilon_{\mu\nu\rho\sigma}$ is the antisymmetric tensor satisfying
$\epsilon_{0123}=+1$. In addition, $m_{0}$ denotes the mass of zeroth mode
$a^{(0)}$ of KK axions. Here, each of the KK axion states $a^{(n)}$ couples to
photon with coupling strength $g_{a\gamma}$.

We assume that the photon field $A^{\mu}=(0,\vec{A})$ produced from the
source, travels across its host galaxy and encounters external magnetic field
$\vec{H}_{0}$. Throughout the discussion, we assume that a monochromatic beam
of light is traveling in $u$ direction with the vector component of the photon
field $\vec{A}$ as in Deffayet:2000pr

$\displaystyle\vec{A}(u,t)=i\,(\vec{A}_{\times}(u),\vec{A}_{\plus}(u),0)e^{-i\omega
t}\,,$ (2)

where the coordinate $u$ is along the direction of propagation. The
corresponding electric and magnetic fields are expressed as

$\displaystyle\vec{E}\equiv-\partial_{t}\vec{A}(u,t)=(\omega\vec{A}_{\times}(u),\omega\vec{A}_{\plus}(u),0)e^{-i\omega
t}\,,$ (3)
$\displaystyle\vec{B}\equiv\nabla\times\vec{A}=(-i\partial_{u}A_{\plus}(u),i\partial_{u}A_{\times}(u),0)e^{-i\omega
t}\,.$ (4)

The equations of motions for the pseudoscalar field axion and photon field
that can be determined from the action in Eq. (1) are given by the coupled
Klein-Gordon and Maxwell equations Deffayet:2000pr

$\displaystyle(\Box-
m_{\vec{p}}^{2})a^{(n)}=\dfrac{g_{a\gamma}}{4}F_{\mu\nu}\tilde{F}^{\mu\nu}\,,$
(5)
$\displaystyle\partial_{\alpha}F^{\alpha\beta}=-g_{a\gamma}\,\partial_{\alpha}\Big{[}\sum_{n}a^{(n)}\tilde{F}^{\alpha\beta}\Big{]}\,.$
(6)

Next, we define a set of orthonormal basis

$\displaystyle\vec{e}_{||}=\frac{\vec{u}}{u},\>\>\>\vec{e}_{\times}=\frac{\vec{H}_{0\perp}}{H_{0\perp}},\>\>\>\vec{e}_{\plus}\,,$
(7)

such that $\vec{H}_{0\perp}$ is the orthogonal component of the external
magnetic field $\vec{H}_{0}$ with respect to the direction of the photon
propagation $\vec{u}$. Furthermore, by writing the axion field as
$a^{(\vec{p})}(u,t)=a^{(\vec{p})}(u)e^{-i\omega t}$, Eqs. (5) - (6) can be
simplified into Deffayet:2000pr

$\displaystyle(\Box-m_{n}^{2})a^{(n)}=-g_{a\gamma}\omega
H_{0\perp}A_{\times}\,,$ (8) $\displaystyle\Box A_{\lambda}=-g_{a\gamma}\omega
H_{0\perp}\delta_{\lambda\times}\sum_{n}a^{(n)}\,,$ (9)

where $\lambda$ stands for photon polarization vector $\lambda=+,\times$. Note
that when $\lambda=+$ the axion states decouple from the photon. This is
because only the parallel component of the photon to $\vec{H}_{0\perp}$
contributes to the of axion-photon dynamics. Therefore, in our case, we
neglect the component $A_{\plus}$ of the photon.

Before simplifying Eqs. (8) - (9) directly, let’s consider the dispersion
relation of the photon and external magnetic field. It is assumed that the
external magnetic field varies very much larger than the photon wavelength in
space. Thus, we may write the D’lambert operator $\Box$ as

$\displaystyle\Box=(\omega^{2}+\partial_{u}^{2})=(\omega+i\partial_{u})(\omega-i\partial_{u})\simeq(\omega+k)(\omega-i\partial_{u})=2\omega(\omega-i\partial_{u})\,,$
(10)

having the refractive index $n$ satisfies the form of general dispersion
relation $\omega\simeq nk$ with $|n-1|<<1$ in our case. Thus, Eqs. (8) - (9)
can be further simplified as Deffayet:2000pr

$\displaystyle(\omega-i\partial_{u}-\dfrac{m_{n}^{2}}{2\omega})a^{(n)}=-\dfrac{g_{a\gamma}}{2}H_{0\perp}A_{\times}\,,$
(11)
$\displaystyle(\omega-i\partial_{u})A_{\lambda}=-\dfrac{g_{a\gamma}}{2}H_{0\perp}\delta_{\lambda\times}\sum_{n}a^{(n)}\,.$
(12)

These coupled equations may be written in a compact form

$\displaystyle(\omega-i\partial_{u}+M)\begin{pmatrix}A_{\times}\\\
a^{({0})}\\\ \vdots\\\ a^{({n})}\\\ \vdots\end{pmatrix}=0\,,$ (13)

with the corresponding mixing matrix $M$ given by

$\displaystyle
M=\begin{pmatrix}\Delta_{\times}&\Delta_{M}&\Delta_{M}&\ldots&\\\
\Delta_{M}&\Delta_{a}^{({0})}&0&0&\ldots&\\\
\Delta_{M}&0&\ddots&\ddots&\ddots&\\\
\vdots&\vdots&\vdots&\Delta_{a}^{({n})}&\ddots&\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\end{pmatrix}\,.$ (14)

Here, $\Delta_{M}=\dfrac{g_{a\gamma}H_{0\perp}}{2}$,
$\Delta_{a}^{(n)}=-\dfrac{m_{n}^{2}}{2\omega}$, and
$\Delta_{\times}=-\dfrac{\omega_{pl}^{2}}{2\omega}$ with
$\omega_{pl}=\dfrac{4\pi\alpha n_{e}}{m_{e}}$. The latter expression is the
photon mass with $n_{e}$, $m_{e}$, and $\alpha$ denoting electron density,
electron mass, and the fine structure constant, respectively.

##  III Photon-Axion Conversion Probability

The high-energy photons detected on Earth can be attributed to two potential
astronomical sources: magnetars, which are newly born, rapidly rotating stars
with intense magnetic fields, and supernovae, resulting from the collapse of
exceptionally massive stars Granot:2015xba ; MacFadyen:1998vz . In this
scenario, we consider that the extremely energetic photons may originate from
the mentioned astronomical sources and travel as concentrated gamma-ray jets
throughout the universe. Concerning the mechanism of the conversion
$\gamma\rightarrow a\rightarrow\gamma$ Wang:2023okw , these photons pass
through the host galaxy, they interact with magnetic fields leading to the
conversion of a portion of the photon fluxes into KK axions. The resulting KK
axions remain undisturbed as they journey into the intergalactic region, while
the photons are scattered by the background photons present in the surrounding
space. The KK axions are then converted back into photons as it arrives in the
Milky Way galaxy.

To compute the conversion probability, we need to solve Eq. (13) with the
procedure provided in Raffelt:1987im  by considering the solution of the form

$\displaystyle\vec{V}(u)=e^{-iMu}e^{-i\omega u}\vec{V}(0)\,,$ (15)

where $\vec{V}(u)$ is given by a column vector

$\displaystyle\vec{V}(u)=\begin{pmatrix}A(u)\\\ a^{0}(u)\\\ a^{1}(u)\\\
\vdots\\\ a^{N}(u)\\\ \end{pmatrix}\,,$ (16)

where we use the notation $A(u)$ to be the photon state for $A_{\times}(u)$
and $a^{n}(u)$ are the axion Kaluza-Klein modes for $n=0,1,2,...,N$. Suppose
we consider the eigenvector $\vec{v}_{n}=(E_{n\gamma},E_{n0},...,E_{nN})^{T}$
corresponds to the eigenvalue $\lambda_{n}$ of $M$, where $n=\gamma,0,1,...,N$
which satisfies the eigenvalue equation $(M-\lambda_{n}\cdot I)\vec{v}_{n}=0$

$\displaystyle\begin{pmatrix}\Delta_{\times}-\lambda_{n}&\Delta_{M}&\Delta_{M}&\ldots&\Delta_{M}&\\\
\Delta_{M}&\Delta_{a}^{({0})}-\lambda_{n}&0&\ldots&0&\\\
\Delta_{M}&0&\Delta_{a}^{({1})}-\lambda_{n}&\ldots&0&\\\
\vdots&\vdots&\vdots&\ddots&0&\\\
\Delta_{M}&0&0&0&\Delta_{a}^{({N})}-\lambda_{n}&\end{pmatrix}\begin{pmatrix}E_{n\gamma}\\\
E_{n0}\\\ E_{n1}\\\ \vdots\\\ E_{nN}\\\ \end{pmatrix}=0\,.$ (17)

This eigenvalue equation yields the $N+2$ relations

$\displaystyle
E_{n\gamma}(\Delta_{\times}-\lambda_{n})+\Delta_{M}\sum_{k=0}^{N}E_{nk}$
$\displaystyle=0\,,$ (18)
$\displaystyle\Delta_{M}E_{n\gamma}+(\Delta_{a}^{({0})}-\lambda_{n})E_{n0}$
$\displaystyle=0\,,$ (19)
$\displaystyle\Delta_{M}E_{n\gamma}+(\Delta_{a}^{({1})}-\lambda_{n})E_{n1}$
$\displaystyle=0\,,$ (20) $\displaystyle\vdots$ (21)
$\displaystyle\Delta_{M}E_{n\gamma}+(\Delta_{a}^{({N})}-\lambda_{n})E_{nN}$
$\displaystyle=0\,.$ (22)

More compactly, aside from the first entry in Eq. (18) we have the relation

$\displaystyle
E_{nj}=\dfrac{\Delta_{M}E_{n\gamma}}{(\lambda_{n}-\Delta_{a}^{({j})})}\,,$
(23)

with $j=0,1,...,N$. Thus, the unitary matrix that diagonalizes the matrix $M$
is given by

$\displaystyle
U=\begin{pmatrix}E_{\gamma\gamma}&E_{0\gamma}&E_{1\gamma}&\ldots&E_{N\gamma}&\\\
E_{\gamma 0}&E_{00}&E_{10}&\ldots&E_{N0}&\\\ E_{\gamma
1}&E_{01}&E_{11}&\ldots&E_{N1}&\\\ \vdots&\vdots&\vdots&\ddots&0&\\\ E_{\gamma
N}&E_{0N}&E_{1N}&\ldots&E_{NN}&\\\ \end{pmatrix}\,.$ (24)

From the unitarity condition $U^{\dagger}U=I$ we have

$\displaystyle|E_{n\gamma}|^{2}+|E_{n0}|^{2}+\ldots+|E_{nN}|^{2}=1\,.$ (25)

We substitute Eq. (23) to express the mixing matrix component $E_{n\gamma}$ in
terms of available parameters in $M$ and its respective eigenvalues
$\lambda_{n}$, the unitarity condition yields

$\displaystyle
E_{n\gamma}^{2}+\Delta_{M}^{2}\sum_{j=0}^{N}\frac{E_{n\gamma}^{2}}{(\lambda_{n}-\Delta_{a}^{(j)})^{2}}=1\,,$
$\displaystyle\implies$ $\displaystyle
E_{n\gamma}=\dfrac{1}{\Big{(}1+\sum_{j=0}^{N}{\frac{\Delta_{M}^{2}}{(\lambda_{n}-\Delta_{a}^{(j)})^{2}}}\Big{)}^{1/2}}\,.$
(26)

Hence, the eigenvector of $M$ can be written as

$\displaystyle\vec{v}_{n}=\begin{pmatrix}E_{n\gamma}\\\ E_{n0}\\\ E_{n1}\\\
\vdots\\\ E_{nN}\\\ \end{pmatrix}=\begin{pmatrix}E_{n\gamma}\\\
\frac{\Delta_{M}E_{n\gamma}}{(\lambda_{n}-\Delta_{a}^{(0)})}\\\
\frac{\Delta_{M}E_{n\gamma}}{(\lambda_{n}-\Delta_{a}^{(1)})}\\\ \vdots\\\
\frac{\Delta_{M}E_{n\gamma}}{(\lambda_{n}-\Delta_{a}^{(N)})}\\\
\end{pmatrix}=E_{n\gamma}\begin{pmatrix}1\\\
\frac{\Delta_{M}}{(\lambda_{n}-\Delta_{a}^{(0)})}\\\
\frac{\Delta_{M}}{(\lambda_{n}-\Delta_{a}^{(1)})}\\\ \vdots\\\
\frac{\Delta_{M}}{(\lambda_{n}-\Delta_{a}^{(N)})}\\\
\end{pmatrix}=\dfrac{1}{\Big{(}1+\sum_{j=0}^{N}{\frac{\Delta_{M}^{2}}{(\lambda_{n}-\Delta_{a}^{(j)})^{2}}}\Big{)}^{1/2}}\begin{pmatrix}1\\\
\frac{\Delta_{M}}{(\lambda_{n}-\Delta_{a}^{(0)})}\\\
\frac{\Delta_{M}}{(\lambda_{n}-\Delta_{a}^{(1)})}\\\ \vdots\\\
\frac{\Delta_{M}}{(\lambda_{n}-\Delta_{a}^{(N)})}\\\ \end{pmatrix}\,.$ (27)

If we assume initially that we have photon-axion state
$\vec{V}(0)=(A(0),a^{0}(0),a^{1}(0),...,a^{N}(0))^{T}$, then it can be
expressed as a linear combination of our basis $\vec{v}_{n}$ as a consequence
of diagonalization, i.e.
$\vec{V}(0)=U(\vec{v}_{\gamma},\vec{v}_{0},...,\vec{v}_{N})^{T}$. In this
case, the initial state is given by

$\displaystyle
A(0)=\sum_{k=\gamma,0}^{N}E_{k\gamma}\vec{v}_{k}\>\>\>\>\textrm{and}\>\>\>\>\>a^{n}(0)=\sum_{k=\gamma,0}^{N}E_{kn}\vec{v}_{k}\>\>\>\>\textrm{for}\>\>\>\>n=0,1,...,N$
(28)

From Eq. (15), the photon state initially given by $\vec{A}(0)$ would evolve
into $\vec{A}(u)$ as it traverse a distance $u$

$\displaystyle\vec{A}(u)$ $\displaystyle=e^{-iMu}e^{-i\omega u}\vec{A}(0)\,,$
$\displaystyle=e^{-i\omega
u}\sum_{k=\gamma,0}^{N}E_{k\gamma}e^{-i\lambda_{k}u}\vec{v}_{k}\,,$
$\displaystyle=e^{-i\omega
u}\sum_{k=\gamma,0}^{N}e^{-i\lambda_{k}u}\dfrac{1}{\Big{(}1+\sum_{j=0}^{N}{\frac{\Delta_{M}^{2}}{(\lambda_{k}-\Delta_{a}^{(j)})^{2}}}\Big{)}}(1,\frac{\Delta_{M}}{(\lambda_{k}-\Delta_{a}^{(0)})},\frac{\Delta_{M}}{(\lambda_{k}-\Delta_{a}^{(1)})},\cdots,\frac{\Delta_{M}}{(\lambda_{k}-\Delta_{a}^{(N)})})^{T}\,.$
(29)

Thus, the conversion probability of photon to KK axions is given by
Deffayet:2000pr

$\displaystyle P(A(0)\rightarrow a(u))$ $\displaystyle=1-P(A(0)\rightarrow
A(u))\,,$ $\displaystyle=1-|\langle{A(u)}|{A(0)}\rangle|^{2}\,,$
$\displaystyle=1-|\langle{e^{-i\omega
u}\sum_{k=\gamma,0}^{N}E_{k\gamma}e^{-i\lambda_{k}u}\vec{v}_{k}}|{A(0)}\rangle|^{2}\,,$
$\displaystyle=1-\Big{|}\sum_{k=\gamma,0}^{N}e^{+i\lambda_{k}u}f_{k}^{2}\Big{|}^{2}\,,$
(30)

with the normalization factor $f^{2}_{k}$ given by

$\displaystyle
f_{k}^{2}=\dfrac{1}{\Big{(}1+\sum_{j=0}^{N}{\frac{\Delta_{M}^{2}}{(\lambda_{k}-\Delta_{a}^{(j)})^{2}}}\Big{)}}\,,\text{with}\,\sum^{N}_{k}f^{2}_{k}=1\,,\textrm{for}\>\>\>k=\gamma,0,1,...,N\,.$
(31)

One may write the conversion probability in a more compact form
Deffayet:2000pr

$\displaystyle P(A(0)\rightarrow a(u))\equiv P(\gamma\rightarrow
a)=2\,\sum^{N}_{i,j}f^{2}_{\lambda_{i}}f^{2}_{\lambda_{j}}\,\text{sin}^{2}\frac{(\lambda_{i}-\lambda_{j})u}{2}\,\,\,,\textrm{for}\>\>\>i,j=\gamma,0,1,...,N\,.$
(32)

We solve numerically the eigenvalues of matrix $M$ and further compute the
photon to axion conversion probability given in Eq. (III). To reduce the
computer power, we include up to 10 KK axions in our numerical results, since
our conclusions remain even more KK modes are involved.

##  IV Results and Discussion

In the energy range between 100 GeV and 300 TeV relevant for our study as well
as the extra dimension size less than 30 $\mu$m Lee:2020zjt , the components
of the matrix $M$ exhibit some hierarchies $\Delta_{\times}\approx 0$,
$\Delta^{(n)}_{a}>>\Delta^{(0)}_{a}\,,\Delta_{M}$. Therefore, the eigenvalues
of the matrix $M$ may be well approximated as

$\displaystyle\lambda_{\gamma}$
$\displaystyle\approx\frac{\Delta^{(0)}_{a}}{2}-\Delta_{M}-\frac{(\Delta^{(0)}_{a})^{2}}{8\Delta_{M}}-\frac{\Delta^{2}_{M}}{\Delta^{(n)}_{a}}\,,$
(33) $\displaystyle\lambda_{0}$
$\displaystyle\approx\frac{\Delta^{(0)}_{a}}{2}+\Delta_{M}+\frac{(\Delta^{(0)}_{a})^{2}}{8\Delta_{M}}\,,$
(34) $\displaystyle\lambda_{n}$
$\displaystyle\approx\Delta^{(n)}_{a}+\frac{\Delta^{2}_{M}}{\Delta^{(n)}_{a}},~{}\rm{for}~{}n\neq
0\,.$ (35)

Furthermore, for extra dimension size of $10\,\mu\textrm{m}\leq R\leq
30\mu\textrm{m}$, the conversion probability can be written analytically as

$\displaystyle P(\gamma\rightarrow a)$ $\displaystyle\approx
4f^{2}_{\gamma}\,f^{2}_{0}\,\text{sin}^{2}\left[\left(\Delta_{M}+\frac{(\Delta^{(0)}_{a})^{2}}{8\Delta_{M}}+\frac{\Delta^{2}_{M}}{2\Delta^{(n)}_{a}}\right)u\right]\,,$
(36)

for a fixed $n$, where the expression of the normalization factor
$f^{2}_{\gamma}$ and $f^{2}_{0}$ are given by

$\displaystyle
f^{2}_{\gamma}\approx\frac{1}{1+\frac{\Delta^{2}_{M}}{\left(\Delta_{M}+\frac{\Delta^{(0)}_{a}}{2}\right)^{2}}+\frac{\pi^{4\,\Delta^{2}_{M}}}{90\,(\Delta^{(1)}_{a})^{2}}}\,\,\,,\,f^{2}_{0}\approx\frac{1}{1+\frac{\Delta^{2}_{M}}{\left(\Delta_{M}-\frac{\Delta^{(0)}_{a}}{2}\right)^{2}}+\frac{\pi^{4\,\Delta^{2}_{M}}}{90\,(\Delta^{(1)}_{a})^{2}}}\,.$
(37)

We may neglect $\Delta^{(1)}_{a}$ term in Eq.(37) since it has much larger
value than other parameters. In addition, as energy gets higher,
$\Delta^{(0)}_{a}$ becomes irrelevant and the conversion probability does not
depend on $E$. In this limit, both $f^{2}_{\gamma}$ and $f^{2}_{0}$ approach
$1/2$ and the contribution of other $f^{2}_{k}$ are negligible based on the
normalization condition in Eq. (31).

In our computation, we employ the recent constraint on the size of the extra
dimension which is less than 30 $\mu$m Lee:2020zjt . In addition, we follow
the constraints on axion to photon coupling $g_{a\gamma}$ and the zeroth mode
of axion mass $m_{0}$ from Galanti:2022pbg ; Wang:2023okw :

$\displaystyle 5.0\times 10^{-12}\>\textrm{GeV}^{-1}\leq g_{a\gamma}\leq
2.1\times
10^{-11}\textrm{GeV}^{-1}\>\>\>\textrm{for}\>\>\>10^{-10}\>\textrm{eV}\leq
m_{0}\leq 2.0\times 10^{-7}\textrm{eV}.$ (38)

To get a better understanding of our results, we put some typical values of
the matrix elements $M$

$\displaystyle\Delta_{\times}$ $\displaystyle\simeq\big{(}-1.1\times
10^{-10}\big{)}\Big{(}\dfrac{n_{e}}{10^{-3}\textrm{cm}^{3}}\Big{)}\Big{(}\dfrac{E}{\textrm{TeV}}\Big{)}^{-1}\textrm{kpc}^{-1}\,,$
(39) $\displaystyle\Delta_{M}$ $\displaystyle\simeq\big{(}3.1\times
10^{-2}\big{)}\Big{(}\dfrac{g_{a\gamma}}{2\times
10^{-11}\textrm{GeV}^{-1}}\Big{)}\Big{(}\dfrac{H_{0\perp}}{\mu\textrm{G}}\Big{)}\textrm{kpc}^{-1}\,,$
(40) $\displaystyle\Delta_{0}$ $\displaystyle\simeq(-7.8\times
10^{-3})\Big{(}\dfrac{m_{0}}{10^{-8}\textrm{eV}}\Big{)}^{2}\Big{(}\dfrac{E}{\textrm{TeV}}\Big{)}^{-1}\textrm{kpc}^{-1}\,,$
(41) $\displaystyle\Delta_{n}$ $\displaystyle\simeq(-3.04n^{2}\times
10^{12})\Big{(}\dfrac{\mu\textrm{m}}{R}\Big{)}^{2}\Big{(}\dfrac{E}{\textrm{TeV}}\Big{)}^{-1}\textrm{kpc}^{-1}\,.$
(42)

As an illustration, setting a particular set of parameters corresponds to

$\displaystyle n_{e}=1.0\times
10^{-3}\textrm{cm}^{-3}\>,u=30\textrm{kpc},\>H_{0\perp}=5\mu\textrm{G},\>g_{a\gamma}=2\times
10^{-11}\textrm{GeV}^{-1},\>R=10\mu\textrm{m}\,,$ (43)

the corresponding matrix elements for 2 KK axions are

$\displaystyle M=\textrm{kpc}^{-1}\begin{pmatrix}-1.1\times
10^{-10}\Big{(}\dfrac{E}{\textrm{TeV}}\Big{)}^{-1}&1.55\times
10^{-1}&1.55\times 10^{-1}&1.55\times 10^{-1}\\\ 1.55\times 10^{-1}&-7.8\times
10^{-3}\Big{(}\dfrac{E}{\textrm{TeV}}\Big{)}^{-1}&0&0\\\ 1.55\times
10^{-1}&0&-3.04\times 10^{10}\Big{(}\dfrac{E}{\textrm{TeV}}\Big{)}^{-1}&0\\\
1.55\times 10^{-1}&0&0&-1.22\times
10^{11}\Big{(}\dfrac{E}{\textrm{TeV}}\Big{)}^{-1}\\\ \end{pmatrix}\,.$ (44)

If we set the photon energy $E=10$ TeV, the numerical values of this matrix
elements become

$\displaystyle M=\textrm{kpc}^{-1}\begin{pmatrix}-1.1\times
10^{-11}&1.55\times 10^{-1}&1.55\times 10^{-1}&1.55\times 10^{-1}\\\
1.55\times 10^{-1}&-7.8\times 10^{-4}&0&0\\\ 1.55\times 10^{-1}&0&-3.04\times
10^{9}&0\\\ 1.55\times 10^{-1}&0&0&-1.22\times 10^{10}\\\ \end{pmatrix}\,.$
(45)

Clearly, the mixing matrix $M$ is dominated by the masses of nonzero modes of
KK axion which can be found in the last two diagonal entries. The eigenvalues
of the $4\times 4$ matrix $M$ are given by

$\displaystyle\lambda_{\gamma}=-0.1546\>\textrm{kpc}^{-1},\>\lambda_{0}=0.1554\>\textrm{kpc}^{-1},\>\lambda_{1}=-3.039\times
10^{9}\>\textrm{kpc}^{-1}\>\textrm{and}\>\lambda_{2}=-1.216\times
10^{10}\>\textrm{kpc}^{-1}$ (46)

with $f_{\gamma}^{2}=0.501$, $f_{0}^{2}=0.499$, $f_{1}^{2}\approx 0$ and
$f_{2}^{2}\approx 0$ as calculated using Eq. (31). The corresponding
conversion probability is $P(\gamma\rightarrow a)=0.996$.

In the upper panel of Fig.1, we show the conversion probability as a function
of galactic magnetic field for photon energy $E=18$ TeV reported by LHAASO
collaboration. The value of $P(\gamma\rightarrow a)$ oscillates as the
magnetic field changes. This behaviour is expected as one can see from Eq.(36)
which shows that $P(\gamma\rightarrow a)$ depends on the square of sinusoidal
function of $\Delta_{M}\propto g_{a\gamma}H_{0\perp}$. On the other hand, from
the lower panel of Fig.1, $P(\gamma\rightarrow a)$ is independent of $E$ in
TeV regime since $\Delta^{(0)}_{a}$ is negligible for $m_{0}=10^{-8}$ GeV and
$R=25\,\mu$m. In this case, both $f^{2}_{\gamma}$ and $f^{2}_{0}$ approach
$1/2$ as reflected by Eq.(37).

Figure 1:  The upper panel shows the conversion probability of 18 TeV photons
into KK axions as a function of transverse magnetic field $H_{0\perp}$. The
lower panel shows the plot of photon to KK axions conversion probability with
varying external magnetic field $H_{0\perp}=1\mu$G, $3\mu$G, $5\mu$G, and
$10\mu$G as a photon travels a distance of $u=30$ kpc from the source. The
axion-photon coupling strength used is $g_{a\gamma}=2\times 10^{-11}$GeV-1
while the mass of axion zero mode is set to $m_{0}=10^{-8}$ eV. The size of
the dominant extra dimension in this case is $R=25\mu$m.

Figure 2:  The two figures show the conversion probability on 18 TeV photons
in the presence of external magnetic field $H_{0\perp}=5\mu$G with radius of
extra dimension $R=25\mu$m as the photon travels a distance of $u=30$ kpc from
the source. The upper panel shows the plot of conversion probability as a
function of KK axion-photon coupling $g_{a\gamma}$ with different axion zero
mode masses $m_{0}=10^{-10}$ eV, $10^{-9}$ eV, $10^{-8}$ eV, $10^{-7}$ eV. The
lower panel shows the plot of photon to KK axions conversion probability as a
function of axion zero mode mass $m_{0}$ with different couplings
$g_{a\gamma}=5\times 10^{-12}$GeV${}^{-1},10^{-11}$GeV${}^{-1},5\times
10^{-11}$GeV${}^{-1},10^{-10}$GeV-1.

The oscillatory behaviour of $P(\gamma\rightarrow a)$ with respect to
$g_{a\gamma}$ can also be seen in the upper panel of Fig.2. Here, for
$m_{0}>10^{-8}$ eV at $g_{a\gamma}=7\times 10^{-12}\text{GeV}^{-1}$, the value
of $\Delta^{(0)}_{a}$ becomes larger. Consequently, the value of $f^{2}_{0}$
is decreased ($\Delta^{(0)}_{a}$ has negative value), resulting in lower value
of conversion probability, see Eq.(36) and Eq.(37). As $g_{a\gamma}$
increases, the value of $\Delta_{M}$ dominates $\Delta^{(0)}_{a}$, and the
oscillation amplitude returns to the situation that
$f^{2}_{\gamma}=f^{2}_{0}\approx 1/2$. On the other hand, the lower panel of
the same figure shows that $P(\gamma\rightarrow a)$ drops significantly for
$m_{0}>10^{-7}$ eV. This is because as $m_{0}$ increases, the value of
$\Delta^{(0)}_{a}$ would be very close to $\Delta_{M}$ leading to the
cancellation in the denominator of $f^{2}_{\gamma}$ in Eq.(37)

$\displaystyle
f^{2}_{\gamma}\approx\frac{1}{1+\Delta^{2}_{M}/\epsilon^{2}}\approx\frac{\epsilon^{2}}{\Delta^{2}_{M}}\,,$
(47)

where $\epsilon=\Delta_{M}-|\Delta^{(0)}_{a}|/2$ is close to zero. As a
result, the conversion probability is suppressed while retaining its
oscillatory behaviour controlled by
$\text{sin}^{2}((\Delta^{0}_{a})^{2}\,u/8\Delta_{M})$.

Figure 3: The upper plots show the conversion probability of photon into KK
axions with different sizes of extra dimensions. The black horizontal line
indicates the 95 percent conversion probability. The lower plot shows the
conversion probability of 18 TeV photons into KK axions as a function of the
size of extra dimension $R$. The magnetic field we consider is
$H_{0\perp}=5\mu$G with axion zero mode mass $m_{0}=10^{-8}$eV.

Next, the effects of $R$ on the conversion probability is displayed in Fig.3.
At both low as well as high energy regime, $P(\gamma\rightarrow a)$ oscillates
rapidly when $R$ is less than 10 $\mu$m as one can see from the upper panel.
The same behaviour is also observed in the lower panel of Fig.3 where
$P(\gamma\rightarrow a)$ saturates at $R\geq$ 10 $\mu$m. Here, one can not use
Eq.(36) to describe conversion probability as it is suppressed by
$\Delta^{2}_{M}/2\Delta^{(a)}_{n}$. The rapid oscillation is due to the
difference between two eigenvalues
$\lambda_{m}-\lambda_{n}\approx\Delta^{(a)}_{m}-\Delta^{(a)}_{n}$ appearing in
the sinusoidal argument of Eq.(32) with $m\neq n$ and $m,n\neq\gamma,0$. Thus,
the corresponding conversion probability is dominated by KK higher modes

$\displaystyle P(\gamma\rightarrow
a)=2\,\sum^{N}_{m,n}f^{2}_{m}f^{2}_{n}\,\text{sin}^{2}\frac{(\Delta^{(a)}_{m}-\Delta^{(a)}_{n})u}{2}\,=2\,\sum^{N}_{m,n}f^{2}_{m}f^{2}_{n}\,\text{sin}^{2}\frac{(m^{2}-n^{2})u}{4E\,R^{2}}\,,$
(48)

where $m,n=1,...,N$ and $R\leq 10\mu$m. We can see that the conversion
probability is given by the superposition of highly oscillated function
leading to spiky oscillation pattern. Such behaviour is very different from
the smooth oscillation given by Eq.(36) which is controlled by the photon and
zero mode axion. As an additional remark, the amplitude of conversion
probability $f^{2}_{n}$ in this case reads

$\displaystyle
f^{2}_{n}=\frac{1}{1+\Delta^{2}_{M}\,\left(\frac{1}{(\lambda_{n}-\Delta^{(1)}_{a})^{2}}+...+\frac{1}{(\lambda_{n}-\Delta^{(n)}_{a})^{2}}+...+\frac{1}{(\lambda_{n}-\Delta^{(N)}_{a})^{2}}\right)}\,.$
(49)

Note that the term $1/(\lambda_{n}-\Delta^{(j)}_{a})$ in the denominator of
$f^{2}_{n}$ is maximized at resonance for $j=n$. The value of $f^{2}_{n}$ and
the associated conversion probability at resonance are given by

$\displaystyle f^{2}_{n,res}$
$\displaystyle\approx\Delta^{2}_{M}/(\Delta^{(n)}_{a})^{2}=\frac{4E^{2}\,R^{4}\,\Delta^{2}_{M}}{n^{4}}\,,$
$\displaystyle P_{res}(\gamma\rightarrow a)$
$\displaystyle=8\sum^{N}_{m}f^{2}_{m}\,E^{2}\frac{R^{4}}{n^{4}}\,\Delta^{2}_{M}\,\text{sin}^{2}\frac{(m^{2}-n^{2})u}{4E\,R^{2}}\,,\text{for
fix}\,\,n\,.$ (50)

We see that all conversion probabilities converge to a constant value at high
energy, see Eq.(36), the second line of Eq.(IV), and lower panel of Fig. 1 and
Fig. 3. Although for the latter case, this saturation depends on the size of
extra dimension $R$. For similar $R$, the saturation occurs at higher energy.
To optimize the high flux of axions from the source, we define the critical
energy $E_{c}$ as the energy where $P(\gamma\rightarrow a)\geq 95\%$ for a
given radius of extra dimension $R$. The $E_{c}-R$ plane is depicted in Fig.4
using $g_{a\gamma}=2\times 10^{-11}\text{GeV}^{-1}$ and $g_{a\gamma}=7\times
10^{-12}\text{GeV}^{-1}$. These two couplings correspond to the first two
peaks in the upper panel of Fig.2. Notice that when $R$ is as large as
$10~{}\mu\rm{m}$, the $E_{c}$ becomes insensitive to $R$, since the conversion
probability $P(\gamma\to a)$ is dominated by the conversion between $\gamma$
and zeroth mode of axion. The red dashed line denotes the upper limit of $R$
from Lee:2020zjt . The energy photons observed by LHAASO ($18$ TeV) and
Carpet-2 ($251$ TeV) are shown in vertical green and purple lines,
respectively. For $g_{a\gamma}=2\times 10^{-11}$ GeV-1 ($7\times
10^{-12}\text{GeV}^{-1}$), the LHAASO data implies that, very likely, the
extra dimension radius $R$ could be as small as $0.43~{}\mu\rm{m}$
($0.75~{}\mu\rm{m}$), while Carpet-2 supports $R\gtrsim 0.11~{}\mu\rm{m}$
($0.19~{}\mu\rm{m}$).

Figure 4: The plot shows the relationship between critical energy $E_{c}$ and
extra dimension radius $R$ at 95 percent conversion probability of photons
into KK axions with magnetic field $H_{0\perp}=5\mu$G and mass of axion zero
mode $m_{0}=10^{-8}$eV. The vertical green dashed line is the maximum energy
(18 TeV) detected by LHAASO Y.Huang ; LHAASO:2023kyg  while the vertical
purple dashed line corresponds to 251 TeV as observed by Carpet-2 D.Dzhappuev
caused by GRB221009A. The horizontal red dashed line is the new upper bound
$R=30\mu$m of the size of extra dimension from the test of Newtonian gravity
at short distance Lee:2020zjt .

##  V Conclusion

In this study, we have demonstrated that the detection of high-energy photons
can be explained by KK axions and photon dynamics in the presence of external
magnetic fields. Highly energetic photons, emitted from astronomical sources
such as supernovae and magnetars, traverse the host galaxy while encountering
its magnetic field. Consequently, some of these photons undergo conversion
into axions, seamlessly continuing their journey across cosmological distances
unimpeded. The KK axions then arrive in the Milky Way galaxy and are converted
back into photons the moment they encounter the magnetic field in Milky Way,
subsequently observed on Earth. Throughout our analysis, we used 10 nonzero
modes of the Kaluza-Klein (KK) axions. We found that an external magnetic
field of $5\mu$G is highly efficient in converting highly energetic photons
into KK axions traversing a distance of $u=30$ kpc. For instance, if the
photon’s energy is around 18 TeV (the highest energy detected by the LHAASO
collaborations) the conversion probability to KK axions reaches almost 1 if
$g_{a\gamma}=7\times 10^{-12}$ GeV-1 and $g_{a\gamma}=2\times 10^{-11}$ GeV-1
in the axion zero mode mass regime $m_{0}\leq 10^{-7}$ eV, which explains well
why such high energy photon events can be observed. Additionally, we found
that the size of the extra dimension cannot be arbitrarily small, as it would
diminish the conversion probability of photons to axions, hindering the
detection of high-energy photons on Earth. In fact, for $H_{0\perp}=5\,\mu$G
and $m_{0}=10^{-8}$ eV, achieving efficient conversion of 18 TeV photons into
KK axions at a 95 percent probability or higher implies an extra-dimensional
size window of $R\in[0.43\mu m,30\mu m]$ and $R\in[0.75\mu m,30\mu m]$ for
coupling strengths $g_{a\gamma}=2\times 10^{-11}$ GeV-1 and
$g_{a\gamma}=7\times 10^{-12}$ GeV-1, respectively. Similarly, for 251 TeV
photons (the energy of photons detected by Carpet-2), the corresponding
windows should be likely $R\in[0.11\mu m,30\mu m]$ and $R\in[0.19\mu m,30\mu
m]$, respectively.

## Acknowledgment

We would like to acknowledge the support of National Center for Theoretical
Sciences (NCTS). This work was supported in part by the National Science and
Technology Council (NSTC) of Taiwan under Grant No. NSTC 111-2112-M-003-006,
NSTC 112-2112-M-003-007, and NSTC 112-2811-M-003-004-.

## References

  * (1) S. Dichiara, J. D. Gropp, J. A. Kennea, et al. (for the Neil Gehrels Swift Observatory Team), GCN Circ., No. 32632 (2022).
  * (2) P. Veres, E. Burns, E. Bissaldi, et al. (for the Fermi- GBM team), GCN Circ., No. 32636 (2022).
  * (3) D. Dzhappuev,Yu.Afashokov,I.Dzaparova,etal.(Car- pet-2 Group), Astron. Telegram, No. 15669 (2022).
  * (4) Y. Huang, Sh. Hu, S. Chen, et al. (for the LHAASO Experiment), GCN Circ., No. 32677 (2022)
  * (5) Z. Cao et al. [LHAASO], Science 380, no.6652, adg9328 (2023) doi:10.1126/science.adg9328 [arXiv:2306.06372 [astro-ph.HE]].
  * (6) R. Gould and G. Schréder, Phys. Rev. Lett. 16, no.6, 252-254 (1966) doi:10.1103/PhysRevLett.16.252
  * (7) G. G. Fazio and F. W. Stecker, Nature 226, 135-136 (1970) doi:10.1038/226135a0
  * (8) A.Nikishov R. J. Protheroe and H. Meyer, Phys. Lett. B 493, 1-6 (2000) doi:10.1016/S0370-2693(00)01113-8 [arXiv:astro-ph/0005349 [astro-ph]].
  * (9) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440-1443 (1977) doi:10.1103/PhysRevLett.38.1440
  * (10) L. Maiani, R. Petronzio and E. Zavattini, Phys. Lett. B 175, 359-363 (1986) doi:10.1016/0370-2693(86)90869-5
  * (11) G. Raffelt and L. Stodolsky, Phys. Rev. D 37, 1237 (1988) doi:10.1103/PhysRevD.37.1237
  * (12) G. Zhang and B. Q. Ma, Chin. Phys. Lett. 40, no.1, 011401 (2023) doi:10.1088/0256-307X/40/1/011401 [arXiv:2210.13120 [hep-ph]].
  * (13) H. Davoudiasl, P. Langacker and M. Perelstein, Phys. Rev. D 65, 105015 (2002) doi:10.1103/PhysRevD.65.105015 [arXiv:hep-ph/0201128 [hep-ph]].
  * (14) A. Esmaili, O. L. G. Peres and Z. Tabrizi, JCAP 12, 002 (2014) doi:10.1088/1475-7516/2014/12/002 [arXiv:1409.3502 [hep-ph]].
  * (15) C. Deffayet and J. P. Uzan, Phys. Rev. D 62, 063507 (2000) doi:10.1103/PhysRevD.62.063507 [arXiv:hep-ph/0002129 [hep-ph]].
  * (16) N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Lett. B 429, 263-272 (1998) doi:10.1016/S0370-2693(98)00466-3 [arXiv:hep-ph/9803315 [hep-ph]].
  * (17) I. Antoniadis, N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Lett. B 436, 257-263 (1998) doi:10.1016/S0370-2693(98)00860-0 [arXiv:hep-ph/9804398 [hep-ph]].
  * (18) K. R. Dienes, E. Dudas and T. Gherghetta, Phys. Rev. D 62, 105023 (2000); L. Di Lella, A. Pilaftsis, G. Raffelt and K. Zioutas, Phys. Rev. D 62, 125011 (2000)
  * (19) S. Chang, S. Tazawa and M. Yamaguchi, Phys. Rev. D 61, 084005 (2000) doi:10.1103/PhysRevD.61.084005 [arXiv:hep-ph/9908515 [hep-ph]].
  * (20) G. Galanti, L. Nava, M. Roncadelli, F. Tavecchio and G. Bonnoli, [arXiv:2210.05659 [astro-ph.HE]]; A. Baktash, D. Horns and M. Meyer, [arXiv:2210.07172 [astro-ph.HE]]; S. V. Troitsky, Pisma Zh. Eksp. Teor. Fiz. 116, no.11, 745-746 (2022); S. Nakagawa, F. Takahashi, M. Yamada and W. Yin, Phys. Lett. B 839, 137824 (2023); C. Eckner and F. Calore, Phys. Rev. D 106, no.8, 083020 (2022); C. Dessert, D. Dunsky and B. R. Safdi, Phys. Rev. D 105, no.10, 103034 (2022); L. Dirson and D. Horns, Astron. Astrophys. 671, A67 (2023)
  * (21) A. Esmaili, O. L. G. Peres and Z. Tabrizi, JCAP 12, 002 (2014); R. Barbieri, P. Creminelli and A. Strumia, Nucl. Phys. B 585, 28-44 (2000)
  * (22) J. Granot, T. Piran, O. Bromberg, J. L. Racusin and F. Daigne, Space Sci. Rev. 191, no.1-4, 471-518 (2015) doi:10.1007/s11214-015-0191-6 [arXiv:1507.08671 [astro-ph.HE]].
  * (23) A. MacFadyen and S. E. Woosley, Astrophys. J. 524, 262 (1999) doi:10.1086/307790 [arXiv:astro-ph/9810274 [astro-ph]].
  * (24) L. Wang and B. Q. Ma, Phys. Rev. D 108, no.2, 023002 (2023) doi:10.1103/PhysRevD.108.023002 [arXiv:2304.01819 [astro-ph.HE]].
  * (25) J. G. Lee, E. G. Adelberger, T. S. Cook, S. M. Fleischer and B. R. Heckel, Phys. Rev. Lett. 124, no.10, 101101 (2020) doi:10.1103/PhysRevLett.124.101101 [arXiv:2002.11761 [hep-ex]].


# Toy models of multibranched polymers: opened vs. circular structures

K. Haydukivska1,2 and V. Blavatska1 1Institute for Condensed Matter Physics of
the National Academy of Sciences of Ukraine 79011 Lviv, Ukraine  2 Institute
of Physics, University of Silesia, 75 Pułku Piechoty 1, 41-500 Chorzów, Poland
<EMAIL_ADDRESS>
###### Abstract

We study the conformational properties of complex Gaussian polymers containing
$f_{c}$ linear branches and $f_{r}$ closed loops, periodically tethered at $n$
branching points to either a linear polymer backbone (generalized bottlebrush
structures) or closed polymer ring (decorated ring structure). Applying the
path integration method, based on Edwards continuous chain model, we obtain in
particular the exact values for the size ratios comparing the gyration radii
of considered complex structures and linear chains of the same total molecular
weight, as functions of $n$, $f_{c}$ and $f_{r}$. Compactification of the
overall effective size of branched macromolecules with the increasing number
of loops is quantitatively confirmed. Our results are supported by numerical
estimates obtained by application of Wei’s method.

###### pacs:

36.20.Fz, 33.15.Bh, 87.15.hp

††: Journal of Physics A: Mathematical and Theoretical

##  1 Introduction

Bottlebrush polymers are important representatives of the class of complex
multibranched structures, where the sets of $f_{c}$ side chains are
periodically tethered at $n$ branching points to a linear polymer backbone
[1]. The properties of such molecules are governed by competition of the
steric repulsion between the side chains and the configurational entropy of
the main backbone chain. In particular, this leads to enhanced stiffness of
the backbone which causes the considerable elongation of molecule in solvent

[2, 3]; the side-chain molecules attend the extended cylindrical shapes in
solvents [4, 5]. In turn, the molecular entanglement in melts of such
molecules is significantly suppressed and their conformational degrees of
freedom are lowered [6, 1, 7], which dictates the viscoelastic properties of
such molecules. These peculiarities can be processed and controlled in a wide
range by varying the structural parameters such as grafting density of the
side chains, their lengths and position along the backbone.

Figure 1:  Schematic representation of the open (left) and circular (right)
bottlebrush polymer containing $f_{c}$ linear side chains and $f_{r}$ side
chains in form of closed loops.

Loop formation in macromolecules plays an important role in a number of
biochemical processes such as stabilisation of globular proteins [8, 9], DNA
compactification in the nucleus [10, 11, 12] etc. Cyclic (or ring)
macromolecules in general exhibit a number of unique properties in comparison
to their linear analogues. In particular, they are characterized by a more
compact coil conformation due to their lower conformational degree of freedom
and thus diffuse much faster in melts than linear polymer chains of the same
molecular weight, thus exhibiting very different viscoelastic and phase
transition behavior [13, 14]. Cyclic polymers are thus finding their
application as viscosity modifiers, since their addition in small amounts to
linear matrices increases the viscosity [15, 16]. Numerous analytical and
numerical studies have been conducted to analyse the cyclisation probability
and loop size distributions in long flexible macromolecules [17, 18, 19, 20,
21, 22, 23, 24, 25]. In particular, the so-called rosette polymers, containing
a single branching point with $f_{c}$ linear chains and $f_{r}$ circular loops
eliminating have been studied recently in our previous works [26, 27]. In this
concern, it is worthwhile to consider the generalised multibranched
structures, where some amount $f_{r}$ of side chains at each the branching
point are closed and form loops (see Fig. 1). One more closely related
architecture of a branched polymer is a so-called decorated ring polymer: a
bottlebrush polymer that has its backbone looped (see Fig. 1). The cyclic
polystyrenes with multiple branching have been synthesized and analyzed
recently in Refs. [28, 29].

In our present study, we will focus on conformational properties of the above
mentioned complex polymer structures in the regime of dilute solution. In
particular, the size ratio of gyration radii of a branched polymer $\langle
R_{g}^{2}\rangle_{\rm branched}$ to that of the linear chain $\langle
R_{g}^{2}\rangle_{\rm chain}$ of the the same total molecular weight

$g=\frac{\langle R_{g}^{2}\rangle_{\rm branched}}{\langle
R_{g}^{2}\rangle_{\rm chain}}$ (1)

is often used to characterize the impact of the complex topology of structure
on the effective elongation in solvent (it is also referred to as the
shrinking factor). Here and below, $\langle(\ldots)\rangle$ denotes averaging
over an ensemble of possible polymer conformations. Originally, this value was
introduced by Zimm [30] and evaluated for the so-called star polymer
structures, containing only one branching point, as well as for a single
closed ring. Another important quantity that characterizes the effective size
of a polymer coil is hydrodynamic radius $R_{h}$, which is directly obtained
in dynamic light scattering experiments. To compare $R^{2}_{g}$ and
$R^{-1}_{h}$, it is convenient to introduce the universal size ratio

$\rho=\sqrt{\langle R_{g}^{2}\rangle}/\langle R_{H}\rangle.$ (2)

The exact analytical result for the ratio (2) in $d=3$ dimensions have been
evaluated for the linear chain [30, 31, 32] and single closed ring [31, 33,
34]. Note, that the above mentioned quantities can be estimated exactly in the
case of ideal Gaussian polymers without taking into account an excluded volume
effect. Recently, both the analytical and numerical estimates for the
universal ratios (1), (2) have been obtained in Refs. [26, 27, 35] for so-
called rosette structure, containing single branching point with $f_{r}$
closed rings and $f_{c}$ linear chains, which can be considered as bottlebrush
polymers with $n=2$. Note that the probability distribution functions of the
gyration radii of Gaussian star and rosette structures have been derived
analytically in Ref. [35] as well. In this approach, an estimate for the
Gaussian bottlebrush polymers containing $n$ branching points of functionality
$f$ is found in Ref. [36]. There is a number of studies dedicated to
conformational properties (and in particular the size ratios) of bottlebrush
polymers with linear side chains [37, 38, 40, 39, 41, 42]. In particular, a
thorough comparative study of size properties of bottlebrush polymer
structures with linear, star and ring polymers is performed in Ref. [39]. The
aim of the present study is to continue these investigations to the case when
the side chains form the closed loops. This problem is closely related with
the process called loop extrusion in chromatin [43, 44, 45, 46], serving as
one of the mechanisms of compaction of long DNA molecules in nuclei. Moreover,
we will consider the case when the backbone chain form a closed loop (the case
of “decorated ring”), which is of great interest in particular since DNA often
has a closed ring shape [47, 48].

The layout of the paper is as follows. In Section 2 we introduce the
continuous chain model for both structures that is followed by analytical
calculation in Section 3. The Wei’s method is introduced in Section 4 that is
followed by Results and Discussions is Section 5. We finish this work with
concluding remarks in Section 6.

##  2 The model

An analytical description is performed in the framework of continuous chain
model. Here, the individual polymer chain is presented as a trajectory of
length $L$ parameterised by a radius vector $\vec{r}(s)$ with $s$ changing
from $0$ to $L$. The trajectories connectivity is determined by the
hamiltonian[49]:

$H^{0}=\frac{1}{2}\int_{0}^{L}\left(\frac{d\vec{r}(s)}{ds}\right)^{2}.$ (3)

We consider a generalized case when each of the $n-1$ branching points on the
backbone linear chain can contain $f_{c}$ side chains and $f_{r}$ rings.
Mathematically this topology is defined by the partition function:

$\displaystyle Z_{bb}=\frac{1}{Z_{0}}\int
D\vec{r}\prod_{k=1}^{n-1}\delta(\vec{r_{k}}(L_{b})-\vec{r_{k+1}}(0))$
$\displaystyle\times\prod_{i=1}^{f_{c}+f_{r}}\delta(\vec{r_{k}}(L)-\vec{r_{i}}(0))\prod_{j=1}^{f_{r}}\delta(\vec{r_{j}}(L)-\vec{r_{j}}(0)){\rm
e}^{-\sum_{j=1}^{F}H_{j}^{0}}.$ (4)

Here, the first set of $\delta$-functions describes the backbone that consists
of $n$ segments with each consecutive segment starting at the end of the
previous one. The second and the third sets of $\delta$-functions stands for
the fact that each of the branching points serves as a starting point for
$f_{c}+f_{r}$ outgoing trajectories. Note that here we consider the case when
all the side chains are of the same length $L$, however the segments of the
backbone are of the different length $L_{b}$ such that
$\lim_{L\rightarrow\infty}(L_{b}/L)=a$. In this case the total length of the
backbone is $nL_{b}$. The total number of trajectories for this structure
being $F=n+(n-1)(f_{c}+f_{r})$.

For the case of decorated ring, that is a bottlebrush polymer that has two end
points of its backbone chain looped, the partition function is presented as:

$\displaystyle Z_{dr}=\frac{1}{Z_{0}}\int
D\vec{r}\,\,\delta(\vec{r_{n}}(L_{b})-\vec{r_{0}}(0))$
$\displaystyle\times\prod_{k=0}^{n-1}\delta(\vec{r_{k}}(L_{b})-\vec{r_{k+1}}(0))\prod_{i=1}^{f_{c}+f_{r}}\delta(\vec{r_{k}}(L)-\vec{r_{i}}(0))$
$\displaystyle\times\prod_{j=1}^{f_{r}}\delta(\vec{r_{j}}(L)-\vec{r_{j}}(0)){\rm
e}^{-\sum_{j=1}^{F}H_{j}^{0}}.$ (5)

Here, an additional condition $\delta(\vec{r_{n}}(L)-\vec{r_{0}}(0))$ stands
for its end points looping. The total number of trajectories in this case is
$F=n(f_{c}+f_{r}+1)$.

Within this model, any observable of interest can be calculated as an average
over an ensemble of all possible trajectories:

$\langle(\ldots)\rangle=\frac{1}{Z}\int D\vec{r}{\rm
e}^{-\sum_{j=1}^{F}H_{j}^{0}},$ (6)

where $Z$ is a partition function given either by Eq. (4) or Eq. (5).

##  3 Gaussian approximation for analyzing the size characteristics of
polymers

###  3.1 Radius of gyration and hydrodynamic radius definitions and
calculation strategies

In terms of continuous chain model introduced above, the size measures of
polymer structure can be defined as:

$\displaystyle\langle
R^{2}_{g}\rangle=\frac{1}{2(FL)^{2}}\sum_{i,j=1}^{F}\int_{0}^{L}\int_{0}^{L}\langle(\vec{r_{i}}(s_{2})-\vec{r_{j}}(s_{2}))^{2}\rangle,$
(7) $\displaystyle\langle
R^{-1}_{h}\rangle=\frac{1}{(FL)^{2}}\sum_{i,j=1}^{F}\int_{0}^{L}\int_{0}^{L}\langle|\vec{r_{i}}(s_{2})-\vec{r_{j}}(s_{2})|^{-1}\rangle,$
(8)

here $s_{2}$ and $s_{1}$ stand for positions along and referred as
“restriction points”.

The calculation of the gyration radius within this model is performed by
making use of identity [49]:

$\displaystyle\langle(\vec{r}_{i}(s_{2})-\vec{r}_{j}(s_{1}))^{2}\rangle=-2\frac{d}{d|\vec{k}|^{2}}\xi(\vec{k})_{\vec{k}=0},$
$\displaystyle\xi(\vec{k})\equiv\langle{\rm
e}^{-\iota\vec{k}(\vec{r}_{i}(s_{2})-\vec{r}_{j}(s_{1}))}\rangle.$ (9)

The crucial point in the calculation of the hydrodynamic radius is exploiting
the following equality [50]:

$|\vec{r}|^{-1}=(2\pi)^{-d}\int\,d\vec{k}\,2^{d-1}\pi^{\frac{d-1}{2}}\Gamma\left(\frac{d-1}{2}\right)\,k^{1-d}{\rm
e}^{i\vec{r}\vec{k}}.$ (10)

With this Fourier transformation the average in the integrand in expression
(8) can be rewritten as:

$\displaystyle\langle|\vec{r}_{i}(s_{2})-\vec{r}_{j}(s_{1})|^{-1}\rangle=(2\pi)^{-d}\int\,d\vec{k}\,2^{d-1}\pi^{\frac{d-1}{2}}\Gamma\left(\frac{d-1}{2}\right)\,k^{1-d}\xi(\vec{k}).$
(11)

Note that in both identities the averaging is performed for the same
expression $\langle{\rm
e}^{i\vec{k}(\vec{r}_{i}(s_{2})-\vec{r}_{j}(s_{1}))}\rangle$. In the Gaussian
case, performing the averaging according to (6) for a single polymer chain of
length $L$ leads to an expression:

$\langle{\rm e}^{i\vec{k}(\vec{r}_{i}(s_{2})-\vec{r}_{j}(s_{1}))}\rangle={\rm
e}^{-\frac{\vec{k}^{2}}{2}(s_{2}-s_{1})},$ (12)

and the application of the identities (9) and (11) results in:

$\displaystyle\langle(\vec{r}_{i}(s_{2})-\vec{r}_{j}(s_{1}))^{2}\rangle=d(s_{2}-s_{1}),$
(13)
$\displaystyle\langle|\vec{r}_{i}(s_{2})-\vec{r}_{j}(s_{1})|^{-1}\rangle=\frac{\Gamma\left(\frac{d-1}{2}\right)}{\sqrt{2}\Gamma\left(\frac{d}{2}\right)}(s_{2}-s_{1})^{-1/2}.$
(14)

After performing integrating over $s_{1}$ and $s_{2}$ and dividing by $L^{2}$
the final expressions for a single polymer chain are:

$\displaystyle\langle R^{2}_{g}\rangle_{{\rm chain}}=\frac{dL}{6},$ (15)
$\displaystyle\langle R^{-1}_{h}\rangle_{{\rm
chain}}=\frac{\sqrt{2}\Gamma\left(\frac{d-1}{2}\right)}{\Gamma\left(\frac{d}{2}\right)}\frac{4}{3}L^{-\frac{1}{2}}.$
(16)

The size ratio between those two size characteristics as defined by expression
(2) is given by:

$\rho_{{\rm chain}}=\frac{\sqrt{\langle R^{2}_{g}\rangle}_{{\rm
chain}}}{\langle R^{-1}_{h}\rangle_{{\rm
chain}}^{-1}}=\frac{\sqrt{2d}\Gamma\left(\frac{d-1}{2}\right)}{\Gamma\left(\frac{d}{2}\right)}\frac{4}{3\sqrt{6}},$
(17)

which reproduces a well-known value [51].

###  3.2 Calculation of size characteristics

Figure 2:  Schematic representation of the contributions into the gyration and
hydrodynamic radius of rosette polymer. Here, solid line represents polymer
trajectory and bullets represent restriction points.

It is convenient to use a diagram technique for schematic presentations of
contribution into $\xi(\vec{k})$. In the case of bottlebrush structure, the
set of diagrams consists of five diagrams given on Fig. 2 (those corresponding
to rosette structures and being calculated previously in our works [26, 27])
and three additional ones presented on Fig. 3. Here, the length of the each
segment plays an important role and has to be accounted for. The analytical
expressions corresponding to all the diagrams are provided in the Appendix.

Figure 3:  Schematic representation of the contributions into the gyration and
hydrodynamic radius of bottlebrush polymer. Notations are the same like on
Fig. 2.

The size characteristics of interest can in general be presented as:

$\displaystyle
R_{x}=\frac{L^{p}}{(an+(f_{r}+f_{c})(n-1))^{2}}\left((n-1)\left[f_{c}\xi_{1,x}(1)+f_{r}\xi_{2,x}(1)\right.\right.$
$\displaystyle\left.\left.+\frac{f_{c}(f_{c}-1)}{2}\xi_{3,x}(1,1)+\frac{f_{r}(f_{r}-1)}{2}\xi_{5,x}(1,1)+f_{r}f_{c}\xi_{4,x}(1,1)\right.\right.$
$\displaystyle\left.+2f_{c}\xi_{3,x}(1,a)+\xi_{3,x}(a,a)+2f_{r}\xi_{4,x}(1,a)\right]+n\xi_{1,x}(a)$
$\displaystyle+\left.\sum_{i=1}^{n-2}(n-i-1)\left[f_{c}^{2}\xi_{7,x}(1,1)+2f_{c}\xi_{7,x}(1,a)+\xi_{7,x}(a,a)+2f_{c}f_{r}\xi_{8,x}(1,1)\right.\right.$
$\displaystyle\left.\left.+2f_{r}\xi_{8,x}(1,a)+f_{r}^{2}\xi_{6,x}(1,1)\right]\right).$
(18)

Here, $p$ is equal to $1$ for the gyration radius and $-1/2$ for the
hydrodynamic radius and $x$ stands to the indexes $g$ or $h$. In the case of
gyration radius the summations are easily performed and the following
expression results:

$\displaystyle\langle
R^{2}_{g}\rangle_{bb}=\frac{dL}{(an+(f_{r}+f_{c})(n-1))^{2}}\left(\frac{a^{3}n^{3}}{6}\right.$
$\displaystyle+\left.\frac{n-1}{6}\left((an^{2}-2an+3n-3)f_{c}^{2}+(2a^{2}n^{2}-a^{2}n+3an-2)f_{c}\right.\right.$
$\displaystyle\left.\left.+\frac{1}{2}(4a^{2}n^{2}-2a^{2}n+2an-1)f_{r}+(an^{2}-2an+n-1)f_{r}^{2}\right.\right.$
$\displaystyle\left.\left.+2(an^{2}-2an+2n-2)f_{c}f_{r}\right)\right).$ (19)

With $f_{r}=0$, the above expression restores that for the bottlebrush polymer
with linear side chains, derived in Ref. [36]. To evaluate an estimate for the
size ratio (1), we divide this expression by the value of gyration radius of a
single chain with the same total molecular weight, that is $\langle
R^{2}_{g}\rangle_{chain}=\frac{dL(an+(f_{r}+f_{c})(n-1))}{6}$. The graphic
representation of resulting expression for some values of $n,f_{c},f_{r}$ are
given on Fig. 6.

Correspondingly, we evaluated an expression for hydrodynamic radius:

$\displaystyle\langle
R^{-1}_{h}\rangle_{bb}=\frac{\sqrt{\frac{2}{L}}\Gamma\left(\frac{d-1}{2}\right)}{(f_{c}+f_{r})^{3}\Gamma\left(\frac{d}{2}\right)}\left(\frac{4a^{\frac{3}{2}}n}{3}+\frac{n-1}{12}\left((16(2a^{\frac{3}{2}}+f_{c}^{2}))(\sqrt{2}-1)\right.\right.$
$\displaystyle+(32(1+a)^{\frac{3}{2}}-32a^{\frac{3}{2}}-16\sqrt{2})f_{c}+12\arctan\left(\frac{1}{2\sqrt{a}}\right)(4a+1)f_{r}$
$\displaystyle\left.-6f_{r}\pi(f_{r}-1)(\sqrt{2}-2)+24\sqrt{a}f_{r}+3f_{r}f_{c}(4-\pi+10\arctan\left(\frac{1}{2}\right))\right)$
$\displaystyle+\frac{4}{3}\sum_{i=1}^{n-2}(n-i-1)\left(f_{c}^{2}(ai+2)^{\frac{3}{2}}+2f_{c}(ai+a+1)^{\frac{3}{2}}-2f_{c}(ai+1)^{\frac{3}{2}}(f_{c}+1)\right.$
$\displaystyle+(ai+2a)^{\frac{3}{2}}-2(a(i+1))^{\frac{3}{2}}(f_{c}+1)\left.+(ai)^{\frac{3}{2}}(f_{c}-1)^{2}\right)$
$\displaystyle+2f_{c}f_{r}\sum_{i=1}^{n-2}(n-i-1)\left(\frac{1}{2}\arctan\left(\frac{1}{2\sqrt{ai+1}}\right)(4ai+5)\right.$
$\displaystyle\left.-\frac{1}{2}\arctan\left(\frac{1}{2\sqrt{ai}}\right)(4ai+1)\right.\left.+\sqrt{ai+1}-\sqrt{ai}\right)$
$\displaystyle+2f_{r}\sum_{i=1}^{n-2}(n-i-1)\left(\frac{1}{2}\arctan\left(\frac{1}{2\sqrt{ai+a}}\right)(4ai+4a+1)\right.$
$\displaystyle\left.-\frac{1}{2}\arctan\left(\frac{1}{2\sqrt{ai}}\right)(4ai+1)+\sqrt{ai+a}-\sqrt{ai}\right)$
$\displaystyle-
f_{r}^{2}\sum_{i=1}^{n-2}\frac{n-i-1}{\sqrt{4ai+2}}\left((4ai+2)\left(\arctan\left(\frac{4ai+1+\sqrt{4ai+2}}{2\sqrt{ai}}\right)\right.\right.$
$\displaystyle\left.\left.+\arctan\left(\frac{\sqrt{4ai+2}-4ai-1}{2\sqrt{ai}}\right)\right)\right.-2\arcsin\left(\frac{1}{\sqrt{4ai+1}}\right)\sqrt{4ai+2}$
$\displaystyle\left.-2\arctan\left(\frac{1}{2\sqrt{ai}}\right)\sqrt{4ai+2}\right).$
(20)

Here $L$ is the length of the side chain and $L_{b}=aL$.

The ratio (2) for the bottlebrush structures can be evaluated next on the base
of (19) and (20). The graphic representation of resulting expression for some
values of $n,f_{c},f_{r}$ are presented on Fig. 8.

For the case of decorated ring structures, we again need to consider three
additional diagrams (see fig. 4) besides of those already considered for
bottlebrush structures.

Figure 4:  Schematic representation of contributions into the gyration and hydrodynamic radii of decorated ring polymer. Notations are the same like on Fig. 2. $f_{c}$ | $f_{r}$ | $n$ | $\rho$
---|---|---|---
$0$ | $1$ | $2$ | $\frac{\sqrt{22}}{32\sqrt{\pi}}\left(3\pi\sqrt{2}-2\pi\sqrt{6}-4\arctan(2\sqrt{2})+2\arcsin\left(\frac{\sqrt{3}}{3}\right)\right.$
|  |  | $\left.+9\pi-6\arctan\left(\frac{5}{\sqrt{2}}\right)\right)$
$0$ | $2$ | $2$ | $\frac{\sqrt{6}}{36\sqrt{\pi}}\left(3\pi\sqrt{2}-4\pi\sqrt{6}+30\pi-16\arctan(2\sqrt{2})\right.$
|  |  | $\left.-24\arctan\left(\frac{5}{\sqrt{2}}\right)+8\arcsin\left(\frac{\sqrt{3}}{3}\right)\right)$
$0$ | $1$ | $3$ | $\frac{\sqrt{69}}{216\sqrt{\pi}}\left(21\pi-3\pi\sqrt{3}-2\arctan\left(\frac{24\sqrt{7}}{103}\right)\sqrt{42}+12\arctan\left(\frac{\sqrt{6}}{4}\right)\right.$
|  |  | $\left.+12\arcsin\left(\frac{\sqrt{3}}{\sqrt{11}}\right)\right)$
$0$ | $2$ | $3$ | $\frac{\sqrt{74/\pi}}{486}\left(3\pi\left(18-2^{\frac{3}{2}}-3^{\frac{3}{2}}\right)-8\arctan\left(\frac{24\sqrt{7}}{103}\right)\sqrt{42}\right.$
|  |  | $\left.+48\left(\arctan\left(\frac{\sqrt{6}}{4}\right)+\arcsin\left(\frac{\sqrt{3}}{\sqrt{11}}\right)\right)\right)$
$1$ | $0$ | $2$ | $\frac{\sqrt{35}}{96\sqrt{\pi}}\left(5\sqrt{10}+\sqrt{2}+20-6\sqrt{6}+18\sqrt{2}\arcsin\left(\frac{\sqrt{3}}{3}\right)\right)$
Table 1: Size ratio (2) for a set of decorated rings with fixed parameters
$n$, $f_{r}$, $f_{c}$, assuming $a=1$.

Similarly to Eq.(18), the general expression for size characteristics reads:

$\displaystyle
R_{x}=\frac{L^{p}}{((f_{r}+f_{c}+a)n)^{2}}\left(n\left(f_{c}\xi_{1,x}(1)+f_{r}\xi_{2,x}(1)\right.\right.$
$\displaystyle\left.+\frac{f_{c}(f_{c}-1)}{2}\xi_{3,x}(1,1)+\frac{f_{r}(f_{r}-1)}{2}\xi_{5,x}(1,1)+f_{r}f_{c}\xi_{4,x}(1,1)\right)$
$\displaystyle+\xi_{1,x}(na)+f_{c}n\xi_{4,x}(1,an)+f_{r}n\xi_{5,x}(1,an)$
$\displaystyle+\left.\sum_{i=1}^{n-1}(n-i-1)\left(f_{c}^{2}\xi_{9,x}(1,1)+2f_{r}f_{c}\xi_{11,x}(1,1)+f_{r}^{2}\xi_{10,x}(1,1)\right)\right).$
(21)

For the gyration radius of decorated ring structure we thus have

$\displaystyle\langle
R^{2}_{g}\rangle_{dr}=\frac{dL}{(an+(f_{r}+f_{c})n)^{2}}\left(\frac{a^{3}n^{3}}{12}+\frac{n}{12}\left((an^{2}-a+6n)f_{c}^{2}\right.\right.$
$\displaystyle\left.\left.+(2a^{2}n^{2}+6an-4)f_{c}+(2a^{2}n^{2}+2an-1)f_{r}\right.\right.$
$\displaystyle\left.\left.+(an^{2}-a+2n)f_{r}^{2}+(2an^{2}-2a+8n-8)f_{c}f_{r}\right)\right).$
(22)

Again, to evaluate the size ratio (1), we divide this expression by $\langle
R^{2}_{g}\rangle_{chain}$ of the same total molecular weight, and give some
graphical representations on Fig. 6.

The expressions for hydrodynamic radius of decorated ring structure is much
more cumbersome and we give here examples of this expression for several
simple cases of $n,f_{c}$ and $f_{r}$ in Table 1.

##  4 Numerical results: Wei’s method

To evaluate the size ratio (1) of multibranched structures as compared with
linear chain, we also applied the method originally developed by Wei [53]. As
the first step, we present the polymer structure as a mathematical graph,
containing $M$ junction points (e.g. for a bottlebrush polymer with $n$
branching points of functionality $f_{c}$ one has $M=nf_{c}+2$, where two free
ends of backbone are included into consideration as junction points with
functionality $1$). The Kirchhoff matrix ${\bf K}$ of the size $M\times M$ of
the given graph is then defined as follows. Its diagonal elements $K_{ii}$
equal the degree of vertex $i$, whereas the non-diagonal elements $K_{ij}$
equal $-1$ when the vertices $i$ and $j$ are adjacent and $0$ otherwise. Let
$\lambda_{2},\ldots,\lambda_{M}$ be $(M-1)$ non-zero eigenvalues of the
$M\times M$ Kirchhoff matrix ($\lambda_{1}$ is $0$). The $g$-ratio of the
radii of gyration of a given branched structure and that of a linear chain
with the same length is given by:

$g=\frac{\sum_{j=2}^{M}1/\lambda_{j}^{{\rm
network}}}{\sum_{j=2}^{M}1/\lambda_{j}^{{\rm chain}}},$ (23)

where $\lambda_{j}^{{\rm network}}$ and $\lambda_{j}^{{\rm chain}}$ are
network and chain Kirchhoff matrix eigenvalues.

Figure 5: The size ratio of the bottlebrush (filled squares) and decorated
ring (open squares) at $n=3$, $f_{c}=2$, $f_{r}=0$, obtained in numerical
simulations with application of Wei formula (23), as functions of a single
link length $l$.

This quantity is determined within the Wei’s approach by extrapolating to
infinite length of links $l$, so that the total number of junction points in
resulting graph (and the size of corresponding Kirchhof matrix) is
$M=l((f_{c}+2f_{r})n+n+2)$. In our analysis, we considered the structures with
number of $l$ up to 20. The set of eigenvalues of corresponding Kirchhof
matrix was evaluated for the structures presented on Fig. 1 at different $n$
and values of $g$ were estimated according to Eq. (23). Figure 5 presents
examples of simulation data for the case $n=3$ and $f_{c}=2$, $f_{r}=0$. For
the finite link length $l$, the values of parameters differ from those for
infinitely long polymer chains and the finite-size deviation can be fitted by:

$\displaystyle g(l)=g+b/l,$ (24)

where $a$, $b$ are constants. The size ratio estimates were obtained by us by
least-square fitting of (24). These are presented in Fig. 6 and compared with
analytical results. Note also, that we found this method to be less reliable
with increasing the number of internal closed loops in structures, due to
emergence of negative eigenvalues in corresponding graphs [55].

##  5 Results and Discussion

Figure 6:  The size ratio (1) for bottlebrush polymers (left) and decorated
rings (right) as functions of $n$ with two different sets of side chain types.

Figure 7:  The size ratios (1) on the left and (2) on the right as functions
of $n$ for toothbrushes and decorated rings with only one closed side chain.

We start by considering a size ratio (1) of the gyration radius of branched
polymers with $f_{c}$ linear branches and $f_{r}$ branches in form of closed
loops at each of branching points, as given by (19) and (22) for open and
closed backbones correspondingly, and that of a linear chain of the same total
molecular weight. Here we limit the discussion to the case with $a=1$ that
corresponds to the backbone segment between the branching points being of the
same length as the side chains. Note that at $n=3$ and $f=f_{c}+1,f_{r}=0$ we
restore the corresponding expression for the so-called pom-pom polymer
architecture [30, 52], containing two branching points each of functionality
$f$. As examples, the values at some fixed parameters are presented on Fig. 6
in comparison with our numerical data, obtained by application of Wei’s
method.

At any value of $n$, the value of $g$ is smaller than $1$, which describes the
shrinking of the branched structure in comparison with the linear one. This
effect is increasing with increasing $n$. In the case of decorated ring (Fig.
6,right) this compactification of the effective structure size at each $n$ is
more pronounced. Furthermore, presence of side loops in turn leads to
shrinking of the effective size of considered structures as compared with
effect caused by linear side chains. Another important point to note is that
both analytical and numerical methods provide results that are in a good
agreement, though the path integration method allows to obtained the exact
values of parameters for the case of Gaussian polymers, while the Wei’s method
provides only approximate values.

Separately we considered the case $f_{c}=0$, $f_{r}=1$ (see Fig. 7), which may
be related to a problem of so-called loop extrusion (formation a set of single
loops along the fibre) in chromatin [43, 44, 45, 46]. Loop extrusion is the
compaction process, that organizes DNA. An expected compactification
(decreasing of the radii of gyration as compared with linear fiber) is
observed.

Figure 8:  The size ratio (2) for bottlebrush polymers (left) and decorated
rings (right) as functions of $n$ with two different sets of side chain types.

An estimate for another important size ratio, defined by (2) have been
evaluated on the basis of expressions (19), (22), (20). As examples, the
estimates at some fixed values of parameters are presented on Fig. 8.
Interesting to note, that dependence on the number of branching points $n$ is
very different for bottlebrush structures with open and closed backbones. In
the first case, the size ratio decreases with $n$ at small $n$ and then start
to increase (left panel of 8). For the case of decorated rings, the behaviour
is crucially defined by the type of side branches: it decreases in presence of
linear side chains, while looped side chains stimulate the size ratio to
increase. For small values of $n$ and closed side chains, a much smaller ratio
in comparison with open side chains might be explained by more compact and
thus more dense structure, so that the difference between gyration and
hydrodynamic radii is very small.

##  6 Conclusions

The aim of the present study was to analyze the impact of the structural
parameters like the number of branching points $n$ and their functionalities
(number of linear $f_{c}$ and circular $f_{r}$ side branches) on the effective
size measures of bottlebrush and decorated rings polymers in solutions. In
particular, the averaged size of typical conformation of individual molecule
defines the diffusivity of such molecules in a solvent and thus impact their
viscoelastic and phase transition behavior.

Applying the path integration method based on Edwards model of continious
chain, we evaluated exact expressions for the gyration and hydrodynamic radii
of multibranched structures under considerations. The size ratios of obtained
gyration radii of branched structures (given by Eqs. (19) and (22)) and that
of the linear chain $\langle R_{g}^{2}\rangle_{\rm chain}$ of the the same
total molecular weight were used to characterize the impact of the complex
topology of structure on the effective size in solvent (the shrinking factor).
Our results quantitatively describe the shrinking of the branched structure in
comparison with the linear one; this effect is increasing with increasing $n$
and is more pronounced in the case of decorated ring. Presence of side loops
leads to further shrinking of the effective size of considered structures as
compared with linear side chains. These results are confirmed also by
numerical simulations performed by us within the frame of so-called Wei’s
method, which encounters some limitations however in presence of multiple
closed loops in a structure. Expressions for hydrodynamic radius of considered
structures is evaluated as well.

The subject of the forthcoming study will be to go beyond the simplified
Gaussian approximation by taking into account the excluded volume effect
between monomers. This, in particular, will allow us to evaluate the impact of
steric effects, caused by side branches in open and lopped form, on the
effective size measure of the backbone chain and to analyze the conformational
properties of considered structures more deeply.

##  Acknowledgements

Authors would like to acknowledge the support from the National Academy of
Sciences of Ukraine, Project KPKBK 6541230.

K.H. would like to acknowledge the support from the National Science Center,
Poland (Grant No. 2018/30/E/ST3/00428).

## Appendix

The analytical expressions corresponding to all the diagrams shown on Figs. 2,
3, 4, needed to calculate the gyration radius of considered polymer
structures, are listed below:

$\displaystyle\xi_{1,g}=d\frac{L^{3}}{6},$ (25)
$\displaystyle\xi_{2,g}=d\frac{L^{3}}{12},$ (26)
$\displaystyle\xi_{3,g}=d\frac{L_{1}L_{2}(L_{1}+L_{2})}{2},$ (27)
$\displaystyle\xi_{4,g}=d\frac{L_{1}L_{2}(3L_{1}+L_{2})}{6},$ (28)
$\displaystyle\xi_{5,g}=d\frac{L_{1}L_{2}(L_{1}+L_{2})}{6},$ (29)
$\displaystyle\xi_{6,g}=d\left(\frac{L_{1}L_{2}(L_{1}+L_{2})}{6}+iL_{b}L_{1}L_{2}\right),$
(30)
$\displaystyle\xi_{7,g}=d\left(\frac{L_{1}L_{2}(L_{1}+L_{2})}{2}+iL_{b}L_{1}L_{2}\right),$
(31)
$\displaystyle\xi_{8,g}=d\left(\frac{L_{1}L_{2}(3L_{1}+L_{2})}{6}+iL_{b}L_{1}L_{2}\right),$
(32)
$\displaystyle\xi_{9,g}=d\left(\frac{L_{1}L_{2}(L_{1}+L_{2})}{2}+iL_{b}L_{1}L_{2}+\frac{iL_{b}L_{1}L_{2}}{n}\right),$
(33)
$\displaystyle\xi_{10,g}=d\left(\frac{L_{1}L_{2}(L_{1}+L_{2})}{6}+iL_{b}L_{1}L_{2}+\frac{iL_{b}L_{1}L_{2}}{n}\right),$
(34)
$\displaystyle\xi_{11,g}=d\left(\frac{L_{1}L_{2}(3L_{1}+L_{2})}{6}+iL_{b}L_{1}L_{2}+\frac{iL_{b}L_{1}L_{2}}{n}\right).$
(35)

Here, $L,L_{1},\,L_{2},\,L_{b}$ are the lengths of the segments with $L_{b}$
being the length between the branching points. Note that diagrams $\xi_{1,g}$
and $\xi_{2,g}$ describe contributions from ether a chain or closed ring and
thus they depend on only one parameter; for diagrams $\xi_{4,g},\xi_{8,g}$ and
$\xi_{11,g}$, $L_{1}$ is the length of the looped trajectory with a
restriction point and $L_{2}$ of the open one (unlike the rest of the diagrams
these three are not symmetrical with respect to $L_{1}$ and $L_{2}$). The
number of backbone segments between the considered branching points is denoted
by $i$, thus changing from $1$ to $n$.

In the same way, we can present the list of expressions used to evaluate the
hydrodynamic radius:

$\displaystyle\xi_{1,h}=C_{\xi}\frac{4L^{3/2}}{3},$ (36)
$\displaystyle\xi_{2,h}=C_{\xi}\frac{\pi L^{3/2}}{2},$ (37)
$\displaystyle\xi_{3,h}=C_{\xi}\frac{4(L_{1}+L_{2})^{3/2}-4L_{1}^{3/2}-4L_{2}^{3/2}}{3},$
(38)
$\displaystyle\xi_{4,h}=C_{\xi}\left(L_{2}\sqrt{L_{1}}+\left(\frac{L_{2}^{3/2}}{2}+2\sqrt{L_{2}}L_{1}\right)\arcsin\left(\frac{L_{2}}{\sqrt{4L_{1}L_{2}+L_{2}^{2}}}\right)\right.$
$\displaystyle\left.-\frac{\pi L_{2}^{3/2}}{4}\right),$ (39)
$\displaystyle\xi_{5,h}=C_{\xi}\pi\left(L_{1}\sqrt{L_{2}}+L_{2}\sqrt{L_{1}}-\sqrt{L_{2}}\sqrt{L_{1}}\sqrt{L_{1}+L_{2}}\right),$
(40)
$\displaystyle\xi_{6,h}=\frac{C_{\xi}}{\sqrt{4L_{2}L_{b}i+L_{1}L_{2}+L_{2}^{2}}}\left(L_{2}\sqrt{L_{1}}\right.$
$\displaystyle\left(\arctan\left(\frac{4L_{b}i+L_{2}+\sqrt{L_{2}(4L_{b}i+L_{1}+L_{2})}}{2\sqrt{L_{1}L_{b}i}}\right)\right.$
$\displaystyle\left.-\arctan\left(\frac{4L_{b}i+L_{2}-\sqrt{L_{2}(L_{b}i+L_{1}+L_{2})}}{2\sqrt{L_{1}L_{b}i}}\right)\right)(4L_{b}i+L_{1}+L_{2})$
$\displaystyle\left.-2\sqrt{L_{2}(4L_{b}i+L_{1}+L_{2})}\left(L_{1}\sqrt{L_{2}}\arcsin\left(\frac{L_{2}}{\sqrt{L_{2}(4L_{b}i+L_{2})}}\right)\right.\right.$
$\displaystyle\left.\left.+L_{2}\sqrt{L_{1}}\arctan\left(\frac{\sqrt{L_{1}}}{2\sqrt{L_{b}i}}\right)\right)\right),$
(41)
$\displaystyle\xi_{7,h}=C_{\xi}\frac{4}{3}\left((L_{b}i)^{\frac{3}{2}}-(L_{b}i+L_{1})^{\frac{3}{2}}-(a+L_{2})^{\frac{3}{2}}+(L_{b}i+L_{1}+L_{2})^{\frac{3}{2}}\right),$
(42)
$\displaystyle\xi_{8,h}=C_{\xi}\frac{\sqrt{L_{1}}}{2}\left(\arctan\left(\frac{\sqrt{L_{1}}}{2\sqrt{L_{b}i+L_{2}}}\right)(4L_{b}i+L_{1}+4L_{2})\right.$
$\displaystyle\left.-\arctan\left(\frac{\sqrt{L_{1}}}{2\sqrt{L_{b}i}}\right)(4L_{b}i+L_{1})+2\sqrt{(L_{b}i+L_{2})L_{1}}-2\sqrt{L_{b}iL_{1}}\right),$
(43)
$\displaystyle\xi_{9,h}=\frac{4C_{\xi}}{3n^{\frac{3}{2}}}\left(((L_{1}+L_{2}+L_{b}i)n-L_{b}i^{2})^{\frac{3}{2}}-((L_{2}+L_{b}i)n-L_{b}i^{2})^{\frac{3}{2}}\right.$
$\displaystyle\left.-((L_{1}+L_{b}i)n-L_{b}i^{2})^{\frac{3}{2}}+(L_{b}i(n-i))^{\frac{3}{2}}\right),$
(44)
$\displaystyle\xi_{10,h}=\frac{C_{\xi}}{\sqrt{L_{2}n((4L_{b}i+L_{1}+L_{2})n-4L_{b}i^{2})}}\left(L_{2}\sqrt{(}L_{1})\right.\times$
$\displaystyle\left.\left(\arctan\left(\frac{4L_{b}i(n-i)-\sqrt{L_{2}}\sqrt{(4L_{b}i+L_{1}n+L_{2}n)-4L_{b}i^{2}}\sqrt{n}+L_{2}n}{2\sqrt{n(n-i)L_{1}iL_{b}}}\right)\right.\right.$
$\displaystyle\left.-\arctan\left(\frac{(4L_{b}i(n-i)+L_{2}n+\sqrt{L_{2}n((4L_{b}in+L_{1}+L_{2})-4L_{b}i^{2})}}{2\sqrt{L_{1}iL_{b}(n-i)n}}\right)\right)$
$\displaystyle\times((4L_{b}in+L_{1}+L_{2})-4L_{b}i^{2})$
$\displaystyle+2\sqrt{L_{2}n((4L_{b}i+L_{1}+L_{2})n-4L_{b}i^{2})}L_{1}L_{2}\times$
$\displaystyle\left.\left(\frac{1}{\sqrt{L_{1}}}\arcsin\left(\frac{\sqrt{L_{2}n}}{4L_{b}in+L_{2}n-4L_{b}i^{2}}\right)\right.\right.$
$\displaystyle\left.\left.+\frac{1}{\sqrt{L_{2}}}\arctan\left(\frac{\sqrt{L_{1}n}}{\sqrt{L_{b}i(n-i)}}\right)\right)\right),$
(45)
$\displaystyle\xi_{11,h}=\frac{C_{\xi}\sqrt{L_{2}}}{2n^{\frac{3}{2}}}\left(\sqrt{n}\arcsin\left(\frac{\sqrt{nL_{2}}}{\sqrt{(4L_{b}i+4L_{1}+L_{2})n-4L_{b}i^{2}}}\right)\right.$
$\displaystyle\times((4L_{b}i+4L_{1}+L_{2})n-4L_{b}i^{2})-\sqrt{n}\arcsin\left(\frac{\sqrt{nL_{2}}}{\sqrt{(4L_{b}i+L_{2})n-4L_{b}i^{2}}}\right)$
$\displaystyle\left.\times((4L_{b}i+L_{2})n-4L_{b}i^{2})+2\sqrt{L_{2}L_{b}in+L_{1}L_{2}n-L_{2}L_{b}i^{2}}n\right.$
$\displaystyle\left.-2\sqrt{(}L_{2}L_{b}in-L_{2}L_{b}i^{2})n\right).$ (46)

Here,
$C_{\xi}=\frac{\sqrt{2}\Gamma\left(\frac{d-1}{2}\right)}{\Gamma\left(\frac{d}{2}\right)}$.

## References

## References

  * [1] S.S. Sheiko, B.S. Sumerlin, and K. Matyjaszewski (2008) Prog. Polym. Sci. 33 759
  * [2] K. Terao, T. Hokajo, Y. Nakamura, and T. Norisuye (1999) Macromolecules 32 3690
  * [3] S. Kawaguchi, K. Akaike, Z.-M. Zhang, H. Matsumoto, and K. Ito (1998) Polym. J. 30 1004
  * [4] H. Liang, Z. Cao, Z. Wang, S.S. Sheiko, and A.V. Dobrynin (2017) Macromolecules 50 3430
  * [5] C.R. López-Barrón, P. Brant, A.P. Eberle, and D.J. Crowther (2015) J. Rheol. 59 865
  * [6] S. Namba, Y. Tsukahara, K. Kaeriyama, K. Okamoto, and M. Takahashi (2000) Polymer 41 5165
  * [7] J. Rzayev (2012) ACS Macro Lett. 1 1146
  * [8] C.N. Pace, G.R. Grimsley, J.A. Thomson, and B.J. Barnett (1988) J. Biol. Chem. 263 11820
  * [9] A.D. Nagi and L. Regan (1997) Folding Des. 2 67
  * [10] P. Fraser (2006) Curr. Opin. Genet. Dev. 16 490
  * [11] M. Simonis, P. Klous, E. Splinter, Y. Moshkin, R. Willemsen, E. de Wit, B. van Steensel, and W. de Laat, (2006) Nat. Genet. 38 1348
  * [12] J. Dorier and A. Stasiak (2009) Nucl. Acids Res. 37 6316
  * [13] M. Kapnistos, M. Lang, D. Vlassopoulos, W. Pyckhout-Hintzen, D. Richter, D. Cho, T. Chang,and M. Rubinstein (2008) Nat. Mat. 7 997\.
  * [14] D. Kawaguchi, K. Masuoka, A. Takano, K. Tanaka, Y. Nagamura, N. Torikai, R.M. Dalgliesh, S. Langridge, and Y. Matsushita (2006) Macromolecules 39 5180\.
  * [15] D. Vlassopoulos, R. Pasquino, and F. Snijkers, In Topological Polymer Chemistry; World Scientific Publishing: Singapore, (2014) 291
  * [16] G.B. Kenna and D.J. Plazek (1986) Polym. Commun. 27 304
  * [17] H.S. Chan and K.A. Dill(1989) J. Chem. Phys. 90 492
  * [18] A.M. Rubio, J.J. Freire, M. Bishop, and J.H.R. Clarke (1993) Macromolecules 26 4018
  * [19] S. Redner (1980) J. Phys. A 13 3525
  * [20] H.P. Hsu, W. Nadler, and P. Grassberger (2004) Macromolecules 37 4658
  * [21] B. Duplantier (1989) J. Stat. Phys. 54 581
  * [22] M. Wittkop, S. Kreitmeier, and D. Göritz (1996) J. Chem. Phys. 104 351
  * [23] H.-X. Zhou, J. Phys. Chem. (2001) 104 6763
  * [24] N.N. Toan, D. Marenduzzo, P.R. Cook, and C. Micheletti, Phys. Rev. Lett. (2006) 97 178302
  * [25] J. Shin, A. G. Cherstvy, and R. Metzler (2015) Soft Matter 11 472
  * [26] V. Blavatska and R. Metzler (2015) J. Phys. A: Math. Theor. 48 135001
  * [27] K. Haydukivska, V. Blavatska, and J. Paturej (2020) Sci. Rep. 10 14127
  * [28] M.D. Hossain, J.C. Reid, D. Lu, Z. Jia, D.J. Searles, and M.J. Monteiro (2018) Biomacromolecules bf 19 616
  * [29] M.D. Hossain, Z.F.Jia, and M.J. Monteiro (2014) Macromolecules 47
  * [30] H. Zimm and W.H. Stockmayer (1949) J. Chem. Phys. 17 1301
  * [31] W. Burchard and M. Schmidt (1980) Polymer 21 745
  * [32] B. Dünweg, D. Reith, M. Steinhauser and K. Kremer (2002) J. Chem. Phys. 117 914
  * [33] M. Fukatsu, M.J. and Kurata (1966) J. Chem. Phys. 44 4539
  * [34] E. Uehara and T. Deguchi (2016) J. Chem. Phys. 145 164905
  * [35] Y. A. Budkov and A. L. Kolesnikov (2021) J. Stat. Mech. 063213
  * [36] Y. Nakamuta et al. (2000) Macromolecules 33 8323
  * [37] H. MA and J. Xu (1994) Polymer J. 26 779
  * [38] Z. Zhiping, X. Jianmin, S. Xubing, and Y. Deyue (1992) Eur. Polym. J. 28 1339
  * [39] A. Chremos and J. F. Douglas (2018) J. Chem. Phys. 149 044904
  * [40] J. Paturej, S. Sheiko, S. Panyukov, and M. Rubinstein (2016) Sci. Adv. 2 e1601478
  * [41] T. Pan, S. Dutta and C. E. Sing (2022) J. Chem. Phys. 156 014903
  * [42] P. E. Theodorakis, N. G. Fytas (2012) American Journal of Condensed Matter Physics 2 101-108
  * [43] E. Alipour and J. F. Marko (2012) Nucleic Acids Res. 40 11202
  * [44] G. Fudenberg et al. (2016) Cell Rep. 15 2038\.
  * [45] A. Goloborodko, M.V. Imakaev, J.F. Marko, and L. Mirny (2016) eLife 5 e14864.
  * [46] L.A. Mirny (2021) Nature 590 554
  * [47] W. Fiers and R.L. Sinsheimer (1962) J. Mol. Biol. 5 424\.
  * [48] H.-X. Zhou (2003) J. Am. Chem. Soc. 125 9280\.
  * [49] des Cloizeaux J and Jannink G (1990) Polymers in Solutions: Their Modelling and Structure (Oxford: Clarendon Press)
  * [50] K. Haydukivska and V. Blavatska (2014) J. Chem. Phys. 141 094906
  * [51] I. Teraoka (2002) Polymer Solutions: An Introduction to Physical Properties (John Wiley & Sons Inc, New York)
  * [52] W. Radke and A. H. E. M/”uller (1996) Macromol. Theory and Simul. 5 759
  * [53] G. Wei (1995) Physica A 222 152; (1995) Physica A 222 155
  * [54] C. von Ferber, M. Bishop, T. Forzaglia, C. Reid, and G. Zajac (2015) J. Chem. Phys. 142 024901
  * [55] W.-C. W. Li, (2001) C. R. Acad. Sci. Paris 333 907