\section{Introduction}

Spin torque nano oscillators (STNOs) have attracted considerable attention as next-generation multifunctional spintronic devices \cite{Locatelli2013,Ebels2017} over the last decade.
Potential applications reach from high data transfer rate hard disk reading \cite{Sato2012}  and wide-band high-frequency communication\cite{Choi2014,Ruiz-Calaforra2017,Kreissig2017,Litvinenko2019}
to broadband microwave energy harvesting \cite{Fang2019} or frequency detection \cite{Jenkins2016,Louis2017}. Furthermore, their exploitation in the reconstruction of bio-inspired networks for neuromorphic computing \cite{Torrejon2017,Romera2018} creates new opportunities for novel unconventional computing architectures.
After more than a decade of research, among different realizations of STNOs, vortex based STNOs (STVOs) have been considered a model system for nonlinear dynamics\cite{Pribiag2007,Dussaux2010} exhibiting the best performances in terms of coherence and tunability.
High output emission powers are simultaneously reached by the large magnetoresistive ratio of magnetic tunnel junctions (MTJs).
While different dynamic characteristics of the non-uniform vortex magnetization distribution have been extensively studied -- including specifically higher order modes corresponding to radial or orthoradial spin wave modes in the vortex tail \cite{Novosad2002,Guslienko2005,Ivanov2005,Taurel2016} --, the exploitation in STVOs is most commonly limited to the vortex core gyrotropic (G) motion \cite{Thiele1973}.
Indeed, only few studies report physics invoking the motion of the vortex core beyond this fundamental mode in STVOs.
For instance, Strachan \textit{et al. }  \cite{Strachan2008} report on the current-induced switching process of a uniform magnetization direction mediated by a vortex C- and G-state.
In Refs. \cite{Jin2009,Aranda2010}, the switching of a vortex core's polarity through a transient C-state mode is described, however, at pillar diameters of $d \lesssim 120\,$nm (depending on the material's exchange length) and an out-of-plane magnetized polarizer \cite{Jin2009,Aranda2010,Kawada2014}.
Under the same geometric limitations, an excited S-state \cite{Aranda2010} and C-state \cite{Kawada2014} of weakly \cite{Aranda2010} or significantly  \cite{Kawada2014}  increased frequency compared to the G-state is found in Refs. \cite{Aranda2010,Kawada2014}.
In this work, we demonstrate the characteristics of a spin transfer induced dynamic C-state with slightly lower frequency than the G-state, but likewise with $df/dI> 0$ and large coherence and emission power in a MTJ based STNO with an in-plane polarizer.
Both in experiment and micromagnetic simulation, we investigate  the stochastic nature of the transition between the two modes, facilitating potential random features\cite{Jenkins2019} that could be also exploited in the emerging field of neuromorphic spintronics.

\section{Experiments}
\label{sec:measurements}

The studied vortex STNO devices are circular shaped MTJs of $300\,$nm dot diameter, in-plane polarized reference layer and a $7\,$nm NiFe free layer (see supplementary\cite{supplementary_C-state_2020}).

\begin{figure}[tbh!]
  \centering
  \captionsetup[subfigure]{
  \subfloat[\label{fig:merged_simulation-exp}]
 \includegraphics[width=0.62\columnwidth]{graph-bilder/spectra_experiment_520mT_4.pdf
}
\\%\hfill
\subfloat[\label{fig:basic_properties}]

\includegraphics[width=0.66\textwidth]{graph-bilder/basic_properties_1_with_i_3.pdf}

  \caption[]{Measurement of the dynamic states.  (a) Power emission spectra  and (b) corresponding extracted properties of the oscillation vs. $I_{dc}$: (b.i) emission power, (b.ii) linewidth FWHM, (b.iii) nonlinear parameter $\nu$, and (b.iv) average sample resistance $R_0$.  An out-of-plane  magnetic field of $520\,$mT was applied.  }
  \label{fig:exp-data_1}
\end{figure}

In fig. \ref{fig:exp-data_1}, we show typical measurements presenting different dynamic states:
In fig. \ref{fig:merged_simulation-exp}, the power emission spectra of the current supplied STVO are plotted.
For $I_{dc} \in [7;9.7]\,$mA, a quasi-linear evolution of the oscillation frequency with the current is observed.
As shown in fig. \ref{fig:basic_properties},  the oscillation parameters in this interval evolve as expected for the gyrotropic vortex motion: The linewidth decreases with the injected dc current $I_{dc}$  and the oscillation amplitude increases, as also the vortex core position radius on its stable limit cycle grows with $I_{dc}$.
In the following, we refer to these characteristics in the mentioned interval as the gyrotropic (G) state.

For $I_{dc} \in [9.7;13.5]\,$mA, the frequency decreases rapidly with strongly enhanced linewidth (fig. \ref{fig:basic_properties}.ii) and nonlinearity parameter\cite{Slavin2009,Wittrock2019_PRB} $\nu$ (fig. \ref{fig:basic_properties}.iii)   compared to the low-linewidth G-state.
This regime represents a transition between two  frequency values, whereas the frequency at the higher current value is lower.
It is generally observed and the transition can be more (see supplementary) or less (fig. \ref{fig:exp-data_1}) discrete at different values of $I_{dc}$ depending on the chosen experimental parameters.
Interestingly, the average sample resistance, which usually decreases monotonically with $I_{dc}$ due to the TMR bias dependance, exhibits a relative increase within the transition (fig. \ref{fig:basic_properties}.iv) indicating a more antiparallel magnetization configuration.
A third dynamical regime appears for $I_{dc} \gtrsim 13.5\,$mA, in which another quasi-linear frequency evolution with the current is observed.
The linewidth FWHM stabilizes after the transition down to  $\sim 300\,$kHz, a similar value as obtained as well in the G-state.
Importantly, the emission power is much higher than for the G-state oscillation (factor $3$-$4$ in the presented measurement) and does furthermore  not systematically change with the applied dc current.
Note that however, the power is subject to rather large fluctuations at constant linewidth, potentially caused by an enhanced flicker noise in this regime\cite{Wittrock2019_PRB,Wittrock2020-SciRep}.
For very high current values above $17\,$mA, the linewidth increases again.
As for the nonlinearity parameter $\nu$, it first reduces almost to zero after the transition and consequently, is practically irrelevant in the measured C-state regime. However, measurements indicate that it significantly increases again for even higher currents in the C-state.

 \begin{figure}[!hbt]
 	\centering
 	{\includegraphics[width=0.658\textwidth]{graph-bilder/TimeTraces_4.pdf}  }
 	\caption{Measured heterodyne voltage signal vs. time at different current values in the (a) G-state, (b-d) transition regime and (e) C-state. The given frequencies (blue) are average values over the corresponding marked (b-d) or total (a,e) sequence.  }
 	\label{fig:time-traces}
 \end{figure}

In fig. \ref{fig:time-traces}, we record the oscillation time signal using a heterodyne detection technique (signal down-conversion to $5$-$10\,$MHz via high-side injection and low-pass filtering) in order to resolve  small variations of frequency.
Note that then, one period of the heterodyne signal roughly corresponds to $\sim 50$ periods of the pure signal.
At the beginning of the transition regime at $9.94\,$mA  (fig. \ref{fig:time-traces}.b),  the oscillator is mainly in the  G mode ($336.5\,$MHz in the intervall $t\sim [8.5;9.3]\,\upmu$s).
In addition, for example at $t\approx 8.5\,\upmu$s, the signal shows a kink (orange arrow), which is interpreted as a part of the oscillation in the C-state or potentially pinning at the disk boundary for a short moment of time.
  For a short time in the interval $\sim [9.3;9.7]\,\upmu$s, the system oscillates at a slightly different frequency:
  We assume here, that the oscillator sometimes performs even several oscillation periods in the C-state before switching back to the G-state (see also the supplementary material\cite{supplementary_C-state_2020}).
For a larger dc current $I_{dc} = 10.8\,$mA in the center of the transition regime (fig. \ref{fig:time-traces}.c), the time signal is much less coherent, which we interpret as a  regular switching between the states together with potential pinning events.
For $I_{dc}=12.8\,$mA (fig. \ref{fig:time-traces}.d), at the end of the transition, the vortex core is again only rarely constrained at the potential barrier.
For currents well above or below the transition region, and hence for the G- or C-state (fig. \ref{fig:time-traces}.a and \ref{fig:time-traces}.e, resp.), the time traces show a stable oscillation signal without observed instabilities and state switching events.

\section{Micromagnetic simulation}
\label{sec:simulation}

\begin{figure*}[bth!]
\centering
\floatsetup[subfigure]{}
\ffigbox[\textwidth]
\begin{subfloatrow}[2]
\thisfloatsetup{capbesidesep=none,objectset=centering,%
capbesideframe=yes,capbesideposition={right,bottom}}
\ffigbox[\FBwidth]
\includegraphics[height=0.344\textwidth]{simulation_graphs/T=300_FLT=0,4_Frequ_radius_7.pdf}}{ \caption{ }  \label{fig_C-state:sim_T=300_FLT=0,4_Frequ_radius}}
\thisfloatsetup{capbesidesep=none,objectset=centering,%
capbesideframe=yes,capbesideposition={right,bottom}}
\ffigbox[\FBwidth]
\includegraphics[height=0.352\textwidth]{simulation_graphs/T=300_FLT=0,4_TimeTraces_3_2.pdf}}{ \caption{ } \label{fig_C-state:sim_T=300_FLT=0,4_TimeTraces}}
\end{subfloatrow}

\vspace{2ex}

   \thisfloatsetup{capbesidesep=columnsep,objectset=centering,%
capbesideframe=yes,capbesideposition={left,center}}
    \fcapside[0.5\textwidth]
     \begin{subfloatrow}[2]
      \ffigbox[\FBwidth]{%

    \includegraphics[width=0.92\columnwidth]{simulation_graphs/T=300_FLT=0,4_T=0_TimeTrace_transition_1.pdf} }{\caption{}
    \label{fig:T=0_with_FLT}  }
     \end{subfloatrow}
     }
          { \caption[Parameter space of the PLL+STNO phase noise PSD ]{Micromagnetic simulation of the magnetization dynamics
  of gyrotropic (G) and C-state and their transition. {(a)} Frequency  and vortex core radius evolution with applied dc current. {(b)} Planar mean  magnetization vs. time at different current values.  (c) Discrete transition when slightly increasing $I_{dc}$ with time at $T=0\,$K.
}
 \label{fig_C-state:sim_T=300_FLT=0,4} }
    }
 {}
