# Variation after projection calculations for high-spin states

Zao-Chun Gao<EMAIL_ADDRESS>China Institute of Atomic Energy, P.O. Box 275
(10), Beijing 102413, PR China

###### Abstract

In this paper, an improvement has been made on the variation after projection
(VAP) method, which is crucial in the calculations of high-spin states. It
turns out that, the form of the trial VAP wave function with spin $J$ can be
simplified by adopting just one projected state rather than previously
adopting all $(2J+1)$ angular momentum projected states for each selected
reference state, $|\Phi\rangle$. The present calculations show that such
simplification still minimizes the angular momentum projected energy up to a
very good approximation. In this simplified VAP scheme, one can obtain almost
equivalent VAP wave functions starting from different sets of the projected
basis states. This clearly shows that a nuclear state can not be identified
with a single intrinsic state, while in the traditional nuclear collective
models, an intrinsic state is usually assigned to a well deformed rotational
band.

###### keywords:

Variation after projection, high-spin states, shell model

††journal: Physics Letter B

Projection is a powerful technique that has long been used in various fields
of physics, such as particle physics [1], atomic physics [2] and the physics
of condensed matter [3]. In nuclear physics, the projection method has also
been used by many authors to improve the quality of the nuclear wave functions
[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. This technique
removes the spurious part of the trial wave function and makes the projected
wave function rather suitable in describing various properties of nuclei and
other quantum systems. The projected wave functions with different forms have
been varied by several authors [5, 8, 18, 13, 14, 16], so that they can be as
close as possible to the corresponding eigenfunctions obtained by the shell
model(SM) [20]. Such variation of the projected wave function is generally
called as variation after projection (VAP) [21]. The VAP methods are expected
to be applicable in large model space where full shell model calculation can
not be performed presently or in the near future.

If the VAP includes the angular momentum projection, it is believed that all
the angular momentum projected states should be used to construct the VAP wave
function. That means if one uses $n$ reference states to construct a VAP wave
function with spin $J$, then the number of included projected states is
$n(2J+1)$. In the calculations of high-spin states, this number may be large.
Clearly, too many projected states may increase the complexity of the VAP
calculation. Actually, large number of projected states may increase the
possibility of redundant projected states and may damage the stability of the
VAP iteration.

In this Letter, an improvement is made that only $n$ projected states are
adopted to construct the VAP wave function with arbitrary spin $J$. This makes
the VAP calculations for high-spin states as easy as those for $J=0$ states.
For simplicity, the reference states used in the present VAP are Slater
determinants (SD), so that the particle number projections can be omitted.

Let me start with the simplest variation of a SD, $|\Phi\rangle$, i.e., the
Hartree-Fock (HF) mean-field approximation. One can vary $|\Phi\rangle$ so
that the expectation energy, ${\langle\Phi|\hat{H}|\Phi\rangle}$, reaches a
minimum. However, the minimized HF energy and the corresponding wave function
$|\Phi\rangle$ are usually far away from the exact SM ones [10].

In the general case, the HF mean-field approximation breaks the symmetries of
rotation and reflection. One may project $|\Phi\rangle$ to get a set of
projected states, $P^{J}_{MK}P^{\pi}|\Phi\rangle$, where, $P^{J}_{MK}$ and
$P^{\pi}$ are the angular momentum and parity projection operators,
respectively. The state $P^{J}_{MK}P^{\pi}|\Phi\rangle$ has good quantum
numbers of spin ($J$), the magnetic quantum number in the laboratory frame
($M$) and parity ($\pi$). Notice that, $K$ is not a good quantum number, as
will be discussed later.

In my formalism, the SD $|\Phi\rangle$ can be assumed to be fully symmetry-
unrestricted, this makes it convenient to calculate all kinds of states in
even-even, odd-odd and odd-mass nuclei in the same way. Additionally, the
converged VAP energy can be lower than the one with symmetry-restricted SD.
But for the sake of simplicity, the presently used model spaces include only
one single harmonic oscillator (HO) shell. The SDs in such model spaces have
automatically a good parity. Therefore, I forget about the parity projection
and only the angular momentum projection is considered in the present work.

With given $J$ and $M$, one can project $|\Phi\rangle$ onto $2J+1$ different
projected states, $P^{J}_{MK}|\Phi\rangle$, with $K=-J$, $-J+1$, $\cdots$,
$J$. The nuclear wave function can be expressed in terms of these projected
states

$\displaystyle|\Psi_{JM\alpha}\rangle=\sum_{K=-J}^{J}f^{J\alpha}_{K}P^{J}_{MK}|\Phi\rangle,$
(1)

where $\alpha$ is used to label the states with the same $J$ and $M$. The
coefficients $f^{J\alpha}_{K}$ and the corresponding energy $E^{J}_{\alpha}$
of Eq. (1) can be determined by solving the Hill-Wheeler (HW) equation of
order $2J+1$,

$\displaystyle\sum_{K^{\prime}=-J}^{J}(H^{J}_{KK^{\prime}}-E^{J}_{\alpha}N^{J}_{KK^{\prime}})f^{J\alpha}_{K^{\prime}}=0,$
(2)

where $H^{J}_{KK^{\prime}}=\langle\Phi|\hat{H}P^{J}_{KK^{\prime}}|\Phi\rangle$
and $N^{J}_{KK^{\prime}}=\langle\Phi|P^{J}_{KK^{\prime}}|\Phi\rangle$. For
convenience, I assume $E^{J}_{1}\leq E^{J}_{2}\leq\cdots\leq E^{J}_{2J+1}$.
The coefficients $f^{J\alpha}_{K^{\prime}}$ should satisfy the normalization
condition.

$\displaystyle\sum_{K,K^{\prime}=-J}^{J}f^{J\alpha*}_{K}N^{J}_{KK^{\prime}}f^{J\alpha}_{K^{\prime}}=1.$
(3)

In Eq. (2), $f^{J\alpha}_{K}$ is independent of $M$. Since $M$ does not carry
any information of physics, it will no longer be discussed and simply be
regarded as a fixed quantum number.

To solve Eq. (2), the first step is the diagonalization of $N^{J}$ and one has

$\displaystyle\sum_{K^{\prime}=-J}^{J}N^{J}_{KK^{\prime}}R^{k}_{K^{\prime}}=n_{k}R^{k}_{K},$
(4)

where $n_{k}\geq 0$ and $R^{k}$ with $k=1,2,...2J+1$ are eigenvalues and the
corresponding eigenvectors, respectively. Here, I assume $n_{1}\geq
n_{2}\geq\cdots\geq n_{2J+1}$. Then one can establish a new set of orthonormal
basis states, $|\psi_{k}^{J}\rangle\,(k=1,2,\cdots,2J+1)$,

$\displaystyle|\psi_{k}^{J}\rangle=\frac{1}{\sqrt{n_{k}}}\sum_{K=-J}^{J}R^{k}_{K}P^{J}_{MK}|\Phi\rangle,$
(5)

and the HW equation (2) can be transformed into a normal eigenvalue equation

$\displaystyle\sum_{k^{\prime}=1}^{2J+1}\left[\langle\psi_{k}^{J}|\hat{H}|\psi_{k^{\prime}}^{J}\rangle-E^{J}_{\alpha}\delta_{kk^{\prime}}\right]u^{J\alpha}_{k^{\prime}}=0.$
(6)

Actually, the energies in Eq. (2) are obtained by solving Eq. (6) and the
coefficients of the wave functions, $f^{J\alpha}_{K}$, are obtained from
$u^{J\alpha}_{k}$, i.e.,

$\displaystyle
f^{J\alpha}_{K}=\sum_{k=1}^{2J+1}\frac{R^{k}_{K}u^{J\alpha}_{k}}{\sqrt{n_{k}}}.$
(7)

The problem is, the $2J+1$ projected states are not orthogonal to one another.
Hence in practical calculations, it is possible that some $n_{k}$ values can
be very tiny or even zero. Consequently, the corresponding
$|\psi_{k}^{J}\rangle$ basis states may not be precise enough to guarantee the
stability of the calculated energies and wave functions. To avoid this
trouble, the $|\psi_{k}^{J}\rangle$ states with tiny $n_{k}$ values should be
abandoned. This can be done by setting a cutoff parameter $\epsilon>0$ and
only those $|\psi_{k}^{J}\rangle$ states with $n_{k}>\epsilon$ are taken to
form the nuclear wave function (see, e.g., Ref. [22]).

Such basis cutoff works if $|\Phi\rangle$ remains unchanged. However, if one
tries to vary $|\Phi\rangle$ in attempt to get the best nuclear wave function
$|\Psi_{JM\alpha}\rangle$, in other words, to perform the VAP calculation, new
trouble arises due to the changes of $n_{k}$ values during VAP iteration. When
a fixed cutoff parameter $\epsilon$ is used, the number of selected
$|\psi_{k}^{J}\rangle$ states likely changes as VAP iteration goes on. This
will destroy the smoothness of the energy as a function of the variational
parameters. On the other hand, if one fixes the number of selected
$|\psi_{k}^{J}\rangle$ states, then some of them might have too small $n_{k}$
values at certain VAP iteration. Certainly, both treatments of the
$|\psi_{k}^{J}\rangle$ selection do not guarantee the stability of the VAP
iteration. These problems are exemplified in the figures included in the
supplemental material.

One can imagine that the safest VAP calculation might be the one that only
$|\psi_{1}^{J}\rangle$ with the largest $n_{1}$ is taken. But
$|\psi_{1}^{J}\rangle$ still includes several projected states
$P^{J}_{MK}|\Phi\rangle$. Sometimes, some of them might have very tiny norms,
$\langle\Phi|P^{J}_{KK}|\Phi\rangle$, or two projected states might be almost
the same. For instance,
$P^{J}_{MK}|\Phi\rangle=(-1)^{J-K}P^{J}_{M,-K}|\Phi\rangle$ when
$|\Phi\rangle$ is time-even [24]. These problems may still affect the
precision of the calculated nuclear wave functions. Finally, it looks the only
way of ensuring the VAP stability is that only one projected state
$P^{J}_{MK}|\Phi\rangle$ is taken. Then, the nuclear wave function can be
simplified as

$\displaystyle|\Psi_{JM}(K)\rangle=\frac{P^{J}_{MK}|\Phi\rangle}{\sqrt{\langle\Phi|P^{J}_{KK}|\Phi\rangle}},$
(8)

where, $K$ can be randomly chosen but should satisfy
$\langle\Phi|P^{J}_{KK}|\Phi\rangle\neq 0$. $K$ is put into the bracket since
it is not a good quantum number.

Now, there have two forms of nuclear wave functions, Eq. (1) and Eq. (8), to
be optimized. One may expect Eq. (1) is better than Eq. (8) if they have the
same $|\Phi\rangle$. Actually, in almost all previous VAP works [8, 18, 14,
16], the former Eq. (1) was taken. However, if one do the VAP calculations
with these two different forms of nuclear wave functions, independently, the
results are quite interesting.

Figure 1:  (Color online) Calculated VAP energies, $E^{J}_{1}$ (red lines) and
$E^{J}(K)$(filled dots), for the yrast $2^{+}$ and $7^{+}$ states in 24Mg,
with wave functions Eq.(1) and Eq.(8), respectively. The corresponding exact
shell model energies are shown as black lines. The USDB interaction is
adopted.

As a first example of such calculations, the USDB interaction [23] defined in
the $sd$ shell is used. The lowest $J^{\pi}=2^{+}$ and $7^{+}$ energies in
24Mg are calculated. The energy of Eq. (8) is written as

$\displaystyle E^{J}(K)=\langle\Psi_{JM}(K)|\hat{H}|\Psi_{JM}(K)\rangle.$ (9)

For each $K$, $E^{J}(K)$ is minimized independently and one gets $2J+1$
converged $E^{J}(K)$ energies at spin $J$. The VAP method of Ref. [14] is
applied to minimize $E^{J}_{1}$ in Eq.(2) and $E^{J}(K)$ in Eq(9),
respectively. The results are shown in Fig. 1 (see more details in the
supplemental material). It is seen that the $E^{J}_{1}$ energies are very
close to the exact ones calculated by the full SM which has been discussed in
Ref. [14]. Surprisingly, one can also see that all the $E^{J}(K)$ energies are
almost the same as $E^{J}_{1}$ energies for both $J^{\pi}=2^{+}$ and $7^{+}$
without exception. Correspondingly, the $2J+1$ converged wave functions
$|\Psi_{JM}(K)\rangle$ for $E^{J}(K)$ are obtained. If $K$ is a good quantum
number, these $|\Psi_{JM}(K)\rangle$ states (with common fixed $J$ and $M$)
should be orthogonal to each other. This is clearly not possible because it is
known the norm matrix, $N^{J}$, in Eq. (2) is usually not diagonal. On the
contrary, it is expected they are the same state since all $E^{J}(K)$ energies
converge to the same level, which wave function should be unique, thus one
should have

$\displaystyle|\langle\Psi_{JM}(K)|\Psi_{JM}(K^{\prime})\rangle|\approx 1$
(10)

for arbitrary $K$ and $K^{\prime}$. All possible overlaps,
$|\langle\Psi_{JM}(K)|\Psi_{JM}(K^{\prime})\rangle|$, for both $J^{\pi}=2^{+}$
and $7^{+}$ states have been calculated. Most of them are above 0.99 and the
worst overlap is still over 0.98. This clearly tells us that one may
arbitrarily select a $K$ number and use Eq. (8) to do the VAP calculation
without losing good approximation.

Before performing further calculations, it is necessary to understand the
reason why $E^{J}(K)$ energies are almost the same as $E^{J}_{1}$. At a given
spin $J$, one can assume that the nucleus is rotating. Once one has the
optimized SD $|\Phi\rangle$ for Eq. (1), $E^{J}_{1}$ is then obtained from the
Eq. (2). This also determines the orientation of rotational axis relative to
$|\Phi\rangle$ for $E^{J}_{1}$. Since $E^{J}_{1}$ is the lowest, the moment of
inertia along this rotational axis should be the largest. On the other hand,
one can first fix the orientation of the rotational axis relative to
$|\Phi\rangle$ and then do the energy minimization, which is actually done by
taking Eq.(8). It is expected that the moment of inertia along this fixed
rotational axis can also reach the same maximum as that for $E^{J}_{1}$, and
$E^{J}(K)$ converges to $E^{J}_{1}$, simultaneously.

Thus it is shown that the VAP calculations with wave functions in Eq. (1) and
Eq. (8) are almost equivalent, which seems to be a universal phenomenon. If
this is true, Eq. (8) can be used to simplify the VAP calculation and make the
VAP iteration equally stable at any spin. This is crucial in extending the
present VAP calculations to arbitrary high-spin states.

If one takes Eq. (8) as the nuclear wave function, then it can be fully
obtained from the $K$-projected state, $P^{K}|\Phi\rangle$. Here
$P^{K}=\frac{1}{2\pi}\int_{0}^{2\pi}e^{i(K-\hat{J}_{z})\phi}d\phi$ is the
$K$-projection operator. When two different reference states, $|\Phi\rangle$
and $|\Phi^{\prime}\rangle$, are taken, one can do the $K$-projection on them,
and get $P^{K}|\Phi\rangle$ and $P^{K^{\prime}}|\Phi^{\prime}\rangle$ with
$K\neq K^{\prime}$. These $K$-projected states are strictly orthogonal to one
another, but their total angular momentum ($J$) projected states can be almost
identical according to Eq. (10).

So far, there is only one SD in Eq. (8). To improve the VAP approximation, Eq.
(8) should be generalized by including more SDs, and the new form of VAP wave
function can be written as,

$\displaystyle|\Psi_{JM\alpha}(K)\rangle=\sum_{i=1}^{n}f^{J\alpha}_{i}P^{J}_{MK}|\Phi_{i}\rangle,$
(11)

where, $n$ is the number of included $|\Phi_{i}\rangle$ SDs. $f^{J\alpha}_{i}$
is determined by the following Hill-Wheeler equation,

$\displaystyle\sum_{i^{\prime}=1}^{n}(H^{J}_{ii^{\prime}}-E^{J}_{\alpha}N^{J}_{ii^{\prime}})f^{J\alpha}_{i^{\prime}}=0,$
(12)

where
$H^{J}_{ii^{\prime}}=\langle\Phi_{i}|\hat{H}P^{J}_{KK}|\Phi_{i^{\prime}}\rangle$
and $N^{J}_{ii^{\prime}}=\langle\Phi_{i}|P^{J}_{KK}|\Phi_{i^{\prime}}\rangle$.

One may expect Eq. (11) can also be used to describe the non-yrast states. To
check the validity of this assumption, here a newly developed VAP algorithm
[16] is applied, in which the low-lying state wave functions with the same
spin and parity can be varied on the same footing. This can be safely realized
by minimizing the sum of the corresponding low-lying energies,
$\sum_{\alpha=1}^{m}E_{\alpha}^{J}$. Here $m$ is the number of calculated
states.

It should be reminded there are two more potential problems with the projected
basis in Eq. (11), which may seriously affect the stability of the VAP
calculation. The one is the norm, $N^{J}_{ii}$, could become very tiny as the
VAP iteration proceeds. The other is the possibility of large overlaps among
the projected basis states. Both problems may damage the precision of the
calculated energies and the corresponding wave functions, and may cause the
collapse of the VAP iteration. Fortunately, such troubles do not appear in the
calculations of Fig. 1. However, to ensure the stability of VAP iteration, in
the following calculations, two constraint terms are attached to the energy
sum, and the final minimized quantity is,

$\displaystyle
Q=\sum_{\alpha=1}^{m}E^{J}_{\alpha}+\chi_{1}\sum_{i=1}^{n}\frac{1}{N^{J}_{ii}}+\frac{\chi_{2}}{2}\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{n}\frac{N^{J}_{ij}N^{J}_{ji}}{N^{J}_{ii}N^{J}_{jj}},$
(13)

where the second term tends to push the norms, $N^{J}_{ii}$, to large values,
and the third term tends to guide the projected basis states to be orthogonal
to one another. The values of the last two terms should be as small as
possible provided that the VAP iteration is stable. So the parameters
$\chi_{1}\geq 0$ and $\chi_{2}\geq 0$ should not be too large. In all the
following calculations, I use $\chi_{1}=10^{-4}$MeV and $\chi_{2}=1$MeV. This
makes the included projected states almost orthogonal to one another after the
VAP calculation converges. It turns out that both quantities of the last two
terms in Eq. (13) are within 100keV in the present calculations (see more
detailed calculations in the supplemental material).

Figure 2:  (Color online) Calculated lowest five $J^{\pi}=8^{+}$ eneriges in
some $sd$ shell nuclei with the Shell model and the VAP, respectively. The
five panels from low to high show the energies from the yrast $8^{+}_{1}$ ones
to the forth excited $8^{+}_{5}$ ones, respectively. VAP:$K=0$ and VAP:$K=8$
refer to the results of the VAP calculations with
$|\Psi_{JM\alpha}(K=0)\rangle$ and $|\Psi_{JM\alpha}(K=8)\rangle$,
respectively. For each nucleus, all calculated energies are shifted by the
same quantity so that the yrast SM enery is zero (see the lowest panel). The
USDB interaction is adopted.

Now, let us come to the calculations for the non-yrast states. In the second
example, the USDB interaction is adopted again and the results are shown in
Fig. 2. Here, the wave functions $|\Psi_{JM\alpha}(K)\rangle$ [see Eq. (11)]
with $n=10$ are taken to be varied. The quantitiy of Eq. (13) including the
sum of the lowest five energies ($m=5$) is minimized. To show the generality
of Eq. (11), the calculations include even-even, odd-odd and odd-mass nuclei.
For even-even and odd-odd nuclei, the calculated spin is taken to be $J=8$,
and for the odd-mass ones, $J=17/2$. As we learned from the first example, the
$K$ number in $|\Psi_{JM\alpha}(K)\rangle$ may be arbitrarily chosen within
the range $|K|\leq J$. For the $J=8$ case, $K=0$ and $K=8$ are chosen, and the
corresponding VAP energies are denoted by VAP:$K=0$ (red dots) and VAP:$K=8$
(blue circles), respectively in Fig. 2a. Similar calculations are performed
for the $J=17/2$ states but $K=1/2$ and $K=17/2$ are chosen and the
corresponding energies are shown in Fig. 2b. Let us first look at the Fig. 2a.
All the calculated energies of VAP:$K=8$ perfectly coincide with those of
VAP:$K=0$. This implies that the results of VAP:$K=0$ and VAP:$K=8$ are the
same. One can imagine VAP results with other $K$ should also be the same as
VAP:$K=0$ (or VAP:$K=8$). Comparing with the full shell model, the energies of
both VAP:$K=0$ and VAP:$K=8$ are very close to the same exact ones. Thus it is
expected that both $|\Psi_{JM\alpha}(K=0)\rangle$ and
$|\Psi_{JM\alpha}(K=8)\rangle$ are close to the same shell model wave function
for each calculated state. This means $|\Psi_{JM\alpha}(K=0)\rangle$ and
$|\Psi_{JM\alpha}(K=8)\rangle$ are almost the same. Indeed, for all calculated
states in Fig. 2a, it is found that

$\displaystyle|\langle\Psi_{JM\alpha}(K=0)|\Psi_{JM\alpha}(K=8)\rangle|\geq
0.95.$ (14)

Actually, most of such overlaps are above $0.98$. For the $J=17/2$ states, the
equivalence of $|\Psi_{JM\alpha}(K=1/2)\rangle$ and
$|\Psi_{JM\alpha}(K=17/2)\rangle$ are also numerically confirmed.

One can understand that, for an arbitrary selected $K$ with $|K|\leq J$, if
the number of selected $|\Phi_{i}\rangle$ SDs, $n$, is large enough, then the
space spanned by the $\hat{P}^{J}_{MK}|\Phi_{i}\rangle$ states may fully cover
the whole $J$-scheme shell model configuration space without including other
projected basis states, $\hat{P}^{J}_{MK^{\prime}}|\Phi_{i}\rangle$ with
$K^{\prime}\neq K$. Together with the good approximation of the present
calculations, it seems that there is no need to consider the complicate
$K$-mixing in all VAP calculations.

Figure 3:  (Color online) Calculated results of the yrast states in 48Cr with
the VAP and the shell model. (a) The calculated energies. (b) The B(E2) values
from the wave functions corresponding to (a). (c) The values of the overlap,
$|\langle\Psi_{JM1}(K=0)|\Psi_{JM1}(K=J)\rangle|$. The KB3 interaction is
adopted.

As a more practical application, the high-spin states in 48Cr are calculated,
which has been studied by the shell model [25], the projected shell model
[26], and the projected configuration interaction(PCI) method [27]. Here the
KB3 interaction [28] is taken, as has been used in Ref.s [25, 27]. This time,
the wave functions $|\Psi_{JM\alpha}(K=0)\rangle$ and
$|\Psi_{JM\alpha}(K=J)\rangle$ are used to minimize the yrast
energies($\alpha=1$) in 48Cr, respectively. In practical calculations, the $Q$
quantity of Eq. (13) with m=1 and n=5 is minimized. The calculated results
have been shown in Fig. 3. From Fig. 3(a), the VAP energies with both $K=0$
and $K=J$ are very close to the SM ones. The B(E2) values calculated with the
wave functions corresponding to Fig. 3(a) are shown in Fig. 3(b). All the
B(E2) values are bunched up tightly. The B(E2) values with wave functions
$|\Psi_{JM1}(K=0)\rangle$ and $|\Psi_{JM1}(K=J)\rangle$ are almost coincide
with each other. This again implies $|\Psi_{JM1}(K=0)\rangle$ and
$|\Psi_{JM1}(K=J)\rangle$ are almost the same. The overlap,
$|\langle\Psi_{JM1}(K=0)|\Psi_{JM1}(K=J)\rangle|$, as a function of spin $J$
has been shown in Fig. 3(c). Indeed, the overlaps for all spins are above
$98\%$.

Therefore, all the present calculations have clearly shown that the VAP wave
function can be expressed in terms of the projected states,
$P^{J}_{MK}|\Phi_{i}\rangle$, with the same $K$ values. Without the need of
performing the $K$-mixing, the present VAP method can be easily extended to
high-spin states. If one prefers to use HFB vacua rather than taking the SDs
in VAP, it is likely such extension may still works. However, the reference
states can not be uniquely determined in the VAP calculation at a given spin.
So, the intrinsic state can not be identified for a single nuclear state. In
the traditional nuclear collective models, an intrinsic state is usually
assigned to a well deformed rotational band. Nevertheless, if one wants to use
the simple, but limited, picture where a single intrinsic state is associated
with the ground-state wave function of a rotational band within a VAP scheme,
a possible way is to take a common reference state and generate all angular
momentum ($J$) projected states with the same $K$ values from it, then
minimize the sum of the projected energies of the band members. In that case,
the intrinsic state and the $K$ value might be determined. The details of such
work will be published in the future.

Acknowledgements This work is supported by the National Natural Science
Foundation of China under Grant Nos. 11975314,11575290, and by the Continuous
Basic Scientific Research Project Nos. WDJC-2019-13, BJ20002501.

## References

  * [1] E. G. Lubeck, M. C. Birse, E. M. Henley, and L. Wilets, Phys. Rev. D 33 (1986) 234.
  * [2] M. Bylicki, Phys. Rev. A 45 (1992) 2079.
  * [3] Q. B. Yang and W. D. Wei, Phys. Rev. Lett. 58 (1987) 1020.
  * [4] K. Hara, Y. Sun, Int. J. Mod. Phys. E 4 (1995) 637.
  * [5] Y.Kanada-En’yo, Phys. Rev. Lett. 81 (1998) 5291.
  * [6] T. Otsuka, M. Honma, T. Mizusaki, N. Shimizu, and Y. Utsuno, Prog. Part. Nucl. Phys. 47 (2001)319.
  * [7] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75 (2003) 121.
  * [8] K. W. Schmid, Prog. Part. Nucl. Phys. 52 (2004) 565.
  * [9] M. Bender and P.-H. Heenen, Phys. Rev. C 78 (2008)024309.
  * [10] Z.-C. Gao, and M. Horoi, Phys. Rev. C 79 (2009) 014311.
  * [11] J. M. Yao, J. Meng, P. Ring, and D. Vretenar, Phys. Rev. C 81 (2010) 044311.
  * [12] B. Bally, B. Avez, M. Bender, and P.-H. Heenen, Phys. Rev. Lett. 113(2014)162501.
  * [13] Z.-C. Gao, M. Horoi, Y. S. Chen, Phys. Rev. C 92 (2015) 064310.
  * [14] T. Ya, Y. He, Z.-C. Gao, J.-Q. Wang, and Y. S. Chen, Phys. Rev. C 95 (2017) 064307.
  * [15] L. M. Robledo, T. R. Rodríguez and R. R. Rodríguez-Guzmán, J. Phys. G: Nucl. Part. Phys. 46 (2019) 013001.
  * [16] J.-Q. Wang, Z.-C. Gao, Y.-J. Ma, and Y. S. Chen, Phys. Rev. C 98 (2018) 021301(R).
  * [17] J. M. Yao, B. Bally, J. Engel, R. Wirth, T. R. Rodríguez and H. Hergert, Phys. Rev. Lett. 124(2020)232501.
  * [18] N. Shimizu, Y. Tsunoda, Y. Utsuno, and T. Otsuka, Phys. Rev. C 103 (2021)014312.
  * [19] B. Bally and M. Bender, Phys. Rev. C 103 (2021) 024315.
  * [20] E. Caurier, G. Martínez-Pinedo, F. Nowacki, A. Poves, A. P. Zuker, Rev. Mod. Phys. 77 (2005) 427.
  * [21] P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer Verlag, New York/Heidelberg/Berlin, 1980).
  * [22] H. Zduńczuk, W. Satuła, J. Dobaczewski, and M. Kosmulski, Phys. Rev. C 76 (2007) 044304.
  * [23] B. A. Brown and W. A. Richter, Phys. Rev. C 74 (2006) 034315.
  * [24] Z.-C. Gao, Y. S. Chen and J. Meng, Chin. Phys. Lett. 19 (2002) 650
  * [25] E. Caurier, J. L. Egido, G. Martínez-Pinedo, A. Poves, J. Retamosa, L. M. Robledo, and A. P. Zuker, Phys. Rev. Lett. 75 (1995)2466.
  * [26] K. Hara, Y. Sun, and T. Mizusaki, Phys. Rev. Lett. 83 (1999) 1922.
  * [27] Z.-C. Gao, M. Horoi, Y. S. Chen, Y. J. Chen and Tuya, Phys. Rev. C 83 (2011) 057303.
  * [28] A. Poves and A. P. Zuker, Phys. Rep. 71 (1981) 141.


