Mode selectivity of dynamically induced conformation in many-body chain-like bead-spring model

The Lagrangian of the system is

$$L(\boldsymbol{r},\dot{\boldsymbol{r}})=K(\dot{\boldsymbol{r}})-V(\boldsymbol{r})$$

(1)

in the Cartesian coordinates, where

$$\boldsymbol{r}=\begin{pmatrix}\boldsymbol{r}_{1}\\ \vdots\\ \boldsymbol{r}_{N}\end{pmatrix},\quad\dot{\boldsymbol{r}}=\begin{pmatrix}\dot{\boldsymbol{r}}_{1}\\ \vdots\\ \dot{\boldsymbol{r}}_{N}\\ \end{pmatrix}\quad\in\mathbb{R}^{2N}.$$

(2)

The kinetic energy function $K(\dot{\boldsymbol{r}})$ is

$$K(\dot{\boldsymbol{r}})=\dfrac{1}{2}\sum_{i=1}^{N}m_{i}||\dot{\boldsymbol{r}}_{i}||^{2}.$$

(3)

The potential energy function $V(\boldsymbol{r})$ consists of the spring potential $V_{\rm spring}$ and the bending potential $U_{\rm bend}$ as

$$V(\boldsymbol{r})=V_{\rm spring}(\boldsymbol{r})+U_{\rm bend}(\boldsymbol{r}),$$

(4)

where we assume that $V_{\rm spring}$ depends on only the lengths of the springs. When we consider a molecule, the spring potential represents stronger bonds, and the bending potential corresponds to weaker bonds. A chain-like model is expressed by the nearest neighbor interactions in the spring potential as

$$V_{\rm spring}(\boldsymbol{r})=\sum_{i=1}^{N-1}V_{i}(||\boldsymbol{r}_{i+1}-\boldsymbol{r}_{i}||).$$

(5)

The system described by the Lagrangian (1) has the two-dimensional translational symmetry. To reduce the two degrees of freedom, we introduce $2N-2$ internal coordinates $\boldsymbol{y}\in\mathbb{R}^{2N-2}$ as