\end{figure*}

In order to further investigate the nature of the different dynamic modes, we perform mumax$^3$ micromagnetic simulations \cite{mumax} of the magnetization dynamics  (see supplementary\cite{supplementary_C-state_2020}).
For this purpose, we include the different spin-torque contributions that can control the current driven magnetization dynamics, namely the Oersted field,  Slonczewski (STT) and the field-like torque (FLT)\cite{Ralph}.

Simulation results are shown in fig. \ref{fig_C-state:sim_T=300_FLT=0,4} and   reproduce the frequency characteristics of the measurements as described for the experimental data:
 Up to  $I_{dc} \sim 13\,$mA, we find  in fig. \ref{fig_C-state:sim_T=300_FLT=0,4_Frequ_radius} the gyrotropic mode, i.e. the quasi-circular motion of the vortex core takes place inside the nanodisk as depicted in fig. \ref{fig_C-state:sim_T=300_FLT=0,4_Frequ_radius}'s inset (see also videos in the supplementary\cite{supplementary_C-state_2020}).
 Note that the vortex gyration is not symmetric relative to the disk center due to the FLT\cite{Dussaux2012}, whose efficiency compared to the STT is set to $\xi_{fl} = 0.4$, a reasonable value for MTJ STVOs\cite{Chanthbouala2011}.

For $I_{dc} \gtrsim 13\,$mA, we observe a transition to the C-state which stabilizes for higher dc currents.
At this point, the introduced expression \textit{"C-state"}\ also becomes clear: The simulations reveal that at higher currents $I_{dc} \gtrsim 17\,$mA, the vortex core is not continuously found to be inside the  nanodisk.
Instead,  it establishes that for a part of the oscillation period, the vortex core is inside and for the other part, it is imaginarily outside the disk, i.e. in the latter the in-plane magnetization of the vortex tail forms a \textit{C} which precesses following an imaginary vortex core.
The situation is depicted in fig. \ref{fig_C-state:sim_T=300_FLT=0,4_Frequ_radius}'s inset.
Furthermore, the lower panel of fig. \ref{fig_C-state:sim_T=300_FLT=0,4_Frequ_radius} clearly proves  the vortex core  oscillation radius partly inside, partly outside the nanodisk of radius $R$.
In the supplementary material\cite{supplementary_C-state_2020}, movies of the occurring dynamic situations can be found.
We also perform simulations without FLT, i.e.  $\xi_{fl} = 0$  (see supplementary material\cite{supplementary_C-state_2020}).
The simulation results in this case are fundamentally different and describe after the transition from the G-state  a pure C-state (no renucleation of the vortex core) with significantly higher frequency in the order of $1\,$GHz.

 These results qualitatively reproduce what has been obtained by Kawada \textit{et al. }\cite{Kawada2014} for the experimental case of a GMR based STVO, where indeed the  FLT is negligible.
Thus, the FLT in MTJs is an important parameter in order to control the described C-state characteristics.

In the transition between the two modes in fig. \ref{fig_C-state:sim_T=300_FLT=0,4_Frequ_radius}, different situations can occur:
Firstly,  the system incoherently switches between both states, performing several periods in the G- and several periods in the C-state, as also observed experimentally.
Secondly, the vortex core pins at the disk boundary for a certain period of time, hence stopping oscillator operation until it resumes dynamics.
In fig. \ref{fig_C-state:sim_T=300_FLT=0,4_TimeTraces}, exactly this behavior can be observed through the planar mean magnetization: For  $I_{dc}=9.86\,$mA, the system steadily oscillates in the G-state.
In this case, oscillations of the in-plane mean magnetization can be obtained because
of the asymmetry in the gyrotropic motion (not centered) induced by the FLT,  particularly included  in the simulations.
For $I_{dc}=12.24\,$mA, we observe a stochastic  switching between the states: For example, from $0.6< t < 0.72\,\upmu$s, the system is in the C-state.
It has a lower frequency and the planar magnetization, due to the lack of the vortex core, is slightly enhanced.
From $0.1< t < 0.2\,\upmu$s, it can be observed that the vortex core pins at the disk boundary resulting in a quasi-constant in-plane magnetization until readopting oscillations.
Comparing to experimental measurements (fig. \ref{fig:exp-data_1}), the switching between the two states translates into a broadening of the effective linewidth.
Note that there, the likelihood of vortex core pinning at the boundary and thus, the range of the transition depends on several experimental parameters, such as the applied magnetic field, sample dimensionality, free layer material, etc. resulting in different characteristics of the transition (shown in the supplementary\cite{supplementary_C-state_2020}).
For $I_{dc}=17.18\,$mA in fig. \ref{fig_C-state:sim_T=300_FLT=0,4_TimeTraces}, the system oscillates in the C-state providing coherent oscillations, although the vortex core pins at the boundary at some very short moments from time to time. The number of pinning events decreases at even higher currents.
Notably, figs. \ref{fig_C-state:sim_T=300_FLT=0,4_Frequ_radius} \& \ref{fig_C-state:sim_T=300_FLT=0,4_TimeTraces} show that the transition region between G- and C-state is completely stochastic. The state probability can be continuously tuned within the transition until, at higher current values, the C-state is stabilized, and at lower current values, the G-state.

In fig. \ref{fig:T=0_with_FLT}, we show the transition between G- and C-state when performing the simulation without thermal noise at $T=0\,$K.
When slightly increasing the current to the critical transition current at $\sim 24\,$ns, the transition from G- to C-state oscillations is observed.
Along with the reduced frequency in the C-state, also an asymmetry within one period (equally present for $T=300\,$K) is recognized: When the vortex is expelled (high mean in-plane magnetization) for approximately half a period, the frequency is slightly reduced compared to the other half-period when the vortex core is located inside the nanodisk.
Moreover, the transition between G- and C-state is discrete and both situations exhibit stable oscillations, indicating that the unstable transition region and stochastic state switching described beforehand  is noise induced due to thermal fluctuations.

\section{Discussion}

The emission power in the C-state is considerably higher than in the G-state (fig. \ref{fig:basic_properties}).
Comparing with the performed simulations, we attribute this to a  larger active magnetic volume contributing to the magnetoresistive rf signal.
The quasi-zero nonlinearity parameter $\nu$ in the C-state might be due to a changed confinement when the magnetization distribution is not vortex-like.
Note that $\nu$ might again become $\nu > 0$ at even higher applied current values.
Moreover, the nonlinearity parameter within the transition region is rather an effective value as the system is bistable and the single-mode hypothesis does not hold any more.
However, as the frequency tunability linked to the large nonlinear parameter $\nu$ is strongly enhanced in this regime, we anticipate that the locking ability to an external stimulus is significantly improved here.
Especially interesting from an application point of view is the vicinity of a high $\nu$ to a $\nu\approx 0$ regime: With the locking efficiency proportional to $\sqrt{1+\nu^2}$, the ability to respond and synchronize to an external signal could thus be easily tuned by only changing the applied dc current.

\section{Conclusion}

In this study, we observe and describe a dynamic C-state beyond the vortex gyrotropic (G) motion in experiment and, consistently, in micromagnetic simulations.
Induced by the applied dc current,
the observed G- to C-state transition is found to be facilitated by the field-like torque in MTJ based STVOs and notably completely stochastic in the presence of thermal noise.
The state probability can be continuously tuned within the transition regime until, at higher current values, the C-state is stabilized, and at lower current values, the G-state.
We believe that the stochastic feature might be exploited for cognitive computing,
for which STNOs have recently been implemented in numerous ways. Besides oscillations,
stochasticity is an important feature for energy efficient neuromorphic operations\cite{Grollier2016}.
The combination of
both, nonlinear oscillations and stochasticity, and furthermore the ability to synchronize or respond to external stimuli, in one unique spintronic device can open the path to
more multifunctional neuron-like building blocks necessary for neuromorphic computing\cite{Grollier2016}.

\section*{\footnotesize{Supplementary Material}}

See supplementary material\cite{supplementary_C-state_2020} for details on the samples and  a further experimental characterization of the occurring transition from G- to C-state. This includes the demonstration of
a more discrete transition and its characterization with the applied magnetic field. Furthermore, simulation results for $\text{FLT}=0$ are shown and videos of the different regimes included.

\section*{\footnotesize{Acknowledgments}}

S.W. acknowledges financial support from Labex FIRST-TF under contract number ANR-10-LABX-48-01.
M.R. acknowledges support from Spanish MICINN (PGC2018-099422-A-I00) and Comunidad de Madrid (Atracción de Talento Grant No. 2018-T1/IND-11935).
The work is supported by the French ANR project ''SPINNET'' ANR-18-CE24-0012.

\section*{\footnotesize{Data Availability Statement}}

The data that support the findings of this study are available from the corresponding author upon reasonable request.

\section{Measurement}

\subsection{Sample details}

The experimental measurements presented in this work have been performed on vortex based spin transfer oscillators (STVOs). This type of STNOs exploits the dynamics of a magnetic vortex in a circular shaped nanodisk and converts it into an electric rf signal due to the magnetoresistive effect\cite{Dussaux2010}.