# Topological and non-topological mechanisms of loops formation in
chromosomes: effects on the contact probability

Kirill Polovnikov1,2, Bogdan Slavov2 1 Institut Curie, PSL Research
University, Sorbonne Université, CNRS UMR3664, Paris, France
2 Skolkovo Institute of Science and Technology, 121205 Moscow, Russia

###### Abstract

Chromosomes are crumpled polymer chains further folded into a sequence of
stochastic loops via loop extrusion. While extrusion has been verified
experimentally, the particular means by which the extruding complexes bind DNA
polymer remains controversial. Here we analyze the behaviour of the contact
probability function for a crumpled polymer with loops for the two possible
modes of cohesin binding, topological and non-topological mechanisms. As we
show, in the non-topological model the chain with loops resembles a comb-like
polymer that can be solved analytically using the quenched disorder approach.
In contrast, in the topological binding case the loop constraints are
statistically coupled due to long-range correlations present in a non-ideal
chain, which can be described by the perturbation theory in the limit of small
loop densities. As we show, the quantitative effect of loops on a crumpled
chain in the case of topological binding should be stronger, which is
translated into a larger amplitude of the log-derivative of the contact
probability. Our results highlight a physically different organization of a
crumpled chain with loops by the two mechanisms of loops formation.

##  I Introduction