$$\boldsymbol{y}=\begin{pmatrix}\boldsymbol{l}\\ \boldsymbol{\phi}\\ \end{pmatrix}\in\mathbb{R}^{2N-2},\quad\left\{\begin{array}[]{l}\boldsymbol{l}=(l_{1},\cdots,l_{N-1})^{\rm T}\in\mathbb{R}^{N-1}\\ \boldsymbol{\phi}=(\phi_{1},\cdots,\phi_{N-1})^{\rm T}\in\mathbb{R}^{N-1}\\ \end{array}\right.$$

(6)

where the superscript T represents transposition. The internal coordinates play an important role to extract coupling between bending motion and spring motion (see Ref. yanao-etal-07 for instance).

The vector $\boldsymbol{l}$ contains the lengths of the springs as

$$l_{i}=||\boldsymbol{r}_{i+1}-\boldsymbol{r}_{i}||,\quad(i=1,\cdots,N-1).$$

(7)

The Euclidean norm is defined by $||\boldsymbol{x}||=\sqrt{x^{2}+y^{2}}$ for $\boldsymbol{x}=(x,y)\in\mathbb{R}^{2}$.

The angle $\phi_{i}~{}(i=1,\cdots,N-2)$ represents the bending angle between the vectors $\boldsymbol{r}_{i+2}-\boldsymbol{r}_{i+1}$ and $\boldsymbol{r}_{i+1}-\boldsymbol{r}_{i}$, and satisfies

$$\cos\phi_{i}=\dfrac{(\boldsymbol{r}_{i+2}-\boldsymbol{r}_{i+1})\cdot(\boldsymbol{r}_{i+1}-\boldsymbol{r}_{i})}{||\boldsymbol{r}_{i+2}-\boldsymbol{r}_{i+1}||~{}||\boldsymbol{r}_{i+1}-\boldsymbol{r}_{i}||},\quad(i=1,\cdots,N-2)$$

(8)

where $\boldsymbol{x}\cdot\boldsymbol{y}$ is the Euclidean inner product between $\boldsymbol{x}$ and $\boldsymbol{y}$. The last angle variable $\phi_{N-1}$ associates to the total angular momentum. The system has the rotational symmetry and $\phi_{N-1}$ is a cyclic coordinate, but we keep it for later convenience.

Rewriting the kinetic energy in the variables $\boldsymbol{y}$ and $\dot{\boldsymbol{y}}$, and assuming that the total angular momentum is zero, we have the Lagrangian

$$L(\boldsymbol{y},\dot{\boldsymbol{y}})=\dfrac{1}{2}\sum_{\alpha,\beta=1}^{2N-2}B^{\alpha\beta}(\boldsymbol{y})\dot{y}_{\alpha}\dot{y}_{\beta}-V(\boldsymbol{y}).$$

(9)

We used the same symbol $V(\boldsymbol{y})$ for the potential energy function to simplify the notation. The function $B^{\alpha\beta}(\boldsymbol{y})$ is the $(\alpha,\beta)$ element of the matrix $\boldsymbol{\rm B}(\boldsymbol{y})\in{\rm Mat}(2N-2)$. Here ${\rm Mat}(n)$ represents the set of real square matrices of size $n$. The explicit form of $\boldsymbol{\rm B}(\boldsymbol{y})$ is given in Appendix A. Note that an alphabetic index runs from $1$ to $N-1$, and a Greek index runs from $1$ to $2N-2$.

From now on, we adopt the Einstein notation for the sum: We take the sum over an index if it appears twice in a term. The Euler-Lagrange equation for the variable $y_{\alpha}~{}(\alpha=1,\cdots,2N-2)$ is expressed as

$$B^{\alpha\beta}(\boldsymbol{y})\ddot{y}_{\beta}+D_{\alpha}^{\beta\gamma}(\boldsymbol{y})\dot{y}_{\beta}\dot{y}_{\gamma}+\dfrac{\partial V}{\partial y_{\alpha}}(\boldsymbol{y})=0$$

(10)

where

$$D_{\alpha}^{\beta\gamma}(\boldsymbol{y})=\dfrac{\partial B^{\alpha\beta}}{\partial y_{\gamma}}(\boldsymbol{y})-\dfrac{1}{2}\dfrac{\partial B^{\beta\gamma}}{\partial y_{\alpha}}(\boldsymbol{y}).$$

(11)

Let $\boldsymbol{l}_{\ast}=(l_{1,\ast},\cdots,l_{N-1,\ast})^{\rm T}$ be the natural length vector of the springs, namely $\boldsymbol{l}_{\ast}$ solves

$$\dfrac{\partial V_{\rm spring}}{\partial\boldsymbol{l}}(\boldsymbol{l}_{\ast})=\boldsymbol{0}.$$

(12)

We introduce the two assumptions with a small parameter $\epsilon$:

• (A1)

  The amplitudes of the springs are sufficiently small comparing with the natural lengths. The ratio is of $O(\epsilon)$.

• (A2)

  Large bending motion is sufficiently slow than the spring motion. The ratio of the two timescales is of $O(\epsilon)$.

We express these assumptions by the expansions of $t,\boldsymbol{l}$, and $\boldsymbol{\phi}$ as

$$\dfrac{{\rm d}}{{\rm d}t}=\dfrac{\partial}{\partial t_{0}}+\epsilon\dfrac{\partial}{\partial t_{1}}$$

(13)

and

$$\begin{split}&\boldsymbol{l}(t_{0},t_{1})=\boldsymbol{l}_{\ast}+\epsilon\boldsymbol{l}^{(1)}(t_{0},t_{1}),\\ &\boldsymbol{\phi}(t_{0},t_{1})=\boldsymbol{\phi}^{(0)}(t_{1})+\epsilon\boldsymbol{\phi}^{(1)}(t_{0},t_{1}),\\ \end{split}$$

(14)

which is summarized as

$$\boldsymbol{y}(t_{0},t_{1})=\boldsymbol{y}^{(0)}(t_{1})+\epsilon\boldsymbol{y}^{(1)}(t_{0},t_{1}).$$

(15)

The two timescales $t_{0}=t$ and $t_{1}=\epsilon t$ correspond to the fast spring motion and the slow bending motion, respectively. The vector $\boldsymbol{l}_{\ast}=(l_{1,\ast},\cdots,l_{N-1,\ast})^{\rm T}$ is constant. We are interested in the large and slow motion of the bending angles described by $\boldsymbol{\phi}^{(0)}(t_{1})$.

Two remarks are in order. First, the leading order of the expanded equations of motion is of $O(\epsilon)$, since the leading order of the velocities and the accelerations is of $O(\epsilon)$. Second, the velocities $\dot{l}_{i}$ and $\dot{\phi}_{i}$ are of the same order of $O(\epsilon)$, and the fastness of motion in the assumption (A2) connotes a short period. Indeed, the distance of a normal mode orbit is of $O(\epsilon)$ and the period is of $O(\epsilon^{0})$, while the distance of the large bending motion is of $O(\epsilon^{0})$ and the period is of $O(\epsilon^{-1})$.

We also expand the potential functions $V_{\rm spring}$ and $U_{\rm bend}$ into power series of $\epsilon$. The expansions of the two functions are performed in different ways, as a result of Eq. (14).

The spring potential function $V_{\rm spring}$ is expanded around $\boldsymbol{l}_{\ast}$ as

$$V_{\rm spring}(\boldsymbol{l})=V_{\rm spring}(\boldsymbol{l}_{\ast})+\dfrac{1}{2}(\boldsymbol{l}-\boldsymbol{l}_{\ast})^{\rm T}\boldsymbol{\rm K}_{l}(\boldsymbol{l}-\boldsymbol{l}_{\ast})+O(\epsilon^{3}),$$

(16)

where the $(i,j)$ element $K_{l}^{ij}$ of the matrix $\boldsymbol{\rm K}_{l}\in{\rm Mat}(N-1)$ is

$$K_{l}^{ij}=\dfrac{\partial^{2}V_{\rm spring}}{\partial l_{i}\partial l_{j}}(\boldsymbol{l}_{\ast}),\quad(i,j=1,\cdots,N-1)$$

(17)

and $\boldsymbol{\rm K}_{l}$ is assumed to be positive definite.

The bending potential $U_{\rm bend}$ can be also expanded around a stationary point $\boldsymbol{\phi}_{\ast}$ if it exists, but this expansion is not useful since $||\boldsymbol{\phi}-\boldsymbol{\phi}_{\ast}||$ is not necessarily of $O(\epsilon)$. We thus expand $U_{\rm bend}$ by its amplitude as

$$U_{\rm bend}(\boldsymbol{y})=U_{\rm bend}^{(0)}(\boldsymbol{y})+\epsilon U_{\rm bend}^{(1)}(\boldsymbol{y})+\epsilon^{2}U_{\rm bend}^{(2)}(\boldsymbol{y})+\cdots.$$

(18)

It is important to note that $U_{\rm bend}^{(0)}(\boldsymbol{y})$ and $U_{\rm bend}^{(1)}(\boldsymbol{y})$ solely depend on $\boldsymbol{l}$ with the conditions

$$\dfrac{\partial U_{\rm bend}^{(0)}}{\partial\boldsymbol{l}}(\boldsymbol{l}_{\ast})=\dfrac{\partial U_{\rm bend}^{(1)}}{\partial\boldsymbol{l}}(\boldsymbol{l}_{\ast})=\boldsymbol{0}$$

(19)

for satisfying the assumptions (see Appendix B). Integrating them into the spring potential $V_{\rm spring}(\boldsymbol{l})$, we may set $U_{\rm bend}^{(0)}(\boldsymbol{y}),U_{\rm bend}^{(1)}(\boldsymbol{y})\equiv 0$. The leading order of $U_{\rm bend}$ is hence of $O(\epsilon^{2})$, and this ordering is consistent with weaker bonds derived from the bending potential.

Substituting Eqs. (13), (14), (16), and (18) into Eq. (10), we have the expanded equations of motion order by order of $\epsilon$. As remarked in the end of Sec. III.1, the nontrivial equations of motion start from $O(\epsilon)$.

The equations (23) depend on the fast timescale $t_{0}$ through $\boldsymbol{y}^{(1)}$, and we eliminate it by taking the average. The average in the timescale $t_{0}$ is defined by

$$\left\langle\varphi\right\rangle=\lim_{T\to\infty}\dfrac{1}{T}\int_{0}^{T}\varphi(t_{0})dt_{0}.$$

(25)

The averaged equations are

$$\begin{split}&B^{\alpha\beta}(\boldsymbol{y}^{(0)})\dfrac{{\rm d}{}^{2}y_{\beta}^{(0)}}{{\rm d}t_{1}^{2}}+D_{\alpha}^{\beta\gamma}(\boldsymbol{y}^{(0)})\dfrac{{\rm d}y_{\beta}^{(0)}}{{\rm d}t_{1}}\dfrac{{\rm d}y_{\gamma}^{(0)}}{{\rm d}t_{1}}\\ &+\dfrac{\partial U_{\rm bend}^{(2)}}{\partial y_{\alpha}}(\boldsymbol{y}^{(0)})=\mathcal{A}_{\alpha},\end{split}$$

(26)

where the averaged term $\mathcal{A}_{\alpha}$ is

$$\begin{split}\mathcal{A}_{\alpha}&=\dfrac{1}{2}{\rm Tr}\left[\dfrac{\partial\boldsymbol{\rm B}}{\partial y_{\alpha}}(\boldsymbol{y}^{(0)})\boldsymbol{\rm X}(\boldsymbol{y}^{(0)})\left\langle\boldsymbol{y}^{(1)}\boldsymbol{y}^{(1){\rm T}}\right\rangle\right].\\ \end{split}$$

(27)

The symbol ${\rm Tr}$ represents the matrix trace. We used the relation

$$\left\langle\dfrac{\partial{\boldsymbol{y}}^{(1)}}{\partial t_{0}}\left(\dfrac{\partial\boldsymbol{y}^{(1)}}{\partial t_{0}}\right)^{\rm T}\right\rangle=\boldsymbol{\rm X}(\boldsymbol{y}^{(0)})\left\langle\boldsymbol{y}^{(1)}\boldsymbol{y}^{(1){\rm T}}\right\rangle$$

(28)

proven by performing the integration by parts.

The averaged term $\mathcal{A}_{\alpha}$ depends on the solution $\boldsymbol{y}^{(1)}$ to Eq. (20), which is obtained by diagonalizing the matrix $\boldsymbol{\rm X}(\boldsymbol{y}^{(0)})$ as

$$\boldsymbol{\rm X}(\boldsymbol{y}^{(0)})\boldsymbol{\rm P}(\boldsymbol{y}^{(0)})=\boldsymbol{\rm P}(\boldsymbol{y}^{(0)})\boldsymbol{\rm\Lambda}(\boldsymbol{y}^{(0)}).$$

(29)

$\boldsymbol{\rm P}\in{\rm Mat}(2N-2)$ is a diagonalizing matrix, and the diagonal matrix $\boldsymbol{\rm\Lambda}$ contains the eigenvalues of $\boldsymbol{\rm X}$:

$$\boldsymbol{\rm\Lambda}(\boldsymbol{y}^{(0)})={\rm diag}(\lambda_{1},\cdots,\lambda_{N-1},0,\cdots,0)\in{\rm Diag}(2N-2),$$

(30)

where the symbol ${\rm Diag}(n)$ represents the set of real diagonal matrices. The $N-1$ zeroeigenvalues come from the fact $\dim({\rm Ker}\boldsymbol{\rm K})=N-1$ in general. We assume that the corresponding $N-1$ amplitudes of $\boldsymbol{y}^{(1)}$ are zero. Since the average of the square of a sinusoidal function is $1/2$, we have in general

$$\left\langle\boldsymbol{y}^{(1)}\boldsymbol{y}^{(1){\rm T}}\right\rangle=\dfrac{1}{2}\boldsymbol{\rm P}(\boldsymbol{y}^{(0)})\boldsymbol{\rm W}(t_{1})^{2}\boldsymbol{\rm P}(\boldsymbol{y}^{(0)})^{\rm T},$$

(31)

where

$$\boldsymbol{\rm W}(t_{1})={\rm diag}(w_{1},\cdots,w_{N-1},0,\cdots,0)\in{\rm Diag}(2N-2).$$

(32)

The matrix $\boldsymbol{\rm W}$ contains the $N-1$ nontrivial amplitudes $w_{i}~{}(i=1,\cdots,N-1)$, which evolve in the slow timescale $t_{1}$ through coupling with $\boldsymbol{\phi}^{(0)}(t_{1})$, while the initial phases of $\boldsymbol{y}^{(1)}$ are eliminated by the average.

Summarizing, the averaged term $\mathcal{A}_{\alpha}$ is modified by Eqs. (29) and (31) to

$$\begin{split}&\mathcal{A}_{\alpha}(\boldsymbol{y}^{(0)})=\dfrac{1}{4}{\rm Tr}\left[\boldsymbol{\rm P}^{\rm T}\dfrac{\partial\boldsymbol{\rm B}}{\partial y_{\alpha}}\boldsymbol{\rm P}\boldsymbol{\rm\Lambda}\boldsymbol{\rm W}(t_{1})^{2}\right],\\ &(\alpha=1,\cdots,2N-2).\end{split}$$

(33)

We can show that

$$\mathcal{A}_{1}=\cdots=\mathcal{A}_{N-1}=0.$$

(34)

In other words, the averaged terms survive only in the equations for $\boldsymbol{\phi}^{(0)}$. See Appendix C.

We stress that Eq. (26) is not closed, because the averaged term (33) depends on $\boldsymbol{\rm W}(t_{1})$, which nontrivially evolves in the timescale $t_{1}$. We have to eliminate the $N-1$ nontrivial amplitudes $w_{i}(t_{1})~{}(i=1,\cdots,N-1)$. To this end, we introduce the hypothesis yamaguchi-etal-22 inspired from the adiabatic invariance:

$$w_{i}(t_{1})^{2}=\nu_{i}w(t_{1})^{2}\quad(i=1,\cdots,N-1).$$

(35)

This hypothesis supposes that the normal mode energy ratios are constant of time, and reduces the number of unknown variables from $N-1$ to one, although the averaged term $\mathcal{A}_{\alpha}$ depends on the constants $\nu_{i}~{}(i=1,\cdots,N-1)$. We arrange the constants in the matrix

$$\boldsymbol{\rm N}={\rm diag}(\nu_{1},\cdots,\nu_{N-1},0,\cdots,0)\in{\rm Diag}(2N-2),$$

(36)

and the amplitude matrix $\boldsymbol{\rm W}$ is simplified to

$$\boldsymbol{\rm W}(t_{1})^{2}=w(t_{1})^{2}\boldsymbol{\rm N}.$$

(37)

The averaged term $\mathcal{A}_{\alpha}$ is then modified to

$$\mathcal{A}_{\alpha}(\boldsymbol{y}^{(0)})=\dfrac{w(t_{1})^{2}}{4}{\rm Tr}\left[\boldsymbol{\rm P}^{\rm T}\dfrac{\partial\boldsymbol{\rm B}}{\partial y_{\alpha}}\boldsymbol{\rm P}\boldsymbol{\rm\Lambda}\boldsymbol{\rm N}\right].$$

(38)

We remark that the numbering of the normal modes is important in the application of the hypothesis, since the eigenvalues of $\boldsymbol{\rm X}(\boldsymbol{y}^{(0)})$ depend on the bending angles $\boldsymbol{\phi}^{(0)}$. In $N=3$, the two springs expand and contract simultaneously (alternatively) in the in-phase mode (the antiphase mode), which has the energy ratio $\nu_{\rm in}$ ($\nu_{\rm anti}$). Adopting $\nu_{1}=\nu_{\rm in}$ and $\nu_{2}=\nu_{\rm anti}$, the hypothesis (35) is approximately verified for any value of $\boldsymbol{\phi}^{(0)}$ yamaguchi-etal-22 .

However, the global numbering is not trivial in general, because two eigenvalues of $\boldsymbol{\rm X}(\boldsymbol{y}^{(0)})$ may cross by varying $\boldsymbol{\phi}^{(0)}$ as shown in Fig. 2. Nevertheless, the hypothesis is approximately valid when the bending motion is not large, since the system is close to the integrable system, Eq. (20), and the mode numbers can be identified in a local region of $\boldsymbol{\phi}^{(0)}$. Consequently, the hypothesis is useful to study stationarity and stability of a conformation. From now on, we locally number the modes in the ascending order of the eigenvalues (see Fig. 2) unless there otherwise stated.

The last unknown variable $w(t_{1})$ is eliminated by the energy conservation. Expanding energy $E$ as $E=\epsilon^{2}E^{(2)}+\cdots$, the leading term is written as

$$\begin{split}E^{(2)}&=\dfrac{1}{2}{\rm Tr}\left[\boldsymbol{\rm B}(\boldsymbol{y}^{(0)})\left(\dot{\boldsymbol{y}}\right)^{(1)}\left(\dot{\boldsymbol{y}}\right)^{(1){\rm T}}\right]+U_{\rm bend}^{(2)}(\boldsymbol{y}^{(0)})\\ &+\dfrac{1}{2}{\rm Tr}\left[\boldsymbol{\rm K}\boldsymbol{y}^{(1)}\boldsymbol{y}^{(1){\rm T}}\right].\end{split}$$

(39)

Taking the average, we have

$$\begin{split}\left\langle E^{(2)}\right\rangle&=\dfrac{1}{2}{\rm Tr}\left[\boldsymbol{\rm B}(\boldsymbol{y}^{(0)})\dfrac{{\rm d}\boldsymbol{y}^{(0)}}{{\rm d}t_{1}}\left(\dfrac{{\rm d}\boldsymbol{y}^{(0)}}{{\rm d}t_{1}}\right)^{\rm T}~{}\right]+U_{\rm bend}^{(2)}(\boldsymbol{y}^{(0)})\\ &+\dfrac{1}{2}{\rm Tr}\left[\boldsymbol{\rm P}^{\rm T}\boldsymbol{\rm K}\boldsymbol{\rm P}\boldsymbol{\rm W}^{2}\right]\\ \end{split}$$

(40)

by the relations (28) and (31). In the right-hand side of Eq. (40), the first line represents the bending energy, and the second line the spring energy.

The hypothesis (36) gives the equality

$$\begin{split}w(t_{1})^{2}{\rm Tr}\left[\boldsymbol{\rm P}^{\rm T}\boldsymbol{\rm K}\boldsymbol{\rm P}\boldsymbol{\rm N}\right]&=2\left[E^{(2)}-U_{\rm bend}^{(2)}(\boldsymbol{y}^{(0)})\right]\\ &-{\rm Tr}\left[\boldsymbol{\rm B}(\boldsymbol{y}^{(0)})\dfrac{{\rm d}\boldsymbol{y}^{(0)}}{{\rm d}t_{1}}\left(\dfrac{{\rm d}\boldsymbol{y}^{(0)}}{{\rm d}t_{1}}\right)^{\rm T}~{}\right],\end{split}$$

(41)

where we denoted $\left\langle E^{(2)}\right\rangle$ by $E^{(2)}$ for simplicity. Substituting Eq. (41) into Eq. (38), we have the averaged term $\mathcal{A}_{\alpha}$ as

$$\begin{split}\mathcal{A}_{\alpha}(\boldsymbol{y}^{(0)})&=\left\{\dfrac{E^{(2)}-U_{\rm bend}^{(2)}(\boldsymbol{y}^{(0)})}{2}\right.\\ &\left.-\dfrac{1}{4}{\rm Tr}\left[\boldsymbol{\rm B}(\boldsymbol{y}^{(0)})\dfrac{{\rm d}\boldsymbol{y}^{(0)}}{{\rm d}t_{1}}\left(\dfrac{{\rm d}\boldsymbol{y}^{(0)}}{{\rm d}t_{1}}\right)^{\rm T}~{}\right]\right\}\mathcal{S}_{\alpha}(\boldsymbol{y}^{(0)})\end{split}$$

(42)

where $\mathcal{S}_{\alpha}~{}(\alpha=1,\cdots,2N-2)$ are

$$\mathcal{S}_{\alpha}(\boldsymbol{y})=\dfrac{{\rm Tr}\left[\boldsymbol{\rm P}(\boldsymbol{y})^{\rm T}\dfrac{\partial\boldsymbol{\rm B}}{\partial y_{\alpha}}(\boldsymbol{y})\boldsymbol{\rm P}(\boldsymbol{y})\boldsymbol{\rm\Lambda}(\boldsymbol{y})\boldsymbol{\rm N}\right]}{{\rm Tr}\left[\boldsymbol{\rm P}(\boldsymbol{y})^{\rm T}\boldsymbol{\rm K}\boldsymbol{\rm P}(\boldsymbol{y})\boldsymbol{\rm N}\right]}.$$

(43)

We can show that

$$\mathcal{S}_{1}=\cdots=\mathcal{S}_{N-1}=0.$$

(44)

See Appendix C. We remark that in the inside traces of $\mathcal{S}_{\alpha}$ the size of the matrices can be reduced from $2N-2$ to $N-1$ as derived in Appendix C.

We have eliminated the nontrivial unknown variables $w_{i}(t_{1})~{}(i=1,\cdots,N-1)$ from the averaged term (42), which depends on only the variables $\boldsymbol{\phi}^{(0)}(t_{1})$ and the constants $\boldsymbol{l}_{\ast}$, $\boldsymbol{\rm N}$, and $E^{(2)}$. Therefore, the equations of motion (26) for $\boldsymbol{\phi}^{(0)}$ is now closed, and dynamics depends on the excited normal modes $\boldsymbol{\rm N}$ and energy $E^{(2)}$.

Substituting Eq. (42) into Eq. (26), we have

$$\begin{split}&B^{\alpha\beta}(\boldsymbol{y}^{(0)})\dfrac{{\rm d}{}^{2}y_{\beta}^{(0)}}{{\rm d}t_{1}^{2}}+\left[D_{\alpha}^{\beta\gamma}(\boldsymbol{y}^{(0)})+\dfrac{1}{4}B^{\beta\gamma}(\boldsymbol{y}^{(0)})\mathcal{S}_{\alpha}(\boldsymbol{y}^{(0)})\right]\dfrac{{\rm d}y_{\beta}^{(0)}}{{\rm d}t_{1}}\dfrac{{\rm d}y_{\gamma}^{(0)}}{{\rm d}t_{1}}+\dfrac{\partial U_{\rm bend}^{(2)}}{\partial y_{\alpha}}(\boldsymbol{y}^{(0)})-\dfrac{E^{(2)}-U_{\rm bend}^{(2)}(\boldsymbol{y}^{(0)})}{2}\mathcal{S}_{\alpha}(\boldsymbol{y}^{(0)})=0\end{split}$$

(45)

for $\alpha=1,\cdots,2N-2$. The final result for the bending variables $\boldsymbol{\phi}^{(0)}(t_{1})$ is

$$\dfrac{{\rm d}{}^{2}\phi_{i}^{(0)}}{{\rm d}t_{1}^{2}}+F_{i}^{jn}(\boldsymbol{y}^{(0)})\dfrac{{\rm d}\phi_{j}^{(0)}}{{\rm d}t_{1}}\dfrac{{\rm d}\phi_{n}^{(0)}}{{\rm d}t_{1}}+G_{i}(\boldsymbol{y}^{(0)})=0$$

(46)

for $i=1,\cdots,N-1$. The functions $F_{i}^{jn}$ and $G_{i}$ are

$$F_{i}^{jn}(\boldsymbol{y})=\left[B_{\phi\phi}^{-1}\right]^{is}\left(\dfrac{\partial B_{\phi\phi}^{sj}}{\partial\phi_{n}}-\dfrac{1}{2}\dfrac{\partial B_{\phi\phi}^{jn}}{\partial\phi_{s}}+\dfrac{1}{4}\mathcal{T}_{s}B_{\phi\phi}^{jn}\right)$$

(47)

and

$$G_{i}(\boldsymbol{y})=\left[B_{\phi\phi}^{-1}\right]^{is}\left(\dfrac{\partial U_{\rm bend}^{(2)}}{\partial\phi_{s}}-\dfrac{E^{(2)}-U_{\rm bend}^{(2)}}{2}\mathcal{T}_{s}\right).$$

(48)

Here, we decomposed the matrix $\boldsymbol{\rm B}\in{\rm Mat}(2N-2)$ into four square submatrices of the size $N-1$ as

$$\boldsymbol{\rm B}(\boldsymbol{y})=\begin{pmatrix}\boldsymbol{\rm B}_{ll}(\boldsymbol{\phi})&\boldsymbol{\rm B}_{l\phi}(\boldsymbol{y})\\ \boldsymbol{\rm B}_{\phi l}(\boldsymbol{y})&\boldsymbol{\rm B}_{\phi\phi}(\boldsymbol{y})\\ \end{pmatrix}.$$

(49)

See Eq. (86) for explicit forms of the submatrices. The functions

$$\begin{split}&\mathcal{T}_{i}(\boldsymbol{y})=\mathcal{S}_{i+N-1}(\boldsymbol{y})=\dfrac{{\rm Tr}\left[\boldsymbol{\rm P}^{\rm T}\dfrac{\partial\boldsymbol{\rm B}}{\partial\phi_{i}}\boldsymbol{\rm P}\boldsymbol{\rm\Lambda}\boldsymbol{\rm N}\right]}{{\rm Tr}\left[\boldsymbol{\rm P}^{\rm T}\boldsymbol{\rm K}\boldsymbol{\rm P}\boldsymbol{\rm N}\right]}\\ &(i=1,\cdots,N-1)\end{split}$$

(50)

are introduced to renumber $\mathcal{S}_{\alpha}$ for avoiding zerocomponents shown in Eq. (44). It is worth noting that the spring potential $V_{\rm spring}$ is included in the final equations of motion (46) up to the second order, namely only through the matrix $\boldsymbol{\rm K}_{l}$.

We introduce the one-dimensional conformations whose set is denoted by

$$\mathcal{C}^{1}=\left\{(\phi_{1},\cdots,\phi_{N-2})~{}|~{}\phi_{i}\in\{0,\pi\}~{}(i=1,\cdots,N-2)\right\}.$$

(57)

A conformation in $\mathcal{C}^{1}$ is stationary as proven in Appendix D. Appearance of the bending potential forbids $\phi_{i}=\pi$ in general to avoid collision between beads, but we accept $\phi_{i}=\pi$ in this section to discuss the simplest case. Later we will perform numerical simulations under appearance of a bending potential which forbids $\phi_{i}=\pi$.

A conformation in $\mathcal{C}^{1}$ is symbolized by a sequence of $1$ and $-1$: $1$ represents the straight joint ($\phi=0$), and $-1$ the fully bent joint ($\phi=\pi$). The conformation symbol is denoted by $\boldsymbol{c}=(c_{1},\cdots,c_{N-2})$. We identify two symmetric conformations like $(1,-1,-1)$ and $(-1,-1,1)$, because each of which is mapped to the other by changing the starting end of the chain. All the possible one-dimensional conformations for $N=5$ are illustrated in Fig. 3 with their conformation symbols.

Let the point $\boldsymbol{\Phi}_{\ast}$ be stationary. Stability of this stationary point is obtained from the eigenvalues of the Jacobian matrix for Eq. (51),

$$\boldsymbol{\rm J}(\boldsymbol{\Phi}_{\ast})=\begin{pmatrix}\boldsymbol{\rm O}&\boldsymbol{\rm E}\\ -D\boldsymbol{G}(\boldsymbol{y}_{\ast}^{(0)})&\boldsymbol{\rm O}\\ \end{pmatrix},$$

(58)

where

$$D\boldsymbol{G}=\begin{pmatrix}\dfrac{\partial G_{1}}{\partial\phi_{1}}&\cdots&\dfrac{\partial G_{1}}{\partial\phi_{N-1}}\\ \vdots&\ddots&\vdots\\ \dfrac{\partial G_{N-1}}{\partial\phi_{1}}&\cdots&\dfrac{\partial G_{N-1}}{\partial\phi_{N-1}}\\ \end{pmatrix}$$

(59)

is the Jacobian matrix of the vector field $\boldsymbol{G}(\boldsymbol{y})$. We do not take the derivatives with respect to the variables $\boldsymbol{l}$, because Eq. (51) includes only the constant lengths $\boldsymbol{l}_{\ast}$.

Let us assume that the matrix $D\boldsymbol{G}(\boldsymbol{y}_{\ast}^{(0)})$ is diagonalizable and has the eigenvalues $g_{1},\cdots,g_{N-1}$. One eigenvalue, denoted by $g_{N-1}$, should be zero from the rotational symmetry, and we remove it from the stability criterion. The nontrivial eigenvalues of $\boldsymbol{\rm J}(\boldsymbol{\Phi}_{\ast})$ are $\pm\sqrt{-g_{1}},\cdots,\pm\sqrt{-g_{N-2}}$. See Appendix E for a proof and Fig. 4 for a schematic explanation of the relation between the eigenvalues of $D\boldsymbol{G}$ and $\boldsymbol{\rm J}$.

An eigenvalue $g_{i}$ is called a stable eigenvalue if $g_{i}\in(0,\infty)$, and a zeroeigenvalue if $g_{i}=0$, and an unstable eigenvalue if $g_{i}\not\in[0,\infty)$. The conformation represented by $\boldsymbol{\phi}^{(0)}_{\ast}$ is (neutrally) stable holm-etal-85 if and only if there is no unstable eigenvalue.

We first excite only one normal mode of the springs: All the diagonal elements of $\boldsymbol{\rm N}$ are zero except for one element. Stability of the one-dimensional conformations is summarized in Table 1 with dependency on the excited normal mode, where the normal modes are numbered in the ascending order of the eigenfrequencies around the considering conformation, as mentioned after Eq. (35). Stability is symbolized by S, Z, and U, and the number after S (Z, U) represents the number of the stable (zero, unstable) eigenvalues of $D\boldsymbol{G}$, whose sum is $N-2$. The symbol is omitted if the number of corresponding eigenvalues is zero. For instance, the symbol S2U1 represents that the conformation has $2$ stable eigenvalues and $1$ unstable eigenvalue, and the conformation is unstable.