The studied samples consist of a pinned layer made of a conventional synthetic antiferromagnetic stack (SAF), a MgO tunnel barrier and a NiFe free layer in a magnetic vortex configuration.
The  magnetoresistive ratio related to the tunnel magnetoresistance effect (TMR) lies around $50\,$\% at room temperature.
The SAF is composed of \selectlanguage{ngerman} PtMn($20$)/""Co$_{70}$Fe$_{30}$($2$)/""Ru($0.85$)/""Co$_{40}$Fe$_{40}$B$_{20}$($2.2$)/""Co$_{70}$Fe$_{30}$($0.5$) and the total layer stack is Ta(3)/""CuN(30)/""Ta(5)/""SAF/""MgO($1$)/""Co$_{40}$Fe$_{40}$B$_{20}$($1.5$)/""Ta($0.2$)/""Ni$_{80}$Fe$_{20}$(7)/""Ta(5)/Ru(7), with the nanometer layer thickness in brackets.
\selectlanguage{english}
The layers are deposited on high resistivity SiO$_2$ substrates by ultrahigh vacuum magnetron sputtering, subsequently annealed for $2\,$h at $T=330\,^{\circ}$C at an applied magnetic field of $1\,$T along the SAF's easy axis, and patterned using e-beam lithography and Ar ion etching.
The $0.2\,$nm Ta-layer decouples the growth of CoFeB and NiFe, utilizing the high tunnel magnetoresistance (TMR) ratio of the crystalline CoFeB-junction and the magnetically softer NiFe for the vortex dynamics.
The circular tunnel junctions have an actual diameter of $2R = 300\,$nm.

\subsection{Transition characteristics}

\begin{figure}[bth!]
  \centering
  \subfloat[  \label{fig_C-state:spectra_and_spectrum_685mT_discrete_transition-1}]
  {
  \includegraphics[width=0.415\textwidth]{graph-bilder/spectra_and_spectrum_685mT_discrete_transition_1.pdf}  }

  \hfill
 \subfloat[  \label{fig_C-state:spectra_experiment_RTN-2}]
  {
  \includegraphics[width=0.4\columnwidth]{graph-bilder/spectra_experiment_RTN_1.pdf}   }

  \caption[Frequency spectra and transition vs. $I_{dc}$ in different measurements]{Frequency spectra vs. $I_{dc}$ in different measurements (a) and (b) showing a G$\rightarrow$C-state transition.  In (a.i) the frequency spectra for different current values in or close to the transition, corresponding to the frequency spectra surface plot in (a.ii), are shown.
}
  \label{fig_C-state:other-spectra}
\end{figure}

The transition between the modes in fig. 2 of the main paper is relatively smooth.
In fig. \ref{fig_C-state:other-spectra}, we show measurements carried out on different spin torque vortex oscillator (STVO) devices under different experimental conditions (same wafer, but different STVO).
This second measurement series is shown because a more discrete transition between the two states (fig. \ref{fig_C-state:spectra_and_spectrum_685mT_discrete_transition-1}) or some stochastic switching (fig. \ref{fig_C-state:spectra_experiment_RTN-2}) is observed in the spectra, generally highlighting that the transition characteristics can be changed through the experimental parameters.
Indeed, the measured spectrum in the transition, shown in fig. \ref{fig_C-state:spectra_and_spectrum_685mT_discrete_transition-1}.i, even shows two distinct frequency peaks becoming more or less pronounced depending on the position in the transition regime.
In fig. \ref{fig_C-state:spectra_experiment_RTN-2}, the transition at $I_{dc} \in \sim [9;11.5]\,$mA is also rather discrete and the system furthermore  switches between the two states when sweeping $I_{dc}$ inside the transition.
The shown characteristics indicate that the transition from G- to C-state is stochastic.

\begin{figure}[bth!]
  \centering

 \subfloat[  \label{fig_C-state:TimeTraces-b}]
  {  \raisebox{0.18cm}{ \hspace{0.0cm}
  \includegraphics[width=0.4\textwidth]{graph-bilder/TimeTrace_4,26kG_9mA_2.pdf}   } }
   \subfloat[  \label{fig_C-state:TimeTraces-c}]
  {
  \includegraphics[width=0.545\textwidth]{graph-bilder/TimeTrace_4,26kG_11,8mA_3.pdf}   }

  \caption[Measured heterodyne voltage signal vs. time in the transition regime G$\rightarrow$C-state -- 2]{Measured heterodyne voltage signal vs. time in the transition regime at different current values, corresponding to the spectra shown in fig. \ref{fig_C-state:spectra_experiment_RTN-2}. The given frequencies (blue) are average values over the corresponding marked sequence.
}
  \label{fig_C-state:TimeTraces_different_conditions}
\end{figure}

In fig. \ref{fig_C-state:TimeTraces-b} \& \ref{fig_C-state:TimeTraces-c}, we show time traces recorded corresponding to the spectra in fig. \ref{fig_C-state:spectra_experiment_RTN-2}.
At lower currents in the transition interval of this measurement ($\sim [8.9;10.5]\,$mA), we observe voltage signals as shown in fig. \ref{fig_C-state:TimeTraces-b}.
 Here, the oscillation frequency is broader dispersed.
 Characteristically, we observe kinks in the signal as e.g. at $5$ or at $5.3\,\upmu$s and marked by green arrows in fig. \ref{fig_C-state:TimeTraces-b}, which is interpreted as a frequent state-change and/or pinning of the vortex core at the disk boundary.
At higher currents in the transition interval ($\sim [10.6;12.5]\,$mA), we find voltage time traces as shown in fig. \ref{fig_C-state:TimeTraces-c}.
A stable oscillation with a specific frequency over a long range of time ($\sim 100\,\upmu$s) is observed.
Beyond this time period, the system changes its state to another frequency of $\sim 4\,$MHz difference for a short duration.
In this interval, the system oscillates mainly in the C-state, but rarely changes to the G-state where it performs several periods.
Pinning at the potential barrier in this regime does not play a significant role and the state change is rather discrete.

\subsection{Dependence of the mode transition on the magnetic field}

In this section, we investigate the dependence of the mode transition on the applied out-of-plane magnetic field $H_{\perp}$.
Regarding fig. \ref{fig_C-state:transition_vs_H-a}, it is observed that an applied field $H_{\perp}$ lowers the energy barrier and thus facilitates the described transition between G- and C-state.
We attribute this fact to the increase of the strength of the driving spin torque force acting on the vortex core, which is determined by the out-of-plane component of the STVO's fixed layer which grows with  $H_{\perp}$  \cite{Dussaux2012}.
Furthermore, the applied magnetic field deforms the magnetic vortex distribution (see Ref. \cite{Dussaux2012}). Since the magnetization $z$-component at the boundary increases\cite{Dussaux2012}, the transition energy barrier is additionally lowered.
Note that the results shown in fig. \ref{fig_C-state:dependence_vs_H} are measured for a different sample than the ones shown beforehand, qualitatively showing the same characteristics but with different values.

\begin{figure}[bth!]
  \centering
  \subfloat[  \label{fig_C-state:transition_vs_H-a}]
  {
  \includegraphics[width=0.4425\textwidth]{graph-bilder/transition_vs_H_2.pdf}  }

  \hfill
 \subfloat[  \label{fig_C-state:f_vs_H-b}]
  {
  \includegraphics[width=0.436\columnwidth]{graph-bilder/f_vs_H_Miguel_new_1.pdf}   }

  \caption[State transition dependence on the magnetic field]{(a) Evolution of the transition range between G- and C-state with the  applied magnetic field $\mu_0 H_{\perp}$. (b) Frequency vs. $\mu_0 H_{\perp}$ at constant dc current values (different sample).
  }
  \label{fig_C-state:dependence_vs_H}
\end{figure}

In fig. \ref{fig_C-state:f_vs_H-b}, we plot the frequency evolution vs. $\mu_0 H_{\perp}$ and observe that the different modes can also be triggered by the magnetic field.
The transition between the states occurs for lower fields if $I_{dc}$ is higher, in agreement with the above interpretation of a higher spin torque strength facilitating the transition.
Note that the slope $df/dH_{\perp}$ in the C- is larger than in the G-state, which might be due to a changed influence of the current induced Oersted field on the dynamics.

\section{Simulation}

The micromagnetic simulations are performed with mumax3 \cite{mumax}.
The dimension of the STVO's free layer (nanodisk diameter $300\,$nm, $7\,$nm thick NiFe free layer) is discretized into $86\times 86\times 1$ cells. A NiFe saturation magnetization of $8\cdot 10^5\,$A/m, exchange stiffness $A_{ex}=1.3\cdot 10^{-11}\,$J/m, and damping $\alpha=0.01$ are chosen.
In order to solve the Landau-Lifshitz-Gilbert-Slonczewski equation of magnetization dynamics, a standard Runge-Kutta 4/5 ODE solver has been used.
Temperature is either set to $T=300\,$K or $T=0\,$K.

  \begin{wrapfigure}[17]{r}{0.467\columnwidth}
\centering
\vspace{-0.68cm}
\includegraphics[width=0.906\columnwidth]{simulation_graphs/Simul_Frequ_T=0_FLT=0_2.pdf}
\caption[Micromagnetic simulation of G- and C-state without field-like torque and thermal noise]{Micromagnetic simulation without field-like torque and thermal noise: Frequency evolution and magnetic states as a function of the applied current $I_{dc}$.
 }
\label{fig_C-state:sim_T=0_FLT=0}
\end{wrapfigure}

If the field-like torque is set to zero, the vortex motion does not exhibit an asymmetric component any more.
The situation is shown in fig. \ref{fig_C-state:sim_T=0_FLT=0} for $T=0\,$K: In a small range at lower currents $I_{dc} \in \sim  [9.9; 10.2]\,$mA, above the current threshold for auto-oscillations, the gyrotropic motion is observed until at $I_{dc} \approx 10.2\,$mA, the vortex core is expulsed.
The free layer is then quasi uniformly magnetized between $I_{dc} \approx 10.2\,$mA and $I_{dc} \approx 14\,$mA without considerable dynamics. At $I_{dc} \approx 14\,$mA, the C-state oscillations start, as shown in fig. \ref{fig_C-state:sim_T=0_FLT=0}.
Note that here, if $T\neq 0\,$K in the presence of thermal noise, the transition between the G- and the C-state would be more continuous  with a smaller range without oscillations.
However, the C-state in fig. \ref{fig_C-state:sim_T=0_FLT=0} has a significantly increased frequency of $\sim 900\,$MHz compared to the G-state and can be seen as a \textit{pure} C-state where no renucleation of the vortex core inside the nanodisk takes place.
This observation is in agreement with findings in Ref. \cite{Kawada2014}, where similar measurements were conducted in GMR-based STVOs (i.e. negligible field-like torque) and a highly frequency increased C-state dynamic mode in the GHz range has been observed.
The results highlight the role of the field-like torque, present in TMR based STVOs, to induce vortex C-state dynamics beyond the typical gyrotropic motion.