Chromosomal DNA polymers are organized into loops which serve both structural
and gene regulatory roles inside the cell nucleus [1, 2]. Folding into
chromatin loops provides an effective means of compaction of 10cm-long human
chromosomes inside a micron-sized nucleus. In particular, it was proposed that
upon the entry to mitosis a fluffy interphase chromosome is arranged into a
dense array of loops via the mechanism of lengthwise compaction [3], which
results in the bottlebrush structure long observed in mitosis [4, 5].

Interphase chromosomes are also organized into loops, however, of
significantly lower linear density compared to the mitotic chains [1, 6, 7].
The loops are formed by SMC proteins (e.g, cohesin, condensin), and their
extruding ability has been independently demonstrated in single-molecule
imaging by two groups [8, 9]. In human cells cohesin provides the main source
of short-scales loops on chromosomes. However, the mechanism of cohesin
binding to chromosomes is not known. It has been proposed that cohesin can
topologically embrace two strands of DNA, forming a loop of progressively
increasing size in time [10, 2]. In this topological binding model cohesin
complex resembles a physical ring that can entrap two DNA strands without
chemically binding them. Alternatively, in the non-topological model, cohesin
operates as a cross-linker: in the base of the loop it creates a bond between
two sites of the chain which is moving in the course of extrusion [10, 2].
Which of the models is realized in the cell is still under debate. As we show
here, such a subtle difference in the microscopic structure of cohesin complex
can affect chromosome organization at all scales.