\section{Introduction}

Molecular motors that work out-of-equilibrium are key elements for the viability of cells. Typical examples include myosin and kinesin motors~\cite{Sanchez2012nature}, RNA and DNA polymerases~\cite{Alberts2014}, helicases~\cite{Calladine1997} and condensin~\cite{Le2013,Terakawa2017}. Although many of these protein complexes are biochemically well characterised, we are far from having a full understanding of their collective action \emph{in vivo}. In some cases, the behaviour may be tuned by inter-motor interactions that are mediated by the substrate itself, a germane example being the change in polymerase-DNA binding affinity due to substrate supercoiling~\cite{Brackley2016supercoil,Naughton2013}.

Most previous theoretical and computational work on non-equilibrium forces within live cells has focused on modelling the propagation of stresses on biopolymer networks~\cite{Lenz2014,Broedersz2014} or active fields which mimic cytoskeleton dynamics~\cite{Tjhung2012}.
While molecular motors can now be readily investigated using single-molecule techniques~\cite{Herbert2008,michaelis2013single,Ganji2018}, the emergent properties of systems where a collection of motors interact on, and through, their target polymer substrate remain poorly studied~\cite{bianco_loca,SaintillanMotor}.

Here we propose a coarse-grained computational model to describe the action of generic molecular motors on polymers with the aim of improving our understanding of the effects of ATP-consuming translocating machineries on the conformational and dynamic properties of their target substrates. Within our model, motors undergo free diffusion in 3-D {until they encounter the polymer, at which point they become ``bound'' to it via an attractive interaction.} {While bound, a motor experiences an active force propelling it along the polymer; at the same time the polymer is subjected to an equal but opposite active force. The result is that the motor displays relative motion with respect to the substrate. The motor can unbind from the polymer, for example due to thermal motion or because it reaches the polymer's end.}
While our model is generic, and our aim is to understand its underlying physics, it could be viewed as a model for DNA or RNA polymerases moving along DNA or kinesin molecules stepping on microtubules~\cite{Alberts2014}.

{There has been a great deal of recent simulation work exploring active polymers, and our model differs from these in a number of ways. To our knowledge our work is the first to explicitly treat discrete motors that track along a polymer in a fully 3-D system. A number of recent studies~\cite{Harder2014,Jiang2014,Isele-Holder2015,Eisenstecken2016,bianco_loca,Anand2018,SaintillanMotor} have considered a chain of beads, where each bead experiences an active force, either directed along the chain bonds~\cite{Jiang2014,Isele-Holder2015,Anand2018} or in a direction independent of the chain~\cite{Eisenstecken2016,bianco_loca}. These works have mainly focused on regimes of different behaviour at different levels of activity (e.g. transitions from active translation to rotational motion~\cite{Jiang2014} or to spiral configurations~\cite{Isele-Holder2015} have been observed -- for a review see Ref.~\cite{Winkler2017}). In the present work the polymer only experiences an active force at the points where the explicit motors are bound (and these vary in time), and we focus on the effect on the conformational and dynamical properties of the polymer. Another class of models have more specifically focused on motor proteins found in the cytoskeleton, like myosin and kinesin which have provided a model experimental system for active materials~\cite{Schaller2010,Sanchez2012}. These can be modelled at a field theoretic level, but there have been a number of simulation studies which resolved the motors explicitly~\cite{nedelec2007,Freedman2017,Ravichandran2017,Gupta2019}. Unlike the current work, those studies focused on dense suspensions of motors and their substrates, and/or 2-D systems where the motors are fixed to a surface, and have often considered joined pairs of motors which can crosslink different substrate molecules to drive them into relative motion. In contrast, here we consider a single, semi-flexible polymer which is free to move in 3-D, and is subject to the action of multiple motors. }

Our simulations reveal that motor activity affects the steady state conformations of the polymer substrate; we find that these are more crumpled
with respect to their passive counterparts. We attribute this effective softening of the substrate to the emergence of double-folded branches, or hairpins; the effect can be tuned by changing, e.g., the strength of motor activity. We also find that, for large enough motor activity, the centre of mass of the polymer substrate can display super-diffusive behaviour at short times, with a return to diffusive behaviour at long times. Again this effect can be controlled by varying the activity or number of motors.
We conclude our work by investigating the effect of motors on
the propensity of the polymer to form knots.
This is relevant for many biological processes involving DNA, where knotting would be detrimental (e.g knots and tangles might hinder transcription and replication, repair of double strand breaks, or segregation of chromosomes).
We find that the motor activity effectively \textit{reduces} the steady state probability of finding knotted conformations. This is surprising since at equilibrium, in a good solvent, and in the absence of confinement, softer substrates show an enhanced probability of knotting~\cite{Coronel2017}. Thus, our findings point to an intriguing non-equilibrium effect of molecular motors that could be important in many biological systems.

\begin{figure}[t!]
	\centering
	\includegraphics[width=.45\textwidth]{./Figures/Model.png}

	\caption{Schematic of our model for molecular motors acting on a generic polymer substrate. Diffusing motors (blue) can interact with the polymer backbone (grey); when bound, the motors (red) and the polymer experience opposite forces which cause relative motion. On the right hand side of the figure we show the polymer in more detail: the vector $\mathbf{\hat{u}}$ defines an orientation for each bead and $\phi_i$ is then defined as the angle between $\mathbf{\hat{u}}_i$ and the vector joining beads $i$ and $i+1$; an orientation potential acts to minimise $\phi_i$ (see text and Appendix~\ref{model_details} for details).}
	\label{fig:model}
\end{figure}

\section{The Model}
We perform Langevin dynamics simulations of a semi-flexible bead-and-spring polymer consisting of $L$ beads in solution with $N$ motors, as described schematically in Fig.~\ref{fig:model}. We use a standard polymer model which includes finitely-extensible non-linear elastic (FENE) springs connecting consecutive beads, a Kratky-Porod potential for triplets of beads providing bending rigidity, and a Weeks-Chandler-Anderson (WCA) potential to regulate steric interactions (see Appendix \ref{model_details} for full details). The intrinsic stiffness of the polymer is controlled by the bending energy $K_{\rm BEND}$ which appears as a parameter in the Kratky-Porod potential. Each polymer bead also displays an orientation, which is tracked by a vector $\mathbf{\hat{u}}_i$ as shown in Fig.~\ref{fig:model}, and {we include} a potential to orient this vector along the local direction of the polymer backbone {(if the position of bead $i$ is $\mathbf{r}_i$ and the orientation is $\hat{\mathbf{u}}_i$, then the potential acts to minimise the angle $\phi_i$ between $\hat{\mathbf{u}}_i$ and $|\mathbf{r}_{i+1}-\mathbf{r}_i|$)}.

The motors are represented by freely diffusing beads which interact with each other sterically (via a WCA potential) and via a Morse-like attractive {potential with the polymer beads. The latter is given by}
\begin{align}{
\label{motorattraction}
U_{\rm MOT}(r) = \left\{
\begin{array}{ll}
K_{\rm MOT} \left[ (e^{-4r} - 2e^{-2r}) - \right.&  \\
\left. ~~~~ (e^{-4r_{\rm m\,cut}} - 2e^{-2r_{\rm m\,cut}}) \right], & r\leq r_{\rm m\,cut} \\
0, & \mbox{otherwise},
\end{array} \right.}
\end{align}
{Here $r$ is the the separation of the motor and the polymer bead, and $K_{\rm MOT}$ is the interaction strength.} {The potential} has a minimum when the motor and polymer beads overlap. {We consider a motor to be ``bound'' to a polymer bead when the centre-to-centre separation of the bead and the motor is less than $r_{\rm m\,cut}$, which we set to be the same as the motor diameter.}
In addition to this attraction, when a motor is bound to a polymer bead, it experiences an active force of magnitude $f$ {in the direction of the} orientation vector $\hat{\bold{u}}_i$ of the bead{; the polymer bead experiences an opposite active force $-f\hat{\mathbf{u}}_i$ (in compliance with Newton's third law).} {It is possible for a motor to be within a distance $r_{\rm m\, cut}$ of more than one polymer bead at the same time: in this case each interacting polymer bead will experience an active force, and the total active force on the motor will be the sum of that due to each interacting polymer bead. (In practice motors spend over 95\% of their ``bound'' time interacting with two consecutive polymer beads.)}

The result {of the active force is that the motors and the polymer are driven into relative motion; due to the connectivity of the chain, the polymer tends to move less compared to motors, so the latter effectively ``track'' along the polymer substrate. Motors can become unbound from the polymer due to thermal fluctuations (controlled by the parameter $K_{\rm MOT}$), and also active effects (e.g. tracking off the end of the polymer -- this is discussed in mote detail below).}
Since {each monomer has a preferred orientation, i.e.} the polymer is polar, motors {will} move in the parallel or anti-parallel direction depending on the sign of $f$. In the present work we consider only the case where all motors are the same and $f\geq0$.

For simplicity motors and polymer beads are taken to have the same {diameter}, $\sigma$, mass, $m$, and unless otherwise stated they experience the same friction $\xi$ due to the implicit solvent {(though in reality a motor such as RNA polymerase would be much bigger)}.
{We use the LAMMPs molecular dynamics software~\cite{LAMMPS} to perform the simulations: a Langevin equation for each polymer bead and motor is solved using a velocity-Verlet algorithm. A separate Langevin equation is solved for each monomer orientation $\hat{\mathbf{u}}_i$. Full details are given in Appendix \ref{model_details}. Throughout we quote energies and times in units of $k_BT$ and $\tau_{\rm LJ}=\sqrt{m\sigma^2/k_BT}$ respectively; how these map to physical units will depend on the system of interest, as detailed in Appendix \ref{model_details}.}

\section{Results}

{As discussed above, although a motor and the polymer bead it is bound to experience forces of equal magnitude (but acting in opposite directions), the chain connectivity means that the motors can move more freely, and tend to track along the polymer substrate.}
{The three control parameters which we vary in our simulations are the number of motors $N$, the interaction strength $K_{\rm MOT}$, and the force acting on the motors $f$. However, in what follows we also consider $\langle n \rangle$, corresponding to the mean number of motors bound to the polymer at any one time. We note that these parameters are not independent: for large forces, motors have a shorter residency time due to the finite size of the polymer substrate and if both $N$ and $f$ are large, collisions between motors can also lead to a decrease in residency time. That is to say, changing $f$ at fixed $N$ leads to a change in $\langle n \rangle$, and this relationship depends on $N$ itself -- this is examined in more detail in Appendix \ref{ssec:motors}.}

{As expected, the tracking speed of the motors in general increases with $f$, though the relationship is complex (as noted, changing $f$ leads to a change in $\langle n \rangle$ which can affect motor speed due to collisions). Within the range of parameters investigated, typically a motor will move from one polymer bead to the next in a few $\tau_{\rm LJ}$ (e.g. for $f=20$ and $N=200$ we measure a speed $\sim0.4$~beads/$\tau_{\rm LJ}$). For an $L=500$ bead polymer, motors typically spend on the order 50--100~$\tau_{\rm LJ}$ bound to it (increasing with $K_{\rm MOT}$ and decreasing with increasing $f$ and $\langle n \rangle$).}

{The motor action also changes the dynamics of the polymer. At a global scale this can be quantified by measuring the decorrelation time for the polymer end-to-end vector. For an $L=500$ bead polymer in the absence of motor action ($f=0$) a typical decorrelation time is of the order $10^5~\tau_{\rm LJ}$. Again, there is a complicated dependence on $f$ and $N$, but the motor effect leads to a \textit{decrease} in this decorrelation time by around an order of magnitude. Nevertheless, the polymer relaxation time remains much longer than the motor time-scales detailed above, i.e., the positions of motors on the polymer change quickly with respect to the polymer dynamics. This implies that, after some transient time, the properties of the polymer will settle into a non-equilibrium steady-state (NESS) which is independent of initial conditions (and the system is ergodic). In the following sections we examine these steady state polymer properties by averaging measurements over time and an ensemble of repeat simulations (50 independent simulations unless otherwise stated).}

\subsection*{Motor activity leads to polymer softening}

{In this section we study} the role of motors on the steady state conformations of the polymer substrate.  We examine the effect of different values of $\langle n \rangle$ and $f$, evolving the system until a NESS is reached. By visual inspection of the simulated trajectories we notice that, for large $\langle n \rangle$, the typical polymer conformations become more crumpled  and display kinks which do not appear in passive ($f=0$) systems (Fig.~\ref{fig:snaps}(a-b)).

\begin{figure*}[h!]
	\centering
	\includegraphics{./Figures/fig2.pdf}
	\caption{\small \textbf{(a-b)}  Comparison of conformations for a polymer of length $L =$ 500 with $\langle n \rangle=0$ and $\langle n \rangle=17$ motors attached at any time. In (b) hairpins are indicated by arrows and the purple beads mark the presence of an active motor bound to the polymer backbone. For clarity of visualisation, motors are depicted as having a bigger size with respect to the polymer beads. \textbf{(c)} Persistence length $L_p$ as a function of the average number of motors attached $\langle n \rangle$ for different values of $f$. Inset: the persistence length for different motor forces can be scaled on a master curve $L_p(x=\langle n \rangle f^2)\sim x^\alpha$ with $\alpha = -0.108 \pm  0.007$. \textbf{(d)} Similar plot showing $R_g$ as a function of $\langle n \rangle$. Again the inset shows the collapse of the curves when plotted as a function of $\langle n\rangle f^2$ with coefficient $\alpha = -0.100 \pm 0.009$. Error bars show the standard error in the mean (though these are often smaller than the points), and the lines connecting symbols are always a guide for the eye.}
	\label{fig:snaps}
\end{figure*}

To quantify the extent of polymer crumpling we measure the radius of gyration $R_g$ defined by
\begin{equation}
\label{Rg}
R_g^2 = \left\langle \frac{1}{L}\sum_{i=1}^{L} \left[ \mathbf{r}_i - \mathbf{R}_{\rm cm} \right]^2 \right\rangle ,
\end{equation}
where $\mathbf{r}_i$ is the position of the $i$-th polymer bead, $\mathbf{R}_{\rm cm}$ is the centre of mass of all polymer beads, and $\langle \, \dots \rangle$ denotes ensemble and time average (after the system reaches a steady state). We also measure the persistence length $L_p$ as defined from the bond-bond correlation function
\begin{equation}
\label{Lp}
\langle \mathbf{b}_i\cdot \mathbf{b}_{i+s} \rangle = e^{-s/L_p},
\end{equation}
where $\mathbf{b}_{i}=\left( \mathbf{r}_{i+1}-\mathbf{r}_i \right) / |\mathbf{r}_{i+1}-\mathbf{r}_i|$ represents the normalised tangent to the polymer at bead $i$. Here the average is additionally evaluated over bead index $i$. To calculate $L_p$ we perform an exponential fit of the bond-bond correlation function
(even for cases in which this displays negative values -- see below for more on this point).

As shown in Figs.~\ref{fig:snaps}\textbf{(c)} and \textbf{(d)}, both $R_g$ and $L_p$ exhibit a marked decrease with increasing $\langle n \rangle$, in agreement with our initial qualitative observation. Furthermore, the decrease in effective persistence length $L_p$ is found to have a power-law scaling
\[
L_p \sim \left( \langle n \rangle f^2 \right)^{\alpha},
\]
with $\alpha\approx -0.1$, which can also be used to collapse the different curves for $R_g$ onto a universal master curve (Fig.~\textbf{(d)} inset). In other words we find that the decrease in average size of the polymer can be well described by a scaling function of the parameter $\langle n \rangle f^2$.

\begin{figure}[t]
	\centering

	\hspace{-10pt}\includegraphics{./Figures/fig3.pdf}
	\caption{\small Radius of gyration of the polymer substrate as a function of substrate lengths $L$, for \textbf{(a)} a fixed motor force ($f\sigma/k_BT=$20) and different values for the total number of motors $N$; and \textbf{(b)} a fixed number of motors ($N=$200) and different values of the motor force $f$.
	}
	\label{fig:scaling}
\end{figure}

In order to shed more light onto the effect of the motors on the metric exponent $\nu$ relating the polymer size to its length ($R_g \sim L^{\nu}$), in Fig.~\ref{fig:scaling} we plot $R_g$ a function of $L$ for different combinations of the number of motors $N$ and force $f$. We observe that passive systems ($N=0$ or $f=0$) closely follow the expected self-avoiding statistics ($\nu\approx0.588$). More interestingly, larger values of $N$ or $f$ lead to polymer behaviours that are better fitted by $\nu \simeq$ 0.43. This value lies between that of an ideal chain $\nu = 1/2$ and that of a fully collapsed conformation $\nu=1/3$.  It is intriguing to note that a simple scaling argument can partially explain this finding using the value of $\alpha$ measured in Fig.~\ref{fig:bond-bond-corr} (see Appendix \ref{sec:scaling}).

\begin{figure}[h!]
	\centering
	\includegraphics{./Figures/fig9.pdf}
	\caption{\small \textbf{(a)} Polymer prolateness $S$ as a function of $\langle n \rangle$, for three different values of the motor force $f$. \textbf{(b)} Polymer asphericity $\Delta$ shows a similar trend. Insets show curve collapse when plotted as a function of $\langle n\rangle f^2$ in a double logarithmic scale.}
	\label{fig:aspher}
\end{figure}

To conclude this section we examine the effect of the motors on the global ``shape'' of the polymer by measuring two quantities: the asphericity $\Delta$ and prolateness $S$, which quantify the degree of spherical symmetry and the oblateness of the polymer respectively (see Appendix \ref{sec:shape} for formal definitions). In Fig.~\ref{fig:aspher}, one can readily observe that for increasing $\langle n\rangle$ and $f$ both shape descriptors decrease. This indicates that the coil is becoming more spherical, in agreement with our previous observations that the polymer effectively crumples under the action of the motors.

\subsection*{Effective softening is explained by hairpin formation}

\begin{figure*}[h!]
	\centering
	\includegraphics{./Figures/fig4.pdf}
	\caption{\small\textbf{{(a)}} Polymer bond correlation as a function of the distance along the polymer backbone $s$, at fixed motor force $f = 20$. The arrow indicates the direction of increasing $\langle n \rangle$.  \textbf{(b)}  Plot showing the minimum value of $\beta(s)$ as a function of $\langle n \rangle$. \textbf{Inset} Scaling of the minimum value of this functions, $\beta_{\rm min}$, with $\langle n\rangle f^2$. \textbf{(c)} Polymer bond correlation as a function of the distance along the polymer backbone $s$, at fixed motor force $f = 20$, $\langle n \rangle = 17 $ and different (initial) persistence lengths $L_p$. \textbf{Inset} Position of the bond correlation minimum, $L_{h}$ (defined as the value which satisfies $\beta(s=L_h)=\beta_{\rm min}$), as a function of the polymer persistence length $L_p$.
	\textbf{(d)} Plot of the polymer bond correlation for the case of motor force $f = 20$ and different values of the friction $\xi$ experienced by the motors. \textbf{Inset} Position of the bond correlation minimum $L_{h}$ as a function of the friction $\xi$.}
	\label{fig:bond-bond-corr}
\end{figure*}

Having established a crumpling effect on the global polymer conformation, we now examine the microscopic mechanism driving this process. We measure the polymer bond auto-correlation function
\[
\beta(s) \equiv  \left\langle \textbf{b}_i\cdot \textbf{b}_{i+s} \right\rangle,
\]
where as before angle brackets denote the time and ensemble average, and average over $i$.  For an equilibrium polymer $\beta(s)$ decays exponentially to zero (the bond orientation becomes uncorrelated for contour separations much larger than the persistence length). If $\beta(s)$ is close to one for a given $s$ this indicates that bonds with that separation are on average oriented parallel to each other; whereas if $\beta(s)\rightarrow-1$, then bonds with that separation are on average anti-parallel. As shown in Fig.~\ref{fig:bond-bond-corr}(a), $\beta(s)$ drops below zero and displays a minimum -- this indicates a non-negligible probability that bonds separated by about $75$ beads are anti-parallel. In other words, the polymer conformations display hairpins, or ``branches''~\cite{Michieletto2016softmatter} and the distance $s$ for which $\beta(s)$ has a minimum approximately corresponds to the length of the hairpins that are formed on the polymer. The minimum $\beta_{\rm min}$ becomes deeper as both $\langle n \rangle$ and $f$ increase (Fig.~\ref{fig:bond-bond-corr}(b)). In other words, higher values of the motor force cause a more significant anti-correlation of the bonds, i.e. more pronounced hairpins.
A fit of $\beta_{min}$ to the function $f(x) = a x^b$ gives us its scaling as a function of $\langle n \rangle$, with coefficients $a = -0.047 \pm 0.003$ and $b = 0.27 \pm 0.03$. An analogous scaling law can be found for $\beta_{min}$ as a function of the force $f$, this time with coefficients $a = -0.024 \pm  0.006 $ and $b = 0.39 \pm 0.08$.

{A possible physical mechanism underlying hairpin formation is that the motors exert localised tangential forces on the semiflexible polymer they move on. As a result the portion of the polymer ``behind'' the motor (i.e., the polymer region through which the motor has just moved) is under compression, while the portion ahead is under tension -- the situation is reminiscent of forced translocation through a pore~\cite{Saito2013}. The conformation of the polymer could then be computed by minimising its bending energy in the presence of such a localised force. This problem may be viewed as a generalisation of the Euler buckling problem~\cite{LandauLifschitzElasticity,likos_hairpin}. Although a detailed calculation is outside the scope of this work, it is reasonable to expect the combination of tensile and compressive forces to result in local hairpin formation.}

In Fig.~\ref{fig:bond-bond-corr}(c) we analyse the dependence of the bond correlation function on the intrinsic stiffness of the polymer set by $K_{\rm BEND}$ (recall that in the absence of motor activity, the persistence length $L_p \simeq K_{\rm BEND} \sigma/k_BT$). To this end, we perform three sets of simulations with $K_{\rm BEND} = 10,20$ {and} $30~{k_BT}$, for the case of $\langle n \rangle = 17$ and $f = 20$. Plotting $\beta(s)$ we observe that the \textit{position} of the minimum becomes progressively shifted towards higher values of $s$ as $K_{\rm BEND}$ increases, i.e. the hairpins  become larger (Fig.~\ref{fig:bond-bond-corr}(c) inset). Additionally, we find that the {\textit value} of $\beta(s)$ at the minimum, $\beta_{min}$, gets closer to zero for smaller values of $K_{\rm BEND}$.

In Fig.\ref{fig:bond-bond-corr}(d) we study how changing the relative impact of the active force on the motors and polymer beads affects the hairpin formation. This is achieved by changing the friction experienced by the motors (i.e. we change the viscous drag only for the motors -- this changes their effective hydrodynamic size, or could reflect a difference in interactions with macromolecular crowders compared to the polymer). Increasing the friction for the motors means that they will move less
as a result of the active force {(which is kept at a constant magnitude)}. Figure~\ref{fig:bond-bond-corr}(d) shows that as the friction for the motors is increased relative to that of the polymer the minimum of $\beta(s)$ shifts to shorter $s$, and gets deeper and narrower. This implies that greater motor friction leads to shorter and tighter hairpins. {We can understand this by considering the opposite limit of very low motor friction; in that case the motors will move more quickly relative to the polymer. Since a motor spends less time at a given polymer bead, that bead will experience the active force for less time and we can expect the effect of the motors to be diminished. The polymer will, however, still experience the active force, so we still expect a ballistic component to its motion.}

\subsection*{Motor activity gives rise to anomalous or enhanced diffusion}

We now turn our attention to the role of molecular motors on the overall polymer dynamics. To this end, we measure the mean square displacement of the polymer centre of mass
\[
{\rm MSD}_{\rm CM}(t) = \left\langle \left[ \mathbf{R}_{\rm cm}(t_0+t) - \mathbf{R}_{\rm cm}(t_0) \right]^2 \right\rangle,
\]
where the average is over start times $t_0$, and $\mathbf{R}_{\rm cm}(t)$ is the position of the centre of mass at time $t$.

Figure~\ref{fig:Fig_diff} shows that at short times the system is diffusive as the propulsion exerted by the motors needs a characteristic time in order to provide a displacement comparable to a Brownian displacement. Following the Rouse model~\cite{rouse} ${\rm MSD}_{\rm CM}(t) = (6D_0/L) t$, where $D_0 = k_B T/\xi$ is the diffusion coefficient of a single monomer. This behaviour is recovered at all times for $\langle n \rangle \rightarrow 0$.

\begin{figure}[h!]
	\centering
	\includegraphics[width=.5\textwidth]{./Figures/fig5.pdf}
	\caption{\small \textbf{(a)} Mean squared displacement of the polymer centre of mass as a function of time. Each set of points shows results for a different value of $\langle n \rangle$ at fixed $f=20~k_BT \sigma^{-1}$. Higher values of $\langle n \rangle$ induce a faster diffusivity of the polymer.  \textbf{(b)} Similar plot but now each set of points shows results for a different value of $f$, with the total number of motors fixed at $N=400$.
}
	\label{fig:Fig_diff}
\end{figure}

\begin{table}[h!]
	\begin{center}
		\begin{tabular}{ |c|c|c| }
			\hline
			\multicolumn{3}{|c|}{\textbf{Intermediate times: super-diffusive}}\\
			\hline
			$\langle n \rangle$ & $\alpha$ & $\mu$ \\
			\hline
			17 & $1.720\pm 0.002$ & $3.1\cdot10^{-5}\pm0.5\cdot10^{-6}$ \\
			9 & $1.690 \pm 0.001$ & $1.7\cdot10^{-5} \pm 0.2\cdot10^{-6}$ \\
			5 & $1.542 \pm 0.006 $ & $2.6\cdot10^{-5} \pm 1.3\cdot10^{-6}$\\
			\hline
			\multicolumn{3}{|c|}{\textbf{Long times: diffusive}}\\
			\hline
			$\langle n \rangle$ & $\alpha$ & $D$ \\
			\hline
			17 & $1.025 \pm0.001 $& $2.7\cdot10^{-2} \pm 2.6\cdot 10^{-4}$ \\
			9 & $1.020 \pm 0.002 $ & $2.0\cdot10^{-2} \pm 0.0 6\cdot10^{-4}$\\
			5 & $1.090 \pm 0.0001$& $0.4\cdot10^{-2} \pm 0.08\cdot10^{-4}$\\
			\hline
		\end{tabular}
	\end{center}
	\caption{Diffusion properties of the centre of mass of the polymer are measured by fitting to the MSD for different numbers of motors. Fits are evaluated using the function $f(t) = a\cdot x^{\alpha}$, and the diffusion coefficient (or mobility $\mu$) is then $D = a/6$ following the Rouse model.}
	\label{table:msd}
\end{table}

At intermediate times, the polymer exhibits super-diffusive motion, i.e. ${\rm MSD}_{\rm CM}(t) \sim t^{\alpha}$ with $\alpha \sim 1.6 $.  Finally, at large times we observed enhanced diffusion, i.e. $\alpha\approx1$ but the diffusion constant increases with motor activity ($f$ or $\langle n \rangle$).  In Table~\ref{table:msd} we report the measured diffusion constant $D$ for the different cases considered.
One can readily observe that the exponent $\alpha$ is roughly independent of $\langle n \rangle$ and $f$, whereas the diffusion coefficient increases with $\langle n \rangle$.  Moreover, by comparing fits at short and long times for each value of $\langle n \rangle$, we find that the crossover between the super-diffusive and diffusive behaviour occurs at progressively larger times as $\langle n \rangle$ decreases.
Together, this implies that the motor action leads to a persistent motion of the centre of mass of the polymer; at longer times the velocity becomes uncorrelated and the motion is again diffusive.

{These three regimes are typical of other active systems, e.g. in polymers where an active force directed along the backbone is applied to each polymer bead~\cite{bianco_loca,Anand2018} and in  Active Brownian Particles (ABPs)~\cite{Cates2012}. At intermediate times, the motion of an ABP is dictated by a force which typically maintains its direction for a characteristic time, leading to super-diffusion. In our model, a correlation on the scale of the whole polymer is set by the end-to-end vector.}

\subsection*{Untying Knots with Molecular Motors}

Finally, we discuss a possible application of our model by studying the knotting probability of a polymer in non-equilibrium conditions.  We start from the observation that a freely diffusing linear polymer in a good solvent may display a physical knot -- i.e. a knot within an open segment of the polymer contour -- with probability $P_{\rm knot}$. It has recently been shown~\cite{matthews2012, poier2014, Coronel2017} that $P_{\rm knot}$ displays a non-monotonic dependence on the stiffness of the substrate with a maximum at around $L_p \simeq 5-10 \sigma$. In light of the results presented above -- i.e., motors drive a decrease in the apparent persistence length of the substrate -- we reason that the presence of the motors should lead to a shift in the knotting probability as a function of substrate persistence length. Specifically, we expect that, due to the action of the motors, the maximum in $P_{\rm knot}$ will be located at larger values of the intrinsic persistence length of the substrate, $l_p$ (to be clear, the lower case $l_p$ refers to the persistence length the polymer would have in the absence of motors).

In order to test this hypothesis we perform simulations of a freely diffusing linear polymer 500 beads long, and measure the probability of forming physical knots as a function of $l_p$ by adopting the algorithm described in Refs.~\cite{Tubiana2011,Tubiana2011prl} and publicly available at \url{http://kymoknot.sissa.it/}~\cite{Tubiana2018}. By comparing $P_{\rm knot}(l_p)$ for cases with and without motors we find that, surprisingly, the location of the maximum of $P_{\rm knot}$ is unchanged by the action of the motors, i.e. $l_p \simeq 5-10 \sigma$. The most apparent effect of the motor activity is a global reduction in the knotting probability across the full spectrum of bending rigidities studied (Fig.~\ref{fig:knots}).
In light of this finding we argue that the effect of the motors goes beyond a simple softening of the polymer backbone, and that subtler non-equilibrium effects must be at play.
We also speculate that the reduced knotting probability due to motor activity may be relevant for the case of DNA under the action of motors such as polymerase. According to our model, the action of these motors would naturally have the beneficial consequence of reducing the occurrence of knots and self-entanglements {\textit in vivo}~\cite{Bates2005}.

\begin{figure}[t!]
	\centering
	\includegraphics{./Figures/fig6.pdf}
	\vspace{-0,35 cm}
	\caption{Molecular motors acting on a polymeric substrate reduce the knotting probability $P_{\rm knot}(l_p)$ of a free chain in a good solvent. Here $l_p$ is the intrinsic persistence length which would be measured in the absence of motors. The simulations are performed on a 500 bead long chain and the points are averages over $10^4$ independent configurations (error bars show the standard error in the mean).}
	\label{fig:knots}
\end{figure}

\section{Conclusions}

In this work we performed Langevin dynamics simulations of a coarse grained model to study the effect of molecular motors on a polymer substrate. We examined both the static and dynamic properties of polymers in dilute conditions while varying the number of motors and the magnitude of the force which they generate.

Our main result was that the effect of the motors is to reduce the overall size of the polymer coil by increasing the likelihood of double-folded segments, or hairpins. This ``softening'' of the backbone manifests as a reduction of the measured persistence length of the polymer. We found that the data collapse onto a single curve when plotted as a function of $\langle n\rangle f^2$, revealing a power-law scaling $R_g \sim (\langle n\rangle f^2)^{-0.1}$ and $L_p \sim (\langle n\rangle f^2)^{-0.1}$ (insets of Fig.~\ref{fig:snaps}\textbf{(c)} and \textbf{(d)}).  We also performed simulations of polymers with different lengths and obtained the explicit values of the metric exponent $\nu$ (where $R_g \sim L^{\nu}$) as a function of $\langle n\rangle$ and $f$. These confirm that at equilibrium the polymers obey self-avoiding statistics, whereas motor activity leads to an exponent between those expected for a crumpled globule and an ideal chain (Fig.~\ref{fig:scaling}).
{A reduction in $R_g$ and $L_p$ has also been observed in models without explicit motors, where an active force directed tangentially to the polymer is added directly to each polymer bead~\cite{bianco_loca,Anand2018}. In that case the softening arises because, for an initially straight segment of the polymer (e.g. of length equal to the persistence length), small fluctuations in backbone direction will lead to components of the tangential force which are perpendicular to the end-to-end vector for the segment; this then drives the polymer into large smooth curves~\cite{bianco_loca}. Those studies, however, did not find hairpins (the decay in bond correlations remained exponential~\cite{Anand2018}), suggesting the present case involves a different mechanism. }

We {also} found that the motors affect the dynamics of the polymer. At intermediate times the centre of mass displays super-diffusive behaviour while at longer times it follows an enhanced diffusion. Both regimes are controlled by motor activity (number of bound motors and force applied).  While a super-diffusive regime was to be expected, our simulation have provided a detailed quantification of how that can be tuned through the action of the motors (Fig.~\ref{fig:Fig_diff} and Table~\ref{table:msd}).
The effect of the motors on the knotting probability of the polymer was subtler. Unexpectedly, $P_{\rm knot}(l_p)$ did not simply show a shift due to the change in effective persistence length, but rather there was a global reduction in the steady state knotting probability across a range of substrate rigidities (Fig.~\ref{fig:knots}). {We can speculate that since the motor action promotes polymer configurations with hairpin bends, but suppresses knotted ones, then these two types of configuration are incompatible. An interesting future study would be to examine knotting probabilities in other active polymer models where hairpins are not observed.}

Throughout this paper we have limited ourselves to the simple case in which all of the motors travel in the same direction along the substrate.  While this is the case for some systems, e.g. kinesin on microtubules, in others the motors can travel in either direction.  For example, RNA polymerase follows the direction of genes, depending on the specifically oriented DNA sequences in gene promoters. This could be readily incorporated into our model in the future by either considering a mixture of two species of motors which move in different directions on the substrate, or by setting some directional binding sites. {Also, we have considered the case of motors which are the same size as the polymer beads -- in reality RNA polymerase is a large molecule of around 30~nm in size. However, as our results are mostly affected by the magnitude of the force, we do not expect motor size to have a qualitative affect. Nevertheless it would be interesting to examine systems with larger motors in the future.}

Another interesting scenario would be to consider a mixture of passive and active polymers. Indeed, an activity-driven phase separation has been observed in another non-equilibrium polymer system -- where two species of polymer were held at different temperatures~\cite{Smrek2017}. It would be interesting to see whether such demixing behaviour can be recovered when the system is driven away from equilibrium by a more biologically relevant mechanism, such as that presented here.

Our model naturally lends itself to the study of the response of polymer substrates to the action of molecular motors. We therefore expect that in the future is could be employed to address more specific biological questions, such as understanding the role of motors in the buckling of actin networks~\cite{Broedersz2014,Lenz2014}, or in the bending and active collective motion of microtubules~\cite{Head2011,Pearce2018} by suitably tuning the rigidity of the substrate and processivity of the motors.

Finally we stress that the model presented here is among the first to explicitly account for the presence of physical motors; as a consequence, it faithfully captures the feedback between the action of the motors and the change in substrate-motor interaction due to, for instance, alterations in local polymer conformation. For this reason, we believe that our model will be pivotal in studying complex active collective phenomena where interactions between motors are mediated through the target substrate.

\subsection*{Acknowledgements}

This research was partly funded by the European Training Network COLLDENSE (N2020-MCSA-ITN-2014) grant number 642774.
DMa, DMi and CAB acknowledge funding from ERC (CoG 648050 THREEDCELLPHYSICS). \\
DMi, CNL and EL would also like to acknowledge the contribution and networking support by the ``European Topology Interdisciplinary Action'' (EUTOPIA) CA17139.
\footnotesize{

\providecommand*{\mcitethebibliography}{\thebibliography}
\csname @ifundefined\endcsname{endmcitethebibliography}
{\let\endmcitethebibliography\endthebibliography}{}
\begin{mcitethebibliography}{47}
\providecommand*{\natexlab}[1]{#1}
\providecommand*{\mciteSetBstSublistMode}[1]{}
\providecommand*{\mciteSetBstMaxWidthForm}[2]{}
\providecommand*{\mciteBstWouldAddEndPuncttrue}
  {\def\unskip.}{\unskip.}}
\providecommand*{\mciteBstWouldAddEndPunctfalse}
  {\let\unskip.}\relax}
\providecommand*{\mciteSetBstMidEndSepPunct}[3]{}
\providecommand*{\mciteSetBstSublistLabelBeginEnd}[3]{}
\providecommand*{\unskip.}}{}
\mciteSetBstSublistMode{f}
\mciteSetBstMaxWidthForm{subitem}
{(\emph{\alph{mcitesubitemcount}})}
\mciteSetBstSublistLabelBeginEnd{\mcitemaxwidthsubitemform\space}
{\relax}{\relax}

\bibitem[Sanchez \emph{et~al.}(2012)Sanchez, Chen, Decamp, Heymann, and
  Dogic]{Sanchez2012nature}
T.~Sanchez, D.~T.~N. Chen, S.~J. Decamp, M.~Heymann and Z.~Dogic,
  \emph{Nature}, 2012, \textbf{491}, 1--5\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Alberts \emph{et~al.}(2014)Alberts, Johnson, Lewis, Morgan, and
  Raff]{Alberts2014}
B.~Alberts, A.~Johnson, J.~Lewis, D.~Morgan and M.~Raff, \emph{{Molecular
  Biology of the Cell}}, Taylor {\&} Francis, 2014, p. 1464\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Calladine \emph{et~al.}(1997)Calladine, Drew, Luisi, Travers, and
  Bash]{Calladine1997}
C.~R. Calladine, H.~Drew, F.~B. Luisi, A.~A. Travers and E.~Bash,
  \emph{{Understanding DNA: the molecule and how it works}}, Elsevier Academic
  Press, 1997, vol.~1\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Le \emph{et~al.}(2013)Le, Imakaev, Mirny, and Laub]{Le2013}
T.~B. Le, M.~V. Imakaev, L.~A. Mirny and M.~T. Laub, \emph{Science (New York,
  N.Y.)}, 2013, \textbf{342 VN -}, 731--734\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Terakawa \emph{et~al.}(2017)Terakawa, Bisht, Eeftens, Dekker, Haering,
  and Greene]{Terakawa2017}
T.~Terakawa, S.~Bisht, J.~M. Eeftens, C.~Dekker, C.~H. Haering and E.~C.
  Greene, \emph{Science}, 2017, \textbf{676}, eaan6516\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Brackley \emph{et~al.}(2016)Brackley, Johnson, Bentivoglio, Corless,
  Gilbert, Gonnella, and Marenduzzo]{Brackley2016supercoil}
C.~A. Brackley, J.~Johnson, A.~Bentivoglio, S.~Corless, N.~Gilbert, G.~Gonnella
  and D.~Marenduzzo, \emph{Phys. Rev. Lett.}, 2016, \textbf{117}, 018101\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Naughton \emph{et~al.}(2013)Naughton, Avlonitis, Corless, Prendergast,
  Mati, Eijk, Cockroft, Bradley, Ylstra, and Gilbert]{Naughton2013}
C.~Naughton, N.~Avlonitis, S.~Corless, J.~G. Prendergast, I.~K. Mati, P.~P.
  Eijk, S.~L. Cockroft, M.~Bradley, B.~Ylstra and N.~Gilbert, \emph{Nat.
  Struct. Mol. Biol.}, 2013, \textbf{20}, 387--395\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Lenz(2014)]{Lenz2014}
M.~Lenz, \emph{Phys. Rev. X}, 2014, \textbf{4}, 1--9\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Broedersz and MacKintosh(2014)]{Broedersz2014}
C.~P. Broedersz and F.~C. MacKintosh, \emph{Rev. Mod. Phys.}, 2014,
  \textbf{86}, 995--1036\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Tjhung \emph{et~al.}(2012)Tjhung, Marenduzzo, and Cates]{Tjhung2012}