A crucial difficulty towards the analytical description of a chromosome as a
chain with loops is non-ideal statistics of chromosomal polymer. The physical
state of a chromosome can be characterized by the scaling of of the average
contact probability function $P(s)$ with the contour distance $s=|i-j|$ (an
analogue of the return probability of a random walk to the origin). For
chromosomes this function can be extracted from Hi-C experiments [11]. As
various data suggest, at large scales, in the range of $s=1-5$ Mb, the scaling
of $P(s)\sim s^{-1}$, is different from $\sim s^{-3/2}$ as expected for a
three-dimensional random walk, indicating that chromosomes are folded into
non-ideal states [11, 12, 13]. This $\approx-1$ exponent of the contact
probability corresponds to the fractal dimension $d_{f}\approx 3$, which is an
asymptotic property of topologically crumpled chains. At shorter scales
$s<1$Mb the scaling of $P(s)$ has a characteristic shoulder, which can be
recapitulated by the simplest model of a fractal chain with cohesin-mediated
loops breaking the scale invariance [7]. Importantly, upon elimination of
cohesin loops in experiments [14, 15], the $\sim s^{-1}$ scaling stretches to
shorter distances up to $s\approx 40$kb, suggesting that short-scale cohesin
loops are formed on top of a crumpled polymer with strong (power-law)
correlations along the genomic contour [7].

Here we demonstrate that due to the intrinsic long-range correlations in
folding of crumpled chromosomes the two models of cohesin binding produce
drastically different relaxation behaviour of the chain with loops. In the
non-topological model a chromosome resembles a comb-like polymer with a
backbone and loopy side chains. Conformation properties of the chain at large
scales are controlled by the backbone (the shortest path) in this model. In
contrast, in the topological binding model the backbone is not formed and the
polymer remains linear, while folded into loops. Importantly, in this case the
chain length is not shortened by the addition of the loops. Furthermore, all
the loops are statistically coupled with each other in the topological model,
being the segments of a single crumpled polymer. Note that this physical
difference between the models vanishes in the case of ideal chain statistics,
$d_{f}=2$, due to absence of long-range correlations in the chain.

We investigate the effect of the specific binding mechanism on the resulting
contact probability function, $P(s)$, for the two models. We treat the loops
as frozen disorder on the fractal polymer [7] and assume formation of a well-
defined sequence of loops and gaps on the chain with exponential distribution
of their lengths (no nested or overlapping loops are allowed) fixed on each
conformation. In order to take into account intrinsic correlations in the
crumpled chains, we make use of the Gaussian measure of fractional Brownian
paths (fBm) as recently proposed [16]. We show that while the non-topological
model under these assumptions can be computed analytically [7], the
topological model is not fully analytically tractable. For the latter case it
is possible to construct a perturbation by a single loop (the one-loop
approximation), and analytically compute the first order correction to the
power-law scaling of $P(s)$. Despite of the lack of the full analytical
solution for the topological model, our results show how the two models
generate different statistics of contacts across scales.

##  II A fractal polymer folded into random loops

We exploit a conventional bead-spring model of a flexible polymer chain,
generalized for the arbitrary fractal dimension $d_{f}$, i.e. the average
square of the segment size of length $s$ is

$r^{2}_{0}(s)=b^{2}s^{2/d_{f}},$ (1)

where $b$ is a scale of a single bead. In the case of absence of any
interactions between the beads other than harmonic polymer bonds, the
statistics of chain is ideal with $d_{f}=2$ [17]. A swollen chain in a good
solvent in three dimensions has fractal dimension $d_{f}\approx 1.7$ [18],
however, excluded volume interactions are likely screened at the relevant
scales for chromatin [13]. Topological constraints (in 3D) in the fractal
globule model force the chain to fold with the asymptotic fractal dimension
$d_{f}=3$ [19, 20]. Without excluded volume $d_{f}$ can be larger than the
dimension of the space: e.g., a phantom randomly-branching polymer has
$d_{f}=4$ [21]. In general, in order to describe a polymer with any fractal
dimension $d_{f}\geq 2$ an effective quadratic Hamiltonian was proposed in
[16]. It generates Gaussian polymer conformations, which map onto trajectories
of a fractional Brownian particle (fBm).

Figure 1: Two possible binding mechanisms of cohesin to DNA. (A). Non-
topological model. (B). Topological model.

Next we consider an ensemble of fractal chains folded into sequences of loops
and gaps, see Fig. 1(b). The loops on each chain in the ensemble are randomly
positioned, have a random contour length drawn from exponential distribution
with the average $\lambda$ and are separated by exponentially distributed gaps
(spacers) with the average length $g$. The loops are fixed for each chain.
Similarly to our earlier work [7], we follow the frozen disorder approach,
i.e. each chain in the ensemble is thermally equillibrated together with its
set of the loops. Importantly, the fractal statistics (Eq.1) is considered as
the inherent property of the chain, owned to specific interactions in the
polymer (e.g. large-scale topological constraints), which are not perturbed by
addition of the loops.

In this paper we focus on the difference between the two models of the loop
formation, depicted in Fig. 1. In the non-topological binding model cohesin
complex forms a covalent bond in the loop base. Therefore, a chain with loops
can be partitioned into a fBm backbone (main chain) and a set of fBm bridges
(loops). Notably, the main chain and the loops are fractal (generally, non-
ideal) chains both having fractal dimension $d_{f}$ and statistically
independent from each other. Relaxation of the polymer at large scales takes
place along the shortest path, i.e along the backbone ($s_{eff}$, see Fig.
1A). In contrast, in the case when the ring-shaped binding protein can embrace
two strands of chromatin (topological model of binding, Fig. 1B), the
chromosome backbone is not formed and stress in the chain propagates along the
whole polymer contour. A spacer and its neighbouring loop belong to the same
fBm polymer, which implies strong (power-law in the contour distance) tangent
correlations between them for $d_{f}>2$. Thus, these correlations are not
destroyed by addition of the loops, which physically determines the difference
in the contact statistics of a chain with loops at small and large scales.

Figure 2: Non-topological model. Diagrams contributing to the contact
probability, $P(s)$: (a) spacer-spacer, (b) loop-spacer and (c) loop-loop
(same), (d) loop-loop (other).

###  II.1 Non-topological model

Decoupling of correlations between the backbone and the loops in the non-
topological model allows for the analytical expression for the $P(s)$, which
was derived recently by us in [7]. We briefly describe the main steps and
results below.