E.~Tjhung, D.~Marenduzzo and M.~E. Cates, \emph{Proc. Natl. Acad. Sci. USA},
  2012, \textbf{109}, 12381--12386\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Herbert \emph{et~al.}(2008)Herbert, Greenleaf, and Block]{Herbert2008}
K.~M. Herbert, W.~J. Greenleaf and S.~M. Block, \emph{Annu Rev Biochem}, 2008,
  \textbf{77}, 149--176\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Michaelis and Treutlein(2013)]{michaelis2013single}
J.~Michaelis and B.~Treutlein, \emph{Chemical reviews}, 2013, \textbf{113},
  8377--8399\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Ganji \emph{et~al.}(2018)Ganji, Shaltiel, Bisht, Kim, Kalichava,
  Haering, and Dekker]{Ganji2018}
A.~M. Ganji, I.~A. Shaltiel, S.~Bisht, E.~Kim, A.~Kalichava, C.~H. Haering and
  C.~Dekker, \emph{Science}, 2018, \textbf{7831}, 1--9\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Bianco \emph{et~al.}(2018)Bianco, Locatelli, and
  Malgaretti]{bianco_loca}
V.~Bianco, E.~Locatelli and P.~Malgaretti, \emph{Phys. Rev. Lett.}, 2018,
  \textbf{121}, 217802\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Saintillan \emph{et~al.}(2018)Saintillan, Shelley, and
  Zidovska]{SaintillanMotor}
D.~Saintillan, M.~J. Shelley and A.~Zidovska, \emph{Proceedings of the National
  Academy of Sciences}, 2018, \textbf{115}, 11442--11447\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Harder \emph{et~al.}(2014)Harder, Valeriani, and Cacciuto]{Harder2014}
J.~Harder, C.~Valeriani and A.~Cacciuto, \emph{Phys. Rev. E}, 2014,
  \textbf{90}, 062312\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Jiang and Hou(2014)]{Jiang2014}
H.~Jiang and Z.~Hou, \emph{Soft Matter}, 2014, \textbf{10}, 1012--1017\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Isele-Holder \emph{et~al.}(2015)Isele-Holder, Elgeti, and
  Gompper]{Isele-Holder2015}
R.~E. Isele-Holder, J.~Elgeti and G.~Gompper, \emph{Soft Matter}, 2015,
  \textbf{11}, 7181--7190\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Eisenstecken \emph{et~al.}(2016)Eisenstecken, Gompper, and
  Winkler]{Eisenstecken2016}
T.~Eisenstecken, G.~Gompper and R.~G. Winkler, \emph{Polymers}, 2016,
  \textbf{8}, year\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Anand and Singh(2018)]{Anand2018}
S.~K. Anand and S.~P. Singh, \emph{Phys. Rev. E}, 2018, \textbf{98},
  042501\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Winkler \emph{et~al.}(2017)Winkler, Elgeti, and Gompper]{Winkler2017}