First, one needs to compute the equilibrium contributions of different
diagrams, which classify relative positions of the points $i,j$ ($s=|i-j|$)
with respect to the loop bases: (a) spacer-to-spacer, (b) loop-to-spacer, (c)
loop-loop (same), (d) loop-loop (different). For each diagram one computes the
variance of the spatial distance $\sigma^{2}(i,j)=\langle r^{2}(i,j)\rangle$
between the points $i,j$ of interest and makes use of the Gaussian relation
$P(i,j)\sim\sigma^{-3/2}(i,j)$ for the corresponding equilibrium contact
probability. For that the vector $\vec{r}$ is decomposed into the sum of
independent vectors (see Fig. 2) ${\vec{r}}_{1}$ (loop), ${\vec{r}}_{2}$
(backbone) for the diagram (b); ${\vec{r}}_{1}$ (loop), ${\vec{r}}_{2}$
(backbone), ${\vec{r}}_{3}$ (loop) for the diagram (d). In the diagrams
(b)-(d) a loop resembles a fBm bridge of the same dimension $d_{f}$. The
analytical tractability of this model is owned to the statistical independence
of the vectors ${\vec{r}}_{i},i=1,2,3$, as they correspond to separate
branches (side loops and the backbone) in the comb-like polymer. Second, one
averages the resulting contact probabilities $P(i,j)$ over all possible pairs
of monomers $i,j$ belonging to different diagrams, such that $|i-j|=s$, using
the appropriate weights of the diagrams. The exponential distribution of loops
and gaps sizes allows to make use of the well-known result for the propagators
of the two-state Markov process and properly weigh contributions of different
diagrams [22]. Finally, the remaining averaging over the distribution of
random loops and gaps is performed, which results in the sum of multiple
integrals involving Bessel functions [7].

Let us see how $P(s)$ for different diagrams look like in this model. As noted
above, the variance of the vector connecting any two points at distance $s$ in
the chain can be expressed through the variances for a segment free of loops,
$\sigma^{2}_{\text{free}}(s)$ (the backbone) and for a fBm bridge of size $L$,
$\sigma^{2}_{\text{bridge}}(s,L)$. For a free segment along the backbone the
result is given by (Eq.1), i.e.

$\sigma^{2}_{\text{free}}(s)=b^{2}s^{2/d_{f}},$ (2)

as the backbone constitutes a fBm trajectory, which is statistically
independent of the intervening of loops. The bridge of size $L$ with the
fractal dimension $d_{f}$ can be described as a fBm polymer with the same
dimension conditionally constrained by the loop. Making use of the Gaussian
measure of the fBm paths, one can derive the following result for the fBm
bridge (see [7] and derivation of the conditional propagators for the
topological model below)

$\displaystyle\sigma^{2}_{\text{bridge}}(s,L)=$ (3)
$\displaystyle=\sigma^{2}_{\text{free}}(s)\left(1-\frac{(\sigma^{2}_{\text{free}}(L)+\sigma^{2}_{\text{free}}(s)-\sigma^{2}_{\text{free}}(L-s))^{2}}{4\sigma^{2}_{\text{free}}(s)\sigma^{2}_{\text{free}}(L)}\right),$

Note, that for the ideal chain $d_{f}=2$, one obtains the well-known
expression for the ideal bridge

$\sigma^{2}_{\text{bridge}}(s,L)=\sigma^{2}_{\text{free}}(s)\frac{L-s}{L}=\frac{b^{2}s(L-s)}{s}.$
(4)

Now let us consider the diagram (a). In the presence of arbitrary number of
non-topological loops between two points of interest the contact probability
reads

$P_{(a)}(s,x)=\frac{1}{(2\pi\sigma_{\text{free}}^{2}[(1-x)s])^{3/2}},$ (5)

where $x$ ($0\leq x<1$) denotes the fraction of the subchain length occupied
by the loops. In other words, the intervening loops lead to reduction of the
effective contour distance between points of interest.

Next, let us consider a subchain of length $s$ with one end belonging to the
gap region and another end belonging to the loop. The loop containing one of
the two sites of interest is parametrized by the lengths $l_{1}$ and $l_{2}$
as shown in Fig 2b. Clearly, the separation vector $\vec{r}$ between two sites
can be represented as a sum of mutually independent zero mean Gaussian random
vectors, $\vec{r}_{1}$ and $\vec{r}_{2}$. Thus, we have