R.~G. Winkler, J.~Elgeti and G.~Gompper, \emph{Journal of the Physical Society
  of Japan}, 2017, \textbf{86}, 101014\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Schaller \emph{et~al.}(2010)Schaller, Weber, Semmrich, Frey, and
  Bausch]{Schaller2010}
V.~Schaller, C.~Weber, C.~Semmrich, E.~Frey and A.~R. Bausch, \emph{Nature},
  2010, \textbf{467}, 73\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Sanchez \emph{et~al.}(2012)Sanchez, Chen, Decamp, Heymann, and
  Dogic]{Sanchez2012}
T.~Sanchez, D.~T.~N. Chen, S.~J. Decamp, M.~Heymann and Z.~Dogic,
  \emph{Nature}, 2012, \textbf{491}, 1--5\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Nedelec and Foethke(2007)]{nedelec2007}
F.~Nedelec and D.~Foethke, \emph{New Journal of Physics}, 2007, \textbf{9},
  427\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Freedman \emph{et~al.}(2017)Freedman, Banerjee, Hocky, and
  Dinner]{Freedman2017}
S.~L. Freedman, S.~Banerjee, G.~M. Hocky and A.~R. Dinner, \emph{Biophysical
  Journal}, 2017, \textbf{113}, 448--460\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Ravichandran \emph{et~al.}(2017)Ravichandran, Vliegenthart,
  Saggiorato, Auth, and Gompper]{Ravichandran2017}
A.~Ravichandran, G.~A. Vliegenthart, G.~Saggiorato, T.~Auth and G.~Gompper,
  \emph{Biophysical Journal}, 2017, \textbf{113}, 1121 -- 1132\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Gupta \emph{et~al.}(2019)Gupta, Chaudhuri, and Chaudhuri]{Gupta2019}
N.~Gupta, A.~Chaudhuri and D.~Chaudhuri, \emph{Phys. Rev. E}, 2019,
  \textbf{99}, 042405\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Coronel \emph{et~al.}(2017)Coronel, Orlandini, and
  Micheletti]{Coronel2017}
L.~Coronel, E.~Orlandini and C.~Micheletti, \emph{Soft Matter}, 2017,
  \textbf{13}, 4260--4267\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Plimpton(1995)]{LAMMPS}
S.~Plimpton, \emph{J Comp Phys}, 1995, \textbf{117}, 1--19\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Michieletto(2016)]{Michieletto2016softmatter}
D.~Michieletto, \emph{Soft Matter}, 2016, \textbf{12}, 9485--9500\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Saito and Sakaue(2013)]{Saito2013}
T.~Saito and T.~Sakaue, \emph{Phys. Rev. E}, 2013, \textbf{88}, 042606\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Landau \emph{et~al.}(1986)Landau, Kosevich, Pitaevskii, and
  Lifshitz]{LandauLifschitzElasticity}
L.~D. Landau, A.~Kosevich, L.~Pitaevskii and E.~Lifshitz, \emph{Theory of
  Elasticity}, Butterworth-Heinemann, Oxford, 1986\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Wynveen and Likos(2010)]{likos_hairpin}
A.~Wynveen and C.~N. Likos, \emph{Soft Matter}, 2010, \textbf{6},
  163--171\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Prince E.~Rouse()]{rouse}
J.~Prince E.~Rouse, \emph{The Journal of Chemical Physics}, \textbf{20},
  year\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Cates(2012)]{Cates2012}
M.~E. Cates, \emph{Rep. Prog. Phys.}, 2012, \textbf{75}, 042601\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Matthews \emph{et~al.}(2012)Matthews, Louis, and Likos]{matthews2012}
R.~Matthews, A.~a. Louis and C.~N. Likos, \emph{ACS Macro Letters}, 2012,
  1352--1356\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Poier \emph{et~al.}(2014)Poier, Likos, and Matthews]{poier2014}
P.~Poier, C.~N. Likos and R.~Matthews, \emph{Macromolecules}, 2014,
  \textbf{47}, 3394--3400\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Tubiana \emph{et~al.}(2011)Tubiana, Orlandini, and
  Micheletti]{Tubiana2011}
L.~Tubiana, E.~Orlandini and C.~Micheletti, \emph{Prog. Theor. Phys. Suppl.},
  2011, \textbf{191}, 192--204\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Tubiana \emph{et~al.}(2011)Tubiana, Orlandini, and
  Micheletti]{Tubiana2011prl}
L.~Tubiana, E.~Orlandini and C.~Micheletti, \emph{Phys. Rev. Lett.}, 2011,
  \textbf{107}, 188302\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Tubiana \emph{et~al.}(2018)Tubiana, Polles, Orlandini, and
  Micheletti]{Tubiana2018}
L.~Tubiana, G.~Polles, E.~Orlandini and C.~Micheletti, \emph{EPJ E}, 2018,
  \textbf{41}, 72\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Bates and Maxwell(2005)]{Bates2005}
A.~Bates and A.~Maxwell, \emph{{DNA topology}}, Oxford University Press,
  2005\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Smrek and Kremer(2017)]{Smrek2017}
J.~Smrek and K.~Kremer, \emph{Phys. Rev. Lett.}, 2017, \textbf{118}, 1--5\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Head \emph{et~al.}(2011)Head, Briels, and Gompper]{Head2011}
D.~a. Head, W.~Briels and G.~Gompper, \emph{BMC biophysics}, 2011, \textbf{4},
  18\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Pearce \emph{et~al.}(2018)Pearce, Heil, Jensen, Jones, and
  Prokop]{Pearce2018}
S.~P. Pearce, M.~Heil, O.~E. Jensen, G.~W. Jones and A.~Prokop, \emph{Bulletin
  of Mathematical Biology}, 2018, \textbf{80}, 3002--3022\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Brackley \emph{et~al.}()Brackley, Morozov, and
  Marenduzzo]{BrackleyTwistable}
C.~A. Brackley, A.~N. Morozov and D.~Marenduzzo, \emph{The Journal of Chemical
  Physics}, \textbf{140}, year\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[de~Gennes(1979)]{DeGennes}
P.~G. de~Gennes, \emph{{Scaling Concepts in Polymer Physics }}, Cornell
  University Press, 1979\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\bibitem[Narros \emph{et~al.}(2013)Narros, Moreno, and Likos]{Narros2013}
A.~Narros, A.~J. Moreno and C.~N. Likos, \emph{Macromolecules}, 2013,
  \textbf{46}, 3654--3668\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\unskip.}
\end{mcitethebibliography}

}

\newpage