$\begin{split}&P_{(b)}(s,l_{1},l_{2},x)=\\\
&=\frac{1}{(2\pi(\sigma_{\text{bridge}}^{2}[l_{2},l_{1}+l_{2}]+\sigma_{\text{free}}^{2}[(1-x)(s-l_{2})])^{3/2}},\end{split}$
(6)

where $0\leq x<1$, $l_{1}\geq 0$ and $0\leq l_{2}\leq s$. As before, $x$ is
the fraction of contour length along the backbone occupied by the loops.

Now let us consider a subchain located inside the loop, see Fig. 2c. This
situation corresponds to a pure bridge, therefore

$P_{(c)}(s,l_{1},l_{2})==\frac{1}{(2\pi\sigma_{\text{bridge}}^{2}[s,l_{1}+l_{2}])^{3/2}},$
(7)

with the parameters belonging to the range $l_{1}\geq 0$ and $l_{2}\geq s$.

Finally, for the diagram (d) the vector $\vec{r}$ is decomposed into the sum
of independent Gaussian vectors $\vec{r}_{1}$ in the first bridge,
$\vec{r}_{2}$ along the backbone and $\vec{r}_{3}$ in the second bridge. Thus,
the contact probability for this diagram reads

$\displaystyle
P_{(d)}(s,l_{1},l_{2},h,\tilde{l}_{1},\tilde{l}_{2},x)=\frac{1}{(2\pi\sigma_{d}^{2})^{3/2}},$
(8)

with the following variance

$\displaystyle\sigma^{2}_{d}=\sigma_{\text{bridge}}^{2}[l_{2},l_{1}+l_{2}]+\sigma_{\text{free}}^{2}[(1-x)h]+$
(9)
$\displaystyle+\sigma_{\text{bridge}}^{2}[\tilde{l}_{2},\tilde{l}_{1}+\tilde{l}_{2}].$
(10)

Expressions (Eq.5),(Eq.6),(Eq.7) and (Eq.8) determine the contact
probabilities for all four classes of diagrams for any loop density
$\lambda/g$ and any fractal dimension $d_{f}$ of the chain. The final result
for the contact probability in the non-topological model is given by the
weighted sum

$\displaystyle P^{\text{non-topo}}(s)=\sum_{i=a,b,c,d}\langle
p_{(i)}(s\mid\\{A\\}_{i})\rangle,$ (11)

where $\langle p_{(i)}(s\mid\\{A\\}_{i})\rangle$ are the particular
contributions to the contact probability from each diagram $i$, integrated
over the possible values of diagram-specific parameters $\\{A\\}_{i}$ and
properly weighted. In particular, the exponential distribution of loops and
gaps sizes allows to make use the probabilities for the two-state Markov
process for calculation of the weights (see [7] for details).

###  II.2 Topological model

Figure 3: One-loop approximation for the topological model. Diagrams (a)
spacer-spacer, (b) loop-spacer and (c) loop-loop (same).

In contrast, in this case all loops and spacers between them are statistically
coupled for $d_{f}\neq 2$, which does not allow to derive the full analytical
expression accounting for all the four diagrams. Indeed, in a chain with long-
range correlations the conditional probability of contact between any points
$i,j$ depends on all the constraints in the chain. In the language of random
walks, due to the intrinsic long-range memory in a fBm walk with the Hurst
parameter $H>1/2$ ($H=1/d_{f}$), the conditional return probability at time
$t=T$ depends on all the intermediate returns (loops) at previous times $t<T$.

However, treating the loops as a perturbation to the fractal chain scaling,
one can obtain a self-consistent linear order correction for
$P^{\text{topo}}(s)$ using the Gibbs measure of fBm polymer conformations. For
that the non-Markovian propagators of the single loop diagrams (see Fig. 3)
should be computed.

Propagators of single-loop diagrams. We have to compute the looping
probability of a fBm random walk at time points $i,j>i$ conditioned on the
loop at time points $k,n$, such as $k<n$. There are three one-loop diagrams
corresponding to this situation that we aim to process (Fig. 3):

  * •

(a) ”spacer-spacer”:
(a0) $k<n<i<j$ the points $i,j$ both locate in the gap right next to the loop;
(a1) $i<k<n<j$ the points $i,j$ belong to the consecutive gaps and there is
one loop in between them;

  * •

(b) ”loop-spacer”: $k<i<n<j$ one point belongs to a loop while the other
belongs to the gap next to it;

  * •

(c) ”loop-loop”: $k<i<j<n$ both points are inside the loop.

Figure 4: Difference between the contact statistics of the two models for
different diagrams and $d_{f}=3$. On all panels $t$ is the distance between
the left point and the loop base. (A). For diagram (a) with one loop in
between the points the topological model yields a smaller contact probability
than the non-topological. (B). For diagram (a) without loops in between the
topological model has a slightly larger $P^{topo}$ than $P^{non-topo}$ for the
points located close to the loop base. (C). For diagram (b) the topological
model yields a larger contact probability, especially when the loop size is
larger than the contour distance between the points, $s<l$.

The constrained propagators $P(i,j\mid k,n)$ for the points $i,j$ subject to
the loop $k,n$ can be formally written as follows

$P(i,j\mid k,n)=\frac{\int
P(\bm{X})\delta(\bm{x}_{k}-\bm{x}_{n})\delta(\bm{x}_{i}-\bm{x}_{j})D\bm{X}}{\int
P(\bm{X})\delta(\bm{x}_{k}-\bm{x}_{n})D\bm{X}}$ (12)

with the Gibbs distribution $P(\bm{X})$ corresponding to the effective
Hamiltonian of the fBm polymer. Since all the integrals in (Eq.12) are
Gaussian, they can computed using the Green function (see our previous works
for details [23, 7]). The result of the computation reads:

$P(i,j\mid
k,n)=\frac{b^{3}(n-k)^{3H}}{(2\pi)^{3/2}3^{3/2}}\frac{1}{(\det\Sigma)^{3/2}},$
(13)

where the matrix $\Sigma$ is the covariance matrix built on the vectors
$x_{i}-x_{j}$ and $x_{k}-x_{n}$:

$\begin{cases}\vskip 6.0pt plus 2.0pt minus
2.0pt\sigma_{11}=b^{2}(n-k)^{2H}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\sigma_{22}=b^{2}(j-i)^{2H}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\sigma_{12}=\frac{b^{2}}{2}\left(|i-n|^{2H}+|j-k|^{2H}-|i-k|^{2H}-|j-n|^{2H}\right).\end{cases}$
(14)

The different diagrams depicted in Fig. 3 correspond to different values of
the determinant of the covariance matrix $\Sigma$ in (Eq.13). Introducing
dimensionless variables $x=1-l/s$ and $y=(k-i)/s$, where $s=j-i$ and $l=n-k$,
one can express the propagator as follows

$P(i,j\mid k,n)=P(s,x,y)=\frac{1}{s^{3H}}\times I(x,y)$ (15)

where the loops-induced correction function $I(x,y)$ is

$I(x,y)\propto\frac{1}{\left(1-\frac{1}{4}(1-x)^{-2H}g^{2}(x,y)\right)^{3/2}},$
(16)

and the function $g(x,y)$ has the following universal form

$g(x,y)=(1-y)^{2H}+|y+1-x|^{2H}-|x-y|^{2H}-|y|^{2H}.$ (17)

The absolute values in (Eq.17) are to be expanded depending on relative
positions of the monomers in each particular diagram. Clearly, the function
$g(x,y)$ reflects the additional damping factor to the unconstrained contact
probability $s^{-3H}$ in (Eq.15), which is due to the loop formation.
Importantly, it contains information on the non-Markovian properties of the
fractal polymer, specific to the topological binding model.

The weights of the diagrams in the one-loop approximation. The weights in the
one-loop approximation can be derived in a limit $\lambda/g\ll 1,s/g\ll 1$ of
the general expressions in the full theory, assuming exponential distributions
of loop and gap lengths. In the most general terms the weights of each diagram
can be calculated as a fraction of configurations of two points $i,j$ at the
fixed distance $s=|i-j|$ with respect to the coordinates of the monomers in
the loop base $k,n$ that fall into the corresponding diagrammatic class,
averaged over all possible loop lengths $l=|k-n|$. For the diagram (a1) with
one loop this gives the following result

$W_{a_{1}}=\frac{1}{\lambda+g}\int_{0}^{s}dl\rho(l)\int_{0}^{s-l}dt_{1}$ (18)

and $t_{1}=|k-i|=|y|s$ is the contour distance from the left point to the loop
base. The diagram (a0) with no loops in between the points $i,j$ we have two
symmetric locations of the points with respect to the loop, thus, it should be
taken with the factor of $2$

$W_{a_{0}}=\frac{2}{\lambda+g}\int_{0}^{\infty}dl\rho(l)\int_{0}^{\frac{g-s}{2}}dt_{1}=\frac{g-s}{\lambda+g}\approx
1-\frac{\lambda}{g}-\frac{s}{g}$ (19)

For the diagram (b) we distinguish two cases $l<s$ and $l>s$ and also have two
symmetric cases as for the previous diagram

$W_{b}=\frac{2}{\lambda+g}\left(\int_{0}^{s}dl\rho(l)\int_{0}^{l}dt_{1}+\int_{s}^{\infty}dl\rho(l)\int_{l-s}^{l}dt_{1}\right)$
(20)

And the diagram (c)

$W_{c}=\frac{1}{\lambda+g}\int_{s}^{\infty}dl\rho(l)\int_{0}^{l-s}dt_{1}$ (21)

It is easy to check that in the $\lambda/g\to 0$ all the above weights indeed
sum to $1$. For the exponential distribution of the loop lengths
$\rho(l)=\lambda^{-1}\exp(-l/\lambda)$ the weights can be explicitly
calculated

$\displaystyle
W_{a_{1}}=\frac{s-\lambda}{g}+\frac{\lambda}{g}\exp(-s/\lambda)$ (22)
$\displaystyle W_{a_{0}}=1-\frac{\lambda}{g}-\frac{s}{g}$ (23) $\displaystyle
W_{b}=\frac{2\lambda}{g}\left(1-\exp(-s/\lambda)\right)$ (24) $\displaystyle
W_{c}=\frac{\lambda}{g}\exp(-s/\lambda)$ (25)

Expanding the propagator for each single-loop diagram. All diagrams can be
processed in a similar way, so that the respective weights are effectively
encoded in the limits of integration over universal variables $x$ and $y$. For
example, for the spacer-spacer diagram (a1) with one loop between points $i,j$
the natural coordinates are the loop size $l$ and the distance from the left
point of interest to the left point in the loop base, $t_{1}=k-i$

$\displaystyle
P_{a_{1}}(s)=\frac{1}{s^{3H}}\frac{1}{\lambda+g}\int_{0}^{s}dl\rho(l)\int_{0}^{s-l}dt_{1}I(x(l,s),y(t_{1},s))$
(26)
$\displaystyle\approx\frac{1}{s^{3H}}\frac{s}{g}\int_{0}^{1}dx\;\rho(s(1-x))\int_{0}^{x}dy\;I(x,y)$
(27)

Note that the spacer-spacer diagram (a) without any loops between the points
$i,j$ (”free propagator”) is trivial only for the Markovian case $H=1/2$. For
$H<1/2$ formation of a spatial contact between $i,j$ depends on the fact of
loop formation between $k$ and $n$ nearby. Introducing an auxiliary variable
$t_{1}=i-n$ for the distance from the loop to the left point $i$ we can
express the propagator as

$\begin{split}&P_{a_{0}}(s)=\frac{2}{(\lambda+g)s^{3H}}\int_{-\infty}^{1}sdx\;\rho(s(1-x))\times\\\
&\int_{0}^{\frac{g-s}{2}}dt_{1}I(x,y(t_{1}))\end{split}$ (28)

where the upper bound in the second integral takes care of the fact that the
main contribution comes from the nearest loop at the distance less or equal
than $(g-s)/2$ and the factor $2$ stands for the two symmetric positions of
the loop with respect to the points $i,j$ of interest. Switching to the
integration over the universal variable $y=-t_{1}/s+x-1$ we arrive at the
following expression

$P_{a_{0}}(s)\approx\frac{2s}{g}\frac{g-\lambda}{g\;s^{3H}}\int_{-\infty}^{1}sdx\;\rho(s(1-x))\int_{x-1-\frac{g-s}{2s}}^{x-1}dy\;I(x,y)$
(29)

which can be further decomposed into Markovian and non-Markovian contributions
as follows

$\displaystyle
P_{a_{0}}(s)\approx\frac{1}{s^{3H}}\left\\{1-\frac{\lambda}{g}-\frac{s}{g}+\right.$
(30)
$\displaystyle\left.+\left(1-\frac{\lambda}{g}\right)\int_{-\infty}^{1}sdx\;\rho(s(1-x))\;\overline{\Delta
I_{0}(x)}\right\\}$ (31)

where $\overline{\Delta I_{0}(x)}$ is the typical contribution to the contact
probability of the diagram coming from the interdependence of loops and
spacers

$\overline{\Delta
I_{0}(x)}\approx\frac{2s}{g}\int_{x-1-\frac{g}{2s}}^{x-1}dy\left(I(x,y)-1\right)$
(32)

Note that the Markovian factor in (Eq.30) is just the weight of the free
propagator diagram $\frac{g-s}{g+\lambda}\approx 1-\lambda/g-s/g$ in the one-
loop approximation, see (Eq.19). Since the function under the integral in
(Eq.32) is integrable, in the limit $s/g\to 0$ the non-Markovian correction is
of order of $s/g$. This means we can use $-\infty$ as the lower bound of
(Eq.32) and neglect the term $\lambda/g$ in the non-Markovian correction of
(Eq.30) to the leading order. Analogously, the propagators for the diagrams
(b) and (c) can be written using the corresponding expressions for the
conditional probabilities (Eq.15) and the weights (Eq.20) and (Eq.21).

The perturbative form of the contact probability. In the one-loop
approximation, the contact probability ${\cal{P}}^{\text{topo}}(s)$ for the
topological model is determined by the sum of contributions coming from the
diagrams $(a)$, $(b)$ and $(c)$, i.e.

${\cal{P}}^{\text{topo}}(s)=P_{a_{1}}(s)+P_{a_{0}}(s)+P_{b}(s)+P_{c}(s),$ (33)

so that the resulting expression for the averaged contact probability in the
one-loop approximation can be written as follows

${\cal{P}}^{\text{topo}}(s)\propto\left(1+\frac{\lambda}{g}f(s)\right)\frac{1}{s^{3H}}$
(34)

where the first term reflects the contribution of the free propagator and the
function $f(s)$ in the brackets responds for the perturbative contribution of
the one-loop diagrams

$\displaystyle f(s)=-1+\frac{s}{\lambda}\Biggl{(}-1+$ (35) $\displaystyle
2\int_{-\infty}^{0}sdx\;\rho(s(1-x))\int_{-\infty}^{x-1}dy\left(I(x,y)-1\right)+$
(36)
$\displaystyle+\int_{0}^{1}sdx\;\rho(s(1-x))\left(\int_{0}^{x}dy\;I(x,y)+2\int_{x-1}^{0}dy\;I(x,y)\right)$
(37)
$\displaystyle+\int_{-\infty}^{0}sdx\;\rho(s(1-x))\left(\int_{x}^{0}dy\;I(x,y)+2\int_{x-1}^{x}dy\;I(x,y)\right)\Biggr{)}.$
(38)

Taking a log derivative of (Eq.34) and keeping only linear in $\lambda/g$
terms, we arrive at the following perturbative expansion

$\frac{\partial\log{\cal{P}}^{\text{topo}}(s)}{\partial\log
s}=-3H+\frac{\lambda}{g}\Delta(s),$ (39)

where $\Delta(s)=s\frac{\partial f(s)}{\partial s}$. It is straightforward to
show that $\lim_{s\to 0}f(s)=0$ so that the loop-free result gets recovered at
small scales (short polymer subchains do not feel the loop constraints). In
the opposite limit one obtains $\lim_{s\to\infty}f(s)=3H$.

Figure 5: The contact probability and its log-derivative for the non-
topological binding model $P^{\text{non-topo}}(s)$ (gray, dashed), for the
non-topological in the one-loop approximation ${\cal{P}}^{\text{non-topo}}(s)$
(black, dashed) and for the topological in the one-loop approximation
${\cal{P}}^{\text{topo}}(s)$ (black, solid). Loop density is $\lambda/g=1$.
Different fractal dimensions of the chain: (A) $d_{f}=3$, (B) $d_{f}=2.5$, (C)
$d_{f}=2$. The difference between the topological and non-topological models
vanishes as the chain is becoming more ideal.

###  II.3 Difference in the shape of $P(s)$ for the two models

Having the results for the two models, it is important to see how the residual
correlations in the topological binding mechanism change the contact
probability across scales. To understand it, in what follows we will analyze
their impact for the main diagrams (a) and (b). The diagram (c) is not
relevant for us here, since it provides equivalent results for the both
models. The diagram (d) is not included in the one-loop approximation, which
is used to compute the topological model. However, at large contour distances
$s$ it can be approximated by the behaviour of the diagram (a). Indeed,
contribution of the backbone segment in (Eq.10) becomes dominant at
$s\gg\lambda$.

The diagram (a) consists of two diagrams: (a0) and (a1). Fig. 4A shows that
for the two points located on the backbone and one loop in between them
(diagram (a1)) the contact probability is always larger for the non-
topological binding of cohesin, $P^{\text{non-
topo}}_{a1}>P^{\text{topo}}_{a1}$. Importantly, the difference gets stronger
with increase of the loop size $l$ in between. The effect is the most
pronounced when the two points are located close to the loop base, $l/s\approx
1$, i.e. when most of the contour distance is occupied by the loop. The source
of this difference is clear. From the viewpoint of a random walk, the time
elapsed between the points is smaller for the non-topological model as it
effectively goes along the backbone ($s_{eff}$). In other words, the effective
contour distance between the points in the diagram (a1) is shortened by
formation of the loops in between. In the case of topological binding such
shortening is absent for $H<1/2$ and the shortest contour distance equals to
the original distance $s>s_{eff}$. Therefore, the equilibrium distance between
the points is larger in the topological model and the contact probability is
smaller.

However, for $H=1/2$, the difference in $P(s)$ for the diagram (a1) between
the models vanishes, despite the shortest path along the polymer in the
topological model is larger than in the non-topological. One can demonstrate
this as follows. As in the previous section, let us denote by $i,j$ the two
points of interest and by $k,n$ the base of the loop. The vector
$\vec{r}_{ij}=\vec{r}_{i}-\vec{r}_{j}$ can be decomposed as

$\vec{r}_{ij}=\vec{r}_{ik}+\vec{r}_{kn}+\vec{r}_{nj},$ (40)

where $\vec{r}_{kn}=0$ due to the loop. Since in the case of ideal chain all
the vectors in (Eq.40) are independent, the distance between the points $i,j$
is controlled by the distances along the backbone

$(\vec{r}_{ij})^{2}=(\vec{r}_{ik})^{2}+(\vec{r}_{nj})^{2}=(k-i)+(j-n),$ (41)

i.e. the loop is effectively eaten up. Thus, the shortest distance in the non-
topological model becomes the effective contour distance in the topological
for $H=1/2$ for the diagram (a1).

Interestingly, for diagrams (b) and (a0) the situation is the opposite (see
Fig. 4B,C): in the topological model the contact probability is larger,
$P^{\text{topo}}_{a0,b}>P^{\text{non-topo}}_{a0,b}$, in the case of non-ideal
chain. Similarly to the diagram (a1), upon increase of the loop size $l$ the
effect gets stronger. However, the physics of the difference between the
models in this case is not related with the chain shortening, but with
residual correlations that are not destroyed by formation of a loop via the
topological mechanism. Indeed, for example, for the diagram (b) residual
correlations between a point in the loop and a point in the spacer persist in
the topological mode of loop formation. Since these correlations are negative
(for $H<1/2$), the equilibrium distance between the points is shorter than in
the non-topological case, where the points fluctuate independently.
Accordingly, for the ideal chain statistics ($H=1/2$) the two points in either
model fluctuate independently, and the contact probability is the same.
Interestingly, for sufficiently large loops, $l\gg s$, the maximal fold-
difference between the models corresponds to the symmetric positions of the
points relatively to the loop base, $t/s=1/2$ (see Fig. 4C). However, upon
decrease of the loop size $l$ the optimal relative distance to the base drops
down to $t/s=1/4$ at $l=s$, meaning that the point in the spacer should be
three times further from the loop base than the one in the loop. Such an
asymmetry in the configuration is related with different statistical
properties of spacer and loop regions, which becomes noticeable when their
contour sizes get similar.

Despite the topological model can be solved only in the one-loop approximation
(i.e. for $\lambda/g\ll 1$), one can understand the qualitative difference
between the models for all $\lambda/g$. Indeed, the conformational statistics
of a loopy chain at $s\gg\lambda$ is determined by the distance between the
points located on the main chain, i.e. by the diagram (a), see Fig. 2(a). This
distance is effectively shortened in the case of non-topological binding
model, thus leading to a higher contact probability, similarly to the one-loop
diagram (a1), see Fig. 4A. With increase of the separation $s$ between the
points the fraction of the chain within the backbone (along the shortest path)
increases to the bulk value $\frac{g}{\lambda+g}$. For the reference loop
density $\lambda/g=1$, this corresponds to half of the contour distance
effectively shortened. For the case of a single loop diagram (a1), this yields
around $\approx 15\%$ for the difference in the absolute values of $P(s)$ at
large scales (Fig. 4a). Therefore, with increase of loop density $\lambda/g$
the drop of $P^{\text{topo}}$ compared to $P^{\text{non-topo}}$ at large
scales increases, in particular, resulting in a deeper dip of the log-
derivative.

On the other hand, at short scales $s\approx\lambda$ all diagrams depicted in
Fig. 4 contribute similarly to $P(s)$. While the diagram (a1) with one loop in
between yields a larger value of $P(s)$ in case of non-topological binding,
this is not the case for diagrams (b) and (a0). Thus, we expect the difference
between the models at short scales to be less pronounced.

This analysis is confirmed by Fig. 5, where for the sake of accuracy we
compare the two models in the one-loop approximation
${\cal{P}}^{\text{topo}}(s)$ and ${\cal{P}}^{\text{non-topo}}(s)$, as well as
the full theory for the non-topological model, $P^{\text{non-topo}}(s)$. For
that we use the expressions (Eq.5)-(Eq.7) for the single-loop diagrams in the
non-topological case and insert them as the particular contributions for
calculation of ${\cal{P}}^{\text{non-topo}}(s)$. The averaging of the single-
loops diagrams, as for the topological model, is conducted using the
linearized weights (Eq.18)-(Eq.21). First we note that at intermediate loop
densities ($\lambda/g=1$) the one-loop approximation produces larger
amplitudes of the peak and the dip on the log-derivative, for all fractal
dimensions $d_{f}$. Interestingly, the effect of the topological binding at
the loop base is qualitatively similar. Indeed, in the full agreement with our
analysis above, if one compares ${\cal{P}}^{\text{topo}}(s)$ and
${\cal{P}}^{\text{non-topo}}(s)$ (black lines in Fig. 5), it is evident that
the topological model yields a smaller contact probability at larger scales
and, thus, larger amplitude of the dip. The corresponding difference at short
scales $s\leq\lambda$ is much less pronounced than at large scales. Still one
sees that the topological model at short scales produces somewhat larger and
slower decaying ${\cal{P}}^{\text{topo}}$, especially for larger values of the
fractal dimension. Upon decrease of the fractal dimension from $d_{f}=3$ to
$d_{f}=2$ the two models in the one-loop approximation give more similar
results, consistently with our analysis.

All together, our analysis shows that the main difference between the binding
models concerns large scales. With increase of the loop density, the
topological model produces smaller $P^{\text{topo}}(s)$ than the non-
topological $P^{\text{non-topo}}(s)$ at $s\gg\lambda$, resulting in a deeper
dip of the log-derivative. Physically, this difference arises from the
effective chain shortening that is established by the loops in the non-
topological binding model.

##  III Discussion

Chromosomes are polymers organized into loops by SMC complexes (e.g.
condensins, cohesins) [1]. These loops are transient, as they are produced in
a process of loop extrusion of a SMC complex that stochastically binds and
unbinds DNA [6]. Therefore, these random loops are not directly visible in the
data (Hi-C or microscopy). Here and in our previous paper [7] we have
suggested an analytically tractable polymer model, in which the presence of
random loops, as well as their statistical characteristics, are encoded in the
specific shape of the average contact probability curve. At the same time,
this $P(s)$ function can be also computed from the experimental Hi-C contact
maps [11], providing a robust framework to characterize the loopy state of a
chromosome.

Despite the loop extrusion by cohesin (and, generally, by SMC complexes) has
been confirmed in in vitro experiments [8, 9, 24, 25, 26], the microscopic
mechanism of binding and reeling on DNA remains contradictory. As other SMC
complexes, cohesin has a ring-shaped structure, but whether this ring is
actually needed to topologically embrace two strands of DNA in the course of
extrusion is a subject of ongoing experimental and computational research [10,
2]. It has been proposed three possible binding mechanisms [10]: (i)
topological, i.e. the cohesin is bound to one string chemically, while the
ring opens and closes when cohesin needs to encircle the second strand of DNA;
(ii) pseudo-topological, when both strands are embraced by the ring without
chemical binding in the loop base, so that cohesin can be easily pulled away
from the chromosome; (iii) non-topological, when the two DNA strands are
chemically bound to the ring, which implies the cohesin ring does not play an
important role in binding the DNA. In particular, in the mechanisms (i) and
(ii) the chromatin strand is reeled into the cohesin ring in the course of
extrusion together with the proteins residing on the chain. Clearly,
translocation of obstacles larger than the ring size ($\approx 35$nm) into the
loop would not be possible. However, a recent single-molecule study has shown
that nanoparticles as large as $\approx 200$nm, as well as big protein
machines (e.g., RNA polymerase), can be easily translocated into the loop by
cohesin [10]. This experimentally observed ability of cohesin to avoid large
obstacles on DNA points to the non-topological binding mechanism. However,
other authors have previously noted that the topological binding could explain
atypically long residence times of cohesin on chromatin ($\approx 10$
minutes), as compared to other molecules without the ring structure [2].
Furthermore, cohesin was previously shown to topologically entrap two sister
chromatids [27]. Thus, the role of the cohesin ring structure in DNA binding
and extrusion of loops is not yet fully established.

In this paper we underscore the importance of a particular microscopic
structure of the cohesin loop base on the 3D organization of a chromosome
across scales. Standing on the experimental observation that without cohesin
loops chromosomes follow the statistics of a crumpled chain with fractal
dimension $d_{f}=3$, and not of an ideal chain with $d_{f}=2$ [14, 15, 7], we
demonstrate that the conformational properties of a crumpled chain with loops
would be different for the topological and non-topological structure of the
loop base. While for the non-topological binding the backbone of a loopy
polymer forms with the loops emanating out of it (as in a comb-like polymer),
in the topological case the polymer folds into loops, but its organization
remains linear. This leads to the two main peculiarities of the two models:
(i) effective shortening of the chromosome at large scales in the non-
topological model, and (ii) negatively correlated fluctuations of adjacent
loops and spacers in the topological model. Notably, for the ideal chain
statistics ($d_{f}=2$) the correlations along the chain are absent, resulting
in equivalent organization of the loopy chain in either binding model.
However, for any $d_{f}>2$, due to strong (power-law) negative correlations,
the mode of cohesin binding has a quantitative effect on the three-dimensional
organization of the chains into loops.

In order to quantify this effect we used an analytically tractable polymer
model [7] and compare the behaviour of the contact probability $P(s)$ for
topological and non-topological binding models. We based our theoretical
argument on the quenched disorder approach and made several assumptions that
allowed us to solve the non-topological model analytically for any fractal
dimension $d_{f}$ and any density of loops $\lambda/g$ (see [7]). Under the
same set of assumptions, the topological model is not analytically solvable,
so here we employed the one-loop approximation which treats the loops as a
perturbation to the fractal scaling of the chain. We obtained the perturbative
form of $P(s)$ in the loop density $\lambda/g$ and analytically computed the
linear order correction using the Gaussian measure of fBm trajectories as
suggested recently. This result allowed us to quantitatively compare the
contact probability between the two models for the single-loop diagrams ((a),
(b) and (c)). We found that the binding mechanism in the loop base has indeed
the impact on the $P(s)$, which is qualitatively similar to the difference
between the one-loop approximation and the full theory within the single
binding mechanism.

The revealed quantitative difference for the two binding models allows to
rationalize the organization of chromosomes beyond the minimal analytically
tractable models. We note the central role of the non-ideal (crumpled)
statistics of the chain in the observed difference. Thus, for rare species
where the emergent statistics of the chromosomes is closer to the ideal chain
(e.g. yeast [28]), our analysis implies universal contact probability
behaviour independent of the microscopic details imposed by the SMC machinery.

##  IV Acknowledgements

We thank Leonid Mirny, Sergei Nechaev, Mehran Kardar, Hugo Brandao, Vittore
Scolari and Sergei Belan for valuable discussions on the subject of the paper.
The work is supported by the Russian Science Foundation (Grant No.
21-73-00176).

## References

## References

  * Mirny and Solovei [2021] L. Mirny and I. Solovei, Keeping chromatin in the loop (s), Nature Reviews Molecular Cell Biology 22, 439 (2021).
  * Banigan and Mirny [2020] E. J. Banigan and L. Mirny, Loop extrusion: theory meets single-molecule experiments, Current opinion in cell biology 64, 124 (2020).
  * Goloborodko  _et al._ [2016] A. Goloborodko, J. F. Marko, and L. A. Mirny, Chromosome compaction by active loop extrusion, Biophysical journal 110, 2162 (2016).
  * Paulson and Laemmli [1977] J. Paulson and U. Laemmli, The structure of histone-depleted metaphase chromosomes, Cell 12, 817 (1977).
  * Gibcus  _et al._ [2018] J. Gibcus _et al._ , A pathway for mitotic chromosome formation, Science 359, 2133 (2018).
  * Fudenberg  _et al._ [2017] G. Fudenberg, N. Abdennur, M. Imakaev, A. Goloborodko, and L. Mirny, Emerging evidence of chromosome folding by loop extrusion, Cold Spring Harbor symposia on quantitative biology , 45 (2017).
  * Polovnikov  _et al._ [2022] K. Polovnikov _et al._ , Fractal polymer with loops recapitulates key features of chromosome organization, bioRxiv (2022).
  * Terakawa  _et al._ [2017] T. Terakawa _et al._ , The condensin complex is a mechanochemical motor that translocates along dna, Science 358, 672 (2017).
  * Ganji  _et al._ [2018] M. Ganji _et al._ , Real-time imaging of dna loop extrusion by condensin, Science 360, 102 (2018).
  * Pradhan  _et al._ [2022] B. Pradhan _et al._ , Smc complexes can traverse physical roadblocks bigger than their ring size, Cell Reports 41, 111491 (2022).
  * Lieberman-Aiden  _et al._ [2009] E. Lieberman-Aiden _et al._ , Comprehensive mapping of long-range interactions reveals folding principles of the human genome, Science 326, 289 (2009).
  * Mirny [2011] L. Mirny, The fractal globule as a model of chromatin architecture in the cell, Chromosome Res. 19, 37 (2011).
  * Halverson  _et al._ [2014] J. D. Halverson, J. Smrek, K. Kremer, and A. Y. Grosberg, From a melt of rings to chromosome territories: the role of topological constraints in genome folding, Reports on Progress in Physics 77, 022601 (2014).
  * Hsieh  _et al._ [2021] T. Hsieh _et al._ , Enhancer-promoter interactions and transcription are maintained upon acute loss of CTCF, cohesin, WAPL, and YY1, BioRxiv (2021).
  * Rao  _et al._ [2017] S. Rao _et al._ , Cohesin loss eliminates all loop domains, Cell 171, 305 (2017).
  * Polovnikov  _et al._ [2018] K. Polovnikov, S. Nechaev, and M. V. Tamm, Effective Hamiltonian of topologically stabilized polymer states, Soft Matter 14, 6561 (2018).
  * Grosberg and Khokhlov [1994] A. Y. Grosberg and A. R. Khokhlov,  _Statistical Mechanics of Macromolecules_ (AIP: Woodbury, NY., 1994) p. 350.
  * De Gennes [1979] P. G. De Gennes, Scaling Concepts in Polymer Physics, Cornell university press (1979).
  * Grosberg  _et al._ [1988] A. Grosberg, S. Nechaev, and E. Shakhnovich, The role of topological constraints in the kinetics of collapse of macromolecules, Journal de physique 49, 2095 (1988).
  * Grosberg  _et al._ [1993] A. Grosberg, Y. Rabin, S. Havlin, and A. Neer, Crumpled globule model of the three-dimensional structure of DNA, Europhysics Letters 23, 373 (1993).
  * Nechaev  _et al._ [1987] S. Nechaev, A. Semenov, and M. Koleva, Dynamics of a polymer chain in an array of obstacles, Physica A: Statistical Mechanics and its Applications 140, 506 (1987).
  * Pedler [1971] P. Pedler, Occupation times for two-state markov chains, J.Appl.Prob. 8, 381 (1971).
  * Polovnikov  _et al._ [2019] K. Polovnikov, S. Nechaev, and M. V. Tamm, Many-body contacts in fractal polymer chains and fractional Brownian trajectories, Physical Review E 99, 032501 (2019).
  * Kim  _et al._ [2019] Y. Kim _et al._ , Human cohesin compacts dna by loop extrusion, Science 366, 1345 (2019).
  * Golfier  _et al._ [2020] S. Golfier _et al._ , Cohesin and condensin extrude dna loops in a cell cycle-dependent manner, Elife 9, e53885 (2020).
  * Davidson  _et al._ [2019] I. Davidson _et al._ , Dna loop extrusion by human cohesin, Science 366, 1338 (2019).
  * Ivanov and Nasmyth [2005] D. Ivanov and K. Nasmyth, A topological interaction between cohesin rings and a circular minichromosome, Cell 122, 849 (2005).
  * Wong  _et al._ [2012] H. Wong _et al._ , A predictive computational model of the dynamic 3d interphase yeast nucleus, Current biology 22, 1881 (2012).