# Coherent two-dimensional THz magnetic resonance spectroscopies for molecular
magnets: Analysis of Dzyaloshinskii–Moriya interaction

Jiaji Zhang  0000-0003-2978-274X<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>Yoshitaka Tanimura  0000-0002-7913-054X
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>Department of
Chemistry, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan

(Last updated: )

###### Abstract

To investigate the novel quantum dynamic behaviors of magnetic materials that
arise from complex spin–spin interactions, it is necessary to probe the
magnetic response at a speed greater than the spin-relaxation and dephasing
processes. Recently developed two-dimensional (2D) terahertz magnetic
resonance (THz-MR) spectroscopy techniques use the magnetic components of
laser pulses, and this allows investigation of the details of the ultrafast
dynamics of spin systems. For such investigations, quantum treatment—not only
of the spin system itself but also of the environment surrounding the spin
system—is important. In our method, based on the theory of multidimensional
optical spectroscopy, we formulate nonlinear THz-MR spectra using an approach
based on the numerically rigorous hierarchical equations of motion. We conduct
numerical calculations of both linear (1D) and 2D THz-MR spectra for a linear
chiral spin chain. The pitch and direction of chirality (clockwise or
anticlockwise) are determined by the strength and sign of the
Dzyaloshinskii–Moriya interaction (DMI). We show that not only the strength
but also the sign of the DMI can be evaluated through the use of 2D THz-MR
spectroscopic measurements, while 1D measurements allow us to determine only
the strength.

##  I Introduction

Electron paramagnetic resonance (EPR) and nuclear magnetic resonance (NMR)
have a long history as typical examples of two-dimensional (2D) spectroscopy
techniques that measure the varying time intervals of magnetic pulse trains
applied to electron or nuclear spin systems. [1, 2] Although these techniques
are powerful means for structural analysis of organic and inorganic materials,
it is difficult to apply them to the investigation of spin dynamics because
the time resolution of the magnetic pulses involved is limited to the order of
microseconds.

Recently, coherent terahertz (THz) magnetic resonance (MR) spectroscopy was
developed using magnetic pulses in the subpicosecond range generated by the
magnetic-field component of THz light. [3, 4, 5] Just as ultrafast laser
spectroscopy has made it possible to study the electronic excitation-state
dynamics and intra- and intermolecular vibrational motions of complex
molecular systems, [6] THz-MR spectroscopy has opened up the possibility of
investigating strongly correlated spin dynamics in molecular magnets; this
depends on the complex spin–spin interactions and the configurations of spins
in the molecular environment, which cause dephasing and relaxation of the spin
system.

The quasi-ferromagnetic (FM) and antiferromagnetic (AFM) precession modes in
$\rm{YFeO_{3}}$ (YFO) that arise from the antisymmetric spin–orbit coupling
[7] have been observed as free induction decay (FID) signals. [8, 9] Spin-wave
excitation in AFM $\rm{NiO}$ has been detected with the aid of the Faraday
effect. [10, 11] The skyrmion, which is a quasiparticle composed of vortex-
like spin orientations, has been detected in several FM and AFM materials with
thin-film structures. [12, 13, 14] THz electron spin resonance (ESR) and Hall
conductivity spectroscopy have been employed to investigate the topological
Hall effect and phase transitions. [15, 16, 3, 17]

Although THz-MR and THz-ESR spectroscopic approaches based on linear response
theory have been successful for the classification of complex spin states in
condensed phases, [18, 19, 20] it is unclear whether the spin states that are
investigated are quantum-mechanically entangled, and whether the width of the
peaks that are measured arise from the relaxation or dephasing processes,
because these peaks are usually broad and overlap. Thus, an extension to
coherent 2D THz-MR spectroscopy was developed. The observed 2D THz-MR spectrum
for a YFO crystal made it possible to illustrate characteristic nonlinear spin
responses, such as double-quantum (2Q) coherence and second-harmonic
generation (SHG). [21]

As already illustrated in 2D optical spectroscopy, [22, 23] 2D THz-MR
spectroscopy is not only useful for identifying complex spin interactions but
also for monitoring the complex quantum spin dynamics under the influence of
relaxation and dephasing arising from the environment at the femtosecond
scale. Theoretical input regarding the complex profiles of spin–spin
interactions is important for analyzing these 2D spectra under ultrafast
nonlinear processes. In particular, this includes those in molecular magnets,
which play a central role in spintronics and next-generation information
technologies. [24, 25, 26]

In this paper, we provide a comprehensive theoretical framework for both 1D
and 2D THz-MR measurements based on the response-function theory developed for
nonlinear optical spectroscopy techniques. To illustrate our approach more
closely, we employ a chiral spin chain described as a Heisenberg model with
exchange coupling and Dzyaloshinskii–Moriya interaction (DMI) arising from the
antisymmetric spin–orbit coupling. [27, 28, 29] Such anisotropy of the spin
system leads to a series of unique magnetic and optical properties. [30]
Unlike with conventional FM materials, in addition to the magnetic anisotropy,
a non-centrosymmetric structure arises as the result of the chiral ordering of
spins. As a result, a series of new phenomena arising from the nonlinear
magnetic response, such as magneto-chiral dichroism and magnetization-induced
SHG, were observed. [31, 32, 33, 34]

To investigate the properties of such materials as devices, it is necessary to
investigate ultrafast spin dynamics under time-dependent external fields, in
which quantum coherence and entanglement—not only among spins but also between
spins and the environment—play a significant role. [2, 35] Thus, we include a
harmonic heat bath of finite temperature in the spin system. The number of
degrees of freedom of the heat bath is then reduced to obtain a time-
irreversible equation of motion describing the effects of thermal fluctuations
and dissipation. [36, 37] Because the motion of the spins is much faster than
the thermal noise arising from the heat bath, the heat bath must be treated in
a non-Markovian manner. Thus, the reduced equations of motion derived using
Markovian approximations and other assumptions—such as the Bloch, Lindblad,
and Redfield equations—are not suitable for the description of such ultrafast
spin dynamics. We then employ the numerically “exact” hierarchical equations
of motion (HEOM) approach, which can be used to treat non-Markovian and non-
perturbative system–bath interactions at finite temperatures. [38, 39, 40, 41]

To demonstrate the applicability of the present theory, 1D and 2D THz-MR
spectra were calculated for the chiral spin model with different DMI strengths
describing the pitch and direction of chirality (clockwise or
counterclockwise). We show that in 1D spectra, only the absolute strength of
the DMI can be evaluated, whereas in 2D spectra, the sign of the DMI, which
determines the direction of the chirality, can also be determined. While
neutron-scatting techniques have been used to determine the structures of
chiral materials, the present results indicate the possibility of determining
the DMI through spin-dynamic processes. This finding should be valuable for
the design of spintronic materials.

The rest of this paper is organized as follows. In Sec. II, we formulate
linear (1D) and 2D THz-MR spectra based on the response-function theory. In
Sec. III, we introduce the Heisenberg model with the DMI coupled to the
harmonic heat bath. The HEOM are then presented. Numerical results and some
discussion are presented in Sec. IV. Finally, Sec. V is devoted to our
conclusions.

##  II Magnetic susceptibilities in 2D THz-MR spectroscopy

We consider a spin system coupled to a bath system driven by an external
magnetic field $B(s)$. The total Hamiltonian is expressed as

$\hat{H}^{\prime}(s)=\hat{H}_{tot}+B(s)\,\hat{M},$ (1)

where $\hat{H}_{tot}$ is the Hamiltonian of a composite system and $\hat{M}$
is the polarization operator for magnetic fields.

The observable in a magnetic measurement at time $t$ is expressed as
$M(t)\equiv\langle\hat{M}(t)\rangle-\langle\hat{M}\rangle$, where $\hat{M}(t)$
is the Heisenberg representation of $\hat{M}$ for $\hat{H}^{\prime}(s)$. The
thermal average for any operator $\hat{A}$ is defined as
$\langle\hat{A}\rangle\equiv{\rm{tr}}\\{\hat{A}\hat{\rho}_{tot}^{eq}\\}$, with
the equilibrium density operator expressed as $\hat{\rho}_{tot}^{eq}$. We can
also express this as
${M}(t)={\rm{tr}}\\{\hat{M}\mathcal{G}^{\prime}(t)\hat{\rho}_{tot}^{eq}\\}-{\rm{tr}}\\{\hat{M}\hat{\rho}_{tot}^{eq}\\}$,
where the Liouville operator is defined as

$\displaystyle\mathcal{G}^{\prime}(t)\hat{\rho}_{tot}(0)\equiv{\underleftarrow{\exp}}\left[-\frac{i}{\hbar}\int_{0}^{t}{{\rm{d}}}s^{\prime}\hat{H}^{\prime}(s^{\prime})\right]\hat{\rho}_{tot}(0)$
$\displaystyle\times{\underrightarrow{\exp}}\left[\frac{i}{\hbar}\int_{0}^{t}{{\rm{d}}}s^{\prime}\hat{H}^{\prime}(s^{\prime})\right].$
(2)

Here, the arrows indicate time-ordered exponentials.

Experiments such as 2D NMR and 2D EPS measurements [1, 2] are conducted using
multiple magnetic pulses with finite time widths. In such experiments, the
excitation by the external field is non-perturbative, and the desired spin
dynamics are investigated by designing the profiles of pulse trains, as in the
cases of spin-echo and correlation spectroscopy measurements. Theoretically,
these signals are obtained by integrating the equations of motion for a spin
system, such as the Bloch, [42] Redfield, [43] or stochastic Liouville
equations, [44] or the HEOM, [40, 45, 46] under a sequence of magnetic pulses.

As in the case of coherent optical laser spectroscopies, the excitations of
coherent THz-MR spectroscopy are assumed to be impulsive, which allows us to
measure the signals in different orders of the field-system interactions
separately. Thus, we can employ a response-function theory developed for
ultrafast nonlinear laser spectroscopy. [6] Up to the third order, the signal
is then expressed as

$\begin{split}M(t)&=\int_{0}^{t}{{\rm{d}}}s\,\chi_{1}(s)\,B(t-s)\\\
&+\int_{0}^{t}{{\rm{d}}}s_{1}\int_{0}^{s_{1}}{{\rm{d}}}s_{2}\,\chi_{2}(s_{1},s_{2})\,B(t-s_{1})\,B(t-s_{2})\\\
&+\int_{0}^{t}{{\rm{d}}}s_{1}\int_{0}^{s_{1}}{{\rm{d}}}s_{2}\,\int_{0}^{s_{2}}{{\rm{d}}}s_{3}\,\chi_{3}(s_{1},s_{2},s_{3})\,\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times
B(t-S_{1})\,B(t-s_{2})\,B(t-s_{3}),\end{split}$ (3)

where the linear, second-order, and third-order response functions are defined
as [40]

$\begin{split}\chi_{1}(s)&\equiv{-\frac{i}{\hbar}}\langle[\hat{M}(s),\hat{M}]\rangle\\\
&={-\frac{i}{\hbar}}{\rm{tr}}\left\\{\hat{M}\mathcal{G}(s)\,\hat{M}^{\times}\hat{\rho}_{tot}^{eq}\right\\},\end{split}$
(4)
$\begin{split}\chi_{2}(s_{1},s_{2})&\equiv{-\frac{1}{\hbar^{2}}}\langle[[\hat{M}(s_{1}+s_{2}),\hat{M}(s_{1})],\,\hat{M}]\rangle\\\
&={-\frac{1}{\hbar^{2}}}{\rm{tr}}\left\\{\hat{M}\,\mathcal{G}(s_{2})\,\hat{M}^{\times}\mathcal{G}(s_{1})\hat{M}^{\times}\hat{\rho}_{tot}^{eq}\right\\},\end{split}$
(5)

and

$\begin{split}\chi_{3}&(s_{1},s_{2},s_{3})\\\
&\equiv{\frac{i}{\hbar^{3}}}\langle[[[\hat{M}(s_{1}+s_{2}+s_{3}),\hat{M}(s_{1}+s_{2})],\,\hat{M}(s_{1})],\,\hat{M}]\rangle\\\
&={\frac{i}{\hbar^{3}}}{\rm{tr}}\left\\{\hat{M}\,\mathcal{G}(s_{3})\,\hat{M}^{\times}\,\mathcal{G}(s_{2})\,\hat{M}^{\times}\mathcal{G}(s_{1})\hat{M}^{\times}\hat{\rho}_{tot}^{eq}\right\\}.\end{split}$
(6)

Here, $\mathcal{G}(s)$ is the Liouvillian without magnetic interaction, which
is defined from $\mathcal{G}^{\prime}(s)$ with $B(s)=0$, and we have
introduced the superoperator
$\hat{M}^{\times}{{\hat{\rho}}}\equiv[\hat{M},{\hat{\rho}}]$. The right-hand
side (rhs) of the second line in each of the above equations allows us to
employ the equations of motion to calculate the response functions and give us
an intuitive picture of the higher-order optical processes. [40, 41] For
example, the rhs of the second line of Eq. (4) can be read from right to left
as follows. The total system is initially in the equilibrium state
$\hat{\rho}_{tot}^{eq}$. The initial state is then modified by the first
magnetic pulses via the dipole operator as
$(\hat{M}^{\times}\hat{\rho}_{tot}^{eq})=[\hat{M},\hat{\rho}_{tot}^{eq}]$ at
$t=0$ and is propagated for time $t=s$ by $\mathcal{G}(s)$. The expectation
value is then obtained by calculating the trace of $\hat{M}$.

Actual 2D experiments have been conducted using a pair of magnetic pulses $a$
and $b$ with inter-pulse delay $\tau$. [21, 47] Under the impulsive
approximation, the magnetic field is expressed as

$B(s)=B_{a}\delta(s)+B_{b}\delta(s-\tau),$ (7)

where $B_{a}$ and $B_{b}$ are the magnetic field strengths. The nonlinear
element of the signal $M_{NL}(t,\tau)$ at $t+\tau$ is then evaluated as

$M_{NL}(t,\tau)=M_{ab}(t+\tau)-M_{a}(t+\tau)-M_{b}(t+\tau),$ (8)

where $M_{ab}(t+\tau)$, $M_{a}(t)$, and $M_{b}(t)$ are the total signal and
the linear elements for pulses $a$ and $b$, respectively.

From Eqs. (3) and (7), the nonlinear element is evaluated as

$\begin{split}M_{NL}(t,\tau)&=B_{a}B_{b}\,\chi_{2}(\tau,t)+B_{a}B_{b}^{2}\,\chi_{3}(\tau,0,t)\\\
&+B_{a}^{2}B_{b}\,\chi_{3}(0,\tau,t).\end{split}$ (9)

Thus, the characteristic features of THz-MR spectroscopy are described by the
nonlinear susceptibilities. Because each term on the rhs has different
proportionality with respect to $B_{a}$ and $B_{b}$, they can be evaluated
separately by changing their respective field strengths.

In nonlinear optical spectroscopies, $\chi_{2}$ is used to analyze 2D THz-
Raman [48, 49] and 2D THz-IR-Raman signals [50, 51] that involve a 2Q
transition process, while $\chi_{3}(s_{3},s_{2},s_{1})$ is used to analyze
nonlinear 2D THz spectroscopies [52, 53] including 2D THz rotational
spectroscopy. [54, 55] As illustrated in nonlinear optical spectroscopies, the
$s_{1}$ and $s_{3}$ periods describe the time evolutions of coherent states,
while $s_{2}$ describes the evolution of population states. [22, 23] While
$\chi_{3}(\tau,0,t)$ is used for 2D spectroscopies with waiting time
$s_{2}=0$, $\chi_{3}(0,\tau,t)$ is used for transient absorption spectra in
which the time evolution of the excited-state population is measured.

Using $\chi_{1}(s)$ in Eq. (4), the linear absorption spectrum is defined as

$\chi_{1}(\omega)={\rm{Im}}\int_{0}^{\infty}{{\rm{d}}}se^{-i\omega
s}\chi_{1}(s),$ (10)

where Im denotes the imaginary part. Following the experimental setup, we
consider two kinds of 2D spectrum for different time configurations of $\tau$
and $t$ in Eq. (6), expressed as: (1) $\chi_{3}(\tau,0,t)$ and (2)
$\chi_{3}(0,\tau,t)$. They are then expressed in the Fourier translation form
as

$\chi_{3}^{(1)}(\omega_{\tau},\omega_{t})={\rm{Im}}\int_{0}^{\infty}{{\rm{d}}}\tau\int_{o}^{\infty}{{\rm{d}}}t\,e^{-i\omega_{\tau}\tau-i\omega_{t}t}\chi_{3}(\tau,0,t),$
(11)
$\chi_{3}^{(2)}(\omega_{\tau},\omega_{t})={\rm{Im}}\int_{0}^{\infty}{{\rm{d}}}\tau\int_{o}^{\infty}{{\rm{d}}}t\,e^{-i\omega_{\tau}\tau-i\omega_{t}t}\chi_{3}(0,\tau,t),$
(12)

where $\omega_{\tau}$ and $\omega_{t}$ represent the excitation and detection
frequencies, respectively. [21]

##  III Chiral spin model and HEOM approach

To describe the novel phenomena of chiral magnets, we consider a linear chain
system consisting of $L$ spins. The system Hamiltonian is defined as

$\begin{split}\hat{H}_{S}&=J\sum_{n=1}^{L}\left[\hat{\sigma}_{n}^{z}\,\hat{\sigma}_{n+1}^{z}+\Delta\left(\hat{\sigma}_{n}^{x}\,\hat{\sigma}_{n+1}^{x}+\hat{\sigma}_{n}^{y}\,\hat{\sigma}_{n+1}^{y}\right)\right]\\\
&+D_{y}\sum_{n=1}^{L}\left[\hat{\sigma}_{n}^{z}\,\hat{\sigma}_{n+1}^{x}-\hat{\sigma}_{n}^{x}\,\hat{\sigma}_{n+1}^{z}\right],\end{split}$
(13)

where $\hat{\sigma}_{n}^{\alpha}$ $(\alpha=x,y,z)$ denotes the Pauli operator
at the $n$-th site in the $\alpha$-th direction. Here, the first term with
$J>0$ represents the AFM Heisenberg exchange coupling, and $\Delta$ is the
anisotropic parameter. The second term in Eq. (13) is the DMI that is
perpendicular to the $xz$ plane, with coupling strength $D_{y}$. The sign of
$D_{y}$ determines the direction (i.e., clockwise or anticlockwise) or
handedness (i.e., right- or left-handed circular configuration) of chirality,
which is not easily determined by experimental measurements. [3, 29] The
magnetization operator is defined as

$\hat{M}=\sum_{l=1}^{L}\hat{\sigma}_{l}^{x}.$ (14)

A very important aspect of investigating ultrafast dynamics in magnetic
materials is the inclusion of a quantum-mechanically consistent relaxation and
dephasing mechanism. This can be achieved by including a harmonic heat bath in
the spin system. The Hamiltonian of this heat bath is defined as

$\hat{H}_{B}=\sum_{j}\left[\frac{\hat{p}_{j}^{2}}{2m_{j}}+\frac{m_{j}\omega_{j}^{2}}{2}\hat{x}_{j}^{2}\right],$
(15)

where $\hat{p}_{j}$, $\hat{x}_{j}$, $\omega_{j}$, and $m_{j}$ represent the
momentum, position, frequency, and mass of the $j$-th oscillator,
respectively. The system–bath interaction is defined by

$\hat{H}_{I}=\hat{V}\,\sum_{j}g_{j}\,\hat{x}_{j},$ (16)

where $\hat{V}$ represents the spin part of the system–bath interaction
function, defined as

$\hat{V}=\sum_{n}^{L}\hat{\sigma}_{n}^{z},$ (17)

and $g_{j}$ is the coupling strength with the $j$-th oscillator. The total
Hamiltonian is then given by

$\hat{H}_{tot}=\hat{H}_{S}+\hat{H}_{I}+\hat{H}_{B}.$ (18)

In open quantum dynamics theory, the time-irreversible process can be
described using a reduced set of equations of motion. After reducing the
number of degrees of freedom of the heat bath, the noise effects are
characterized by the correlation function,

$C(t)=\frac{1}{\pi}\int_{0}^{\infty}{{\rm{d}}}\omega\,J(\omega)\left[\coth\left(\frac{\beta\hbar\omega}{2}\right)\cos(\omega
t)-i\sin(\omega t)\right],$ (19)

where $J(\omega)$ is the spectral density function (SDF), $\beta=1/k_{B}T$ is
the inverse temperature, and $k_{B}$ is the Boltzmann constant. For the heat
bath, we assume the Drude SDF, which is expressed as,

$J(\omega)=\frac{\zeta\gamma^{2}\omega}{\gamma^{2}+\omega^{2}},$ (20)

where $\zeta$ is the coupling strength and $\gamma$ is the inverse correlation
time.

In non-perturbative and non-Markovian conditions, the time evolution can be
described by the HEOM approach. [38, 39] We rewrite Eq. (19) in a linear-
summation form,

$C(t)=\sum_{k}^{K}c_{k}\,e^{-\nu_{k}|t|},$ (21)

where $c_{k}$ and $\nu_{k}$ are complex-valued coefficients. The HEOM can then
be expressed as [40, 41]

$\begin{split}\frac{\partial}{\partial
t}\hat{\rho}_{[\vec{n}]}(t)=&-\left[i\hat{H}_{S}^{\times}+\sum_{k}^{N}n_{k}\nu_{k}\right]\hat{\rho}_{[\vec{n}]}(t)\\\
&-i\hat{V}^{\times}\sum_{k}^{K}\hat{\rho}_{[\vec{n}+\vec{e}_{k}]}(t)\\\
&-i\sum_{k}^{K}\left[c_{k}\hat{V}\hat{\rho}_{[\vec{n}-\vec{e}_{k}]}(t)-c_{k}^{\ast}\hat{\rho}_{[\vec{n}-\vec{e}_{k}]}(t)\hat{V}\right],\end{split}$
(22)

where $\vec{n}=\\{n_{1},n_{2},\cdots,n_{K}\\}$ denotes the index vector, and
$n_{k}$ are non-negative integers. Among all the $\hat{\rho}_{[\vec{n}]}(t)$,
the $0$-th, with $\vec{0}=\\{0,0,\cdots,0\\}$, corresponds to the density
operator of the reduced system, while all the others are introduced for
ancillary purposes.

In the HEOM approach, the response functions Eqs. (4)–(6) are evaluated as the
time evolution of the system under external excitation. The density matrix is
replaced by the HEOM elements, and the Liouvillian $\mathcal{G}(t)$ is
replaced using Eq. (22).

For example, we can evaluate $\chi_{3}(\tau,0,t)$ using the expression of the
second line in Eq. (6) as follows. [40, 41] We first run the HEOM program for
a sufficiently long period from a temporally initial condition at $t=-t_{i}$
(such as the factorized initial condition,
$\hat{\rho}_{[\vec{0}]}(-t_{i})=\exp[-\beta\hat{H}_{S}]$ with all the other
hierarchy elements set to zero) to time $t=0$ to reach a true thermal
equilibrium, denoted by $\hat{\rho}_{[\vec{n}]}^{eq}(0)$. All the hierarchy
elements have to be used to define a correlated initial condition. The system
is next excited by the first magnetic interaction $\hat{M}$ at $t=0$ as
$\hat{\rho}_{[\vec{n}]}^{\prime}(0)=[\hat{M},\hat{\rho}_{[\vec{n}]}^{eq}(0)]$.
The perturbed elements $\hat{\rho}_{[\vec{n}]}^{\prime}$ then evolve in time
by numerically integrating Eq. (22) up to $t$. At $t$, the system is excited
by the second and third magnetic interactions as
$\hat{\rho}_{[\vec{n}]}^{\prime\prime}(t)=[\hat{M},[\hat{M},\hat{\rho}_{[\vec{n}]}^{\prime}(t)]]$.
Then, after $\hat{\rho}_{[\vec{n}]}^{\prime\prime}$ evolves in time with the
initial condition $\hat{\rho}_{[\vec{n}]}^{\prime\prime}(t)$ to $t+\tau$, the
response function is calculated from the expectation value of the magnetic
moment as
$\chi_{3}(\tau,0,t)=tr_{A}\\{\hat{M}\hat{\rho}_{[\vec{0}]}^{\prime\prime}(t+\tau)\\}$.

Notice that to take the system–bath coherence (or bath entanglement [41]) into
account during the external perturbation, it is important to operate $\hat{M}$
to all of the hierarchy elements $\hat{\rho}_{[\vec{n}]}^{\prime\prime}(t)$.
Although we only use $\hat{\rho}_{[\vec{0}]}(t)$ to calculate an expectation
value, the other elements are essential to obtain an echo signal for non-
Markovian noise in 2D spectroscopy. [40, 41]

##  IV Numerical Results

In the following, we set the exchange-coupling strength as the base unit $J=1$
and calculate 1D (or linear absorption) and 2D spectra for $\Delta=0.1$ as a
model for typical AFM metal oxides. [56, 57] The number of spin sites was
$L=6$, with periodic boundary conditions
$\hat{\sigma}_{n}^{\alpha}=\hat{\sigma}_{n+L}^{\alpha}$. Here, we consider
cases without the DMI ($D_{y}=0$) and with the weak DMI ($|D_{y}|\leq 0.05$),
which is appropriate for a small model system. Although the HEOM used in this
study can investigate cases in which the system interacts strongly with the
heat bath, here, we focus on the characterization of the chiral spin system
and keep the coupling weak. Thus, the bath parameters were chosen as
$\zeta=0.05$, $\gamma=1.0$, and $\beta\hbar=10$. The hierarchy parameters were
chosen as $K=50$ and $\sum_{k}^{K}n_{k}\leq 10$, and the time step for
numerical integration was set to $\delta t=0.01$.

###  IV.1 Energy eigenstates and transition elements of $\hat{M}$

Table 1:  Eigenenergies of spin Hamiltonian, $\hat{H}_{S}$ (a) without the DMI and (b) with the DMI.  | (a) $D_{y}=0$ | (b) $D_{y}=\pm 0.05$
---|---|---
$|0_{0}\rangle$ | 0 | 0
$|0_{1}\rangle$ | 0.001 | 0.001
$|1_{0}\rangle$ | $1.02$ | $0.99$
$|1_{0}^{\prime}\rangle$ | $-$ | $1.14$
$|1_{1}\rangle$ | $0.87$ | $0.81$
$|1_{1}^{\prime}\rangle$ | $-$ | $0.98$
$|1_{2}\rangle$ | $0.93$ | $0.93$
$|1_{3}\rangle$ | 1.15 | 1.15
$|2_{0}\rangle$ | 2.02 | 1.99
$|2_{1}\rangle$ | 2.03 | 2.03
$|2_{2}\rangle$ | 1.85 | 1.86
$|2_{3}\rangle$ | $2.02$ | $2.02$
$|2_{4}\rangle$ | $2.12$ | $2.12$
Figure 1: Schematic view of the energy states of the spin Hamiltonian,
$\hat{H}_{S}$, (a) without the DMI $(D_{y}=0)$ and (b) with the DMI
$(D_{y}=\pm 0.05)$. The black, blue, and red lines represent the states that
are (I) independent of the DMI, (II) dependent on the square of $D_{y}$, and
(III) independent of the DMI but the phases of their wave functions change
depending on the sign of $D_{y}$, respectively. The purple arrows denote the
possible transitions that arise from $\hat{M}$.

To illustrate the origin of the peaks in the 1D and 2D spectra that will be
shown later, in Table 1, we present some representative energy eigenstates
$|n_{m}^{(\prime)}\rangle$ that are necessary to explain THz-MR spectra. Here,
$n=0$, $1$, and $2$ represent the ground, first, and second excited states,
respectively, and $m$ and ′ are introduced to signify the energy splitting
arising from the anisotropic coupling $\Delta$ and the DMI, respectively.
Thus, for example, the states $|1_{0}\rangle$ and $|1_{0}^{\prime}\rangle$ are
degenerate when $D_{y}=0$. Note that to conduct numerical simulations, we
employed all of the spin-$z$ basis states to maintain the numerical accuracy.

In Fig. 1, we depict the energy eigenstates (a) without the DMI $(D_{y}=0)$
and (b) with the DMI $(D_{y}=\pm 0.05)$. The possible magnetic transitions
that arise from $\hat{M}$ are denoted by vertical arrows. Here, depending on
the role of $D_{y}$, we classify the eigenstates as (I) independent of the DMI
(black lines), (II) dependent on the square of $D_{y}$ (blue lines), and (III)
independent of the DMI but the phases of their wave functions change depending
on the sign of $D_{y}$ (red lines). For a finite $D_{y}$ in Fig. 1(b), the
degeneracy of eigenenergy is resolved because the inversion symmetry is
broken, and we observe the finite energy gap in (II) (blue dashed lines).

###  IV.2 1D THz-MR spectrum

Figure 2: Linear absorption (1D) spectra, $\chi_{1}(\omega)$ calculated for
$D_{y}=$ (a) $0$, (b) $+0.025$, and (c) $+0.05$. The signal intensities are
normalized with respect to the absolute values of the peak intensities.

In Fig. 2, linear absorption (1D) spectra $\chi_{1}(\omega)$ calculated from
Eq. (10) with Eq. (4) are presented for different $D_{y}$. Here, we only focus
on $D_{y}>0$ because the spectral profiles for $-|D_{y}|$ are identical to
those for $+|D_{y}|$. The tiny symmetrical peaks around the main peak arise as
an artifact of numerical Fourier transformation; these can be suppressed by
increasing the time interval and introducing window functions.

Under current low-temperature conditions, the spin system is almost in the
ground equilibrium states $|0_{0}\rangle$ and $|0_{1}\rangle$. The populations
of other excited states are less than $0.1\%$, and we cannot observe the
transitions from the excited states in the 1D spectrum. The main peak $A_{1}$
arises from the transition $|0_{0}\rangle\to|1_{0}\rangle$. Because the DMI
resolves the degeneracy of states $|1_{0}\rangle$ and $|1_{1}\rangle$ (blue
lines), the energy states $|1_{0}^{\prime}\rangle$ and
$|1_{1}^{\prime}\rangle$ (blue dashed lines) appear for finite DMI, as
illustrated in Fig. 1(b). Accordingly, we observe the adjoint peak $A_{2}$
that arises from the transition $|0_{0}\rangle\to|1_{0}^{\prime}\rangle$. The
energy eigenvalues of $|1_{0}\rangle$ and $|1_{0}^{\prime}\rangle$ are (II)
dependent upon the square of $D_{y}$, and we found that, for our system with
$\Delta/J=0.1$, these peak positions are evaluated as
$\omega_{A_{1}}\approx-8.34\,D_{y}^{2}+1.02$ and $\omega_{A_{2}}\approx
8.34\,D_{y}^{2}+1.12$. Thus, from the position of the $A_{2}$ peak, we can
estimate the magnitude of $D_{y}$, but we cannot determine the direction of
chirality, which is described as the sign of $D_{y}$.

###  IV.3 2D THz-MR spectra

To elucidate the chirality of the system further, we next present numerical
results for 2D THz-MR spectrum evaluated from
$\chi_{3}^{(1)}(\omega_{\tau},\omega_{t})$ and
$\chi_{3}^{(2)}(\omega_{\tau},\omega_{t})$, expressed as Eqs. (11) and (12)
with Eq. (6). In the present case, the rephasing parts ($\omega_{t}<0$ and
$\omega_{\tau}>0$) in both 2D spectra have no signal, indicating that there is
no spin echo. [22, 23] Thus, we only present the non-rephasing parts of the
spectra ($\omega_{t}>0$ and $\omega_{\tau}>0$).

####  IV.3.1 Contribution from $\chi_{3}^{(1)}(\omega_{\tau},\omega_{t})$

Figure 3: 2D THz-MR spectrum calculated from
$\chi_{3}^{(1)}(\omega_{\tau},\omega_{t})$ for three values of $D_{y}$: (a)
$0$, (b) $+0.05$, and (c) $-0.05$. The intensity of each spectrum is
normalized with respect to its maximum peak intensity. The red and blue areas
represent the positive and negative intensities of the spectra.  Figure 4:
Double-sided Feynman diagrams for all labeled peaks. Two different
contributions of $B_{3}$ are shown separately.

In Fig. 3, we present contour maps of
$\chi_{3}^{(1)}(\omega_{\tau},\omega_{t})$ in the (a) $D_{y}=0$, (b) $+0.05$,
and (c) $-0.05$ cases, in which the red and blue areas represent positive and
negative intensities. To analyze the physical origin of each peak, we present
a double-sided Feynman diagram corresponding to each labeled peak in Fig.
4(a). Here, the $B_{3}$ peak consists of two Liouville paths, which we denote
as $B_{3}^{(i)}$ and $B_{3}^{(ii)}$.

These diagrams describe the time evolution of the left and right wave
functions of the density operators involved in the response function from
bottom to top under the magnetic excitations and deexcitations. For example,
diagram $B_{3}^{(i)}$ in Fig. 4(a) describes the time evolution of the density
operator undergoing magnetic interactions described by the purple arrows in
Fig. 1 with the left wave function at time $s=0$, $s=\tau$, and $s=\tau+t$
with the right wave function at time $s=\tau$. The left wave function evolves
with time
$|0_{1}\rangle\rightarrow|2_{1}\rangle\rightarrow|2_{4}\rangle\rightarrow|1_{3}\rangle$,
while the right one evolves with time $\langle 0_{1}|\rightarrow\langle
1_{3}|$. [6, 41]

The peaks labeled as $B_{1}$ and $B_{2}$ along $\omega_{t}=0.98$ correspond to
peaks $A_{1}$ and $A_{2}$ in the 1D spectra, and the frequency difference
$\omega_{A_{2}}-\omega_{A_{1}}$ is equal to $\omega_{B_{2}}-\omega_{B_{1}}$.
As illustrated in Fig. 4(a), the diagrams $B_{1}$ and $B_{2}$ only involve the
eigenstates that are (I) independent of the DMI and (II) dependent on the
square of $D_{y}$. Thus, the intensities of the $B_{1}$ and $B_{2}$ peaks do
not depend on the sign of $D_{y}$.

However, the $B_{3}^{(ii)}$ diagram in Fig. 4(a) involves states (II) and
(III), while the $B_{3}^{(i)}$ diagram contains the state $|2_{4}\rangle$ that
is (III) independent of the DMI but the phase of its wave function changes
depending on the sign of $D_{y}$. As a result, the sign of the peak intensity
changes for $D_{y}=+0.05$ and $-0.05$. However, the peak profiles in Figs.
3(b2) and 3(c2) are not positively and negatively symmetric due to the
contribution from $B_{3}^{(ii)}$.

####  IV.3.2 Contribution from $\chi_{3}^{(2)}(\omega_{\tau},\omega_{t})$

Figure 5: 2D THz-MR spectrum calculated from
$\chi_{3}^{(2)}(\omega_{\tau},\omega_{t})$ for three values of $D_{y}$: (a)
$0$, (b) $+0.05$, and (c) $-0.05$. The intensity of each spectrum is
normalized with respect to its maximum peak intensity. The red and blue areas
represent positive and negative intensities. All the parameter values are the
same as those in Fig. 3.

We next depict $\chi_{3}^{(2)}(\omega_{\tau},\omega_{t})$ in Fig. 5. Double-
sided Feynman diagrams corresponding to each labeled peak are presented in
Figs. 4(b) and 4(c).

Peaks $C_{1}$–$C_{4}$ appear in the case $D_{y}\neq 0$ because they involve
the states classified in (II), as depicted in Fig. 4. Moreover, because peaks
$C_{1}$–$C_{4}$ involve the $|1_{2}\rangle$ and $|2_{3}\rangle$ states
classified in (III) in each diagram in Fig. 4(b), the sign of the peak
intensity changes in correspondence with the sign of $D_{y}$. Note that the
peak profiles of $C_{1}$ for $D_{y}=0$ and $D_{y}>0$ are identical because the
$|1_{2}\rangle$ state does not vanish even for $D_{y}=0$. We also find that
the positive peak near $(\omega_{\tau},\omega_{t})=(1.05,1)$ in Fig. 5(a1)
arises as the consequence of the heat-bath-induced coherence
$|0_{0}\rangle\langle 1_{0}|$.

In Figs. 4(b2) and 4(c2), we observe that peaks $C_{5}$–$C_{9}$ arise from the
2Q coherence. Although the $C_{9}$ diagram involves the $|2_{3}\rangle$ state
in (III), it overlaps with $C_{7}$ for $D_{y}\neq 0$, and the signs of the
peak intensities cannot be easily evaluated. Note that here we chose a small
system ($L=6$) with periodic boundary conditions, so energy states higher than
the second excited state are lowered, and the profiles of peaks
$C_{5}$–$C_{9}$ are distorted. Hence, within the calculations of the present
model, it is more reliable to use the profiles of $C_{1}$–$C_{4}$ to determine
the sign of $D_{y}$.

##  V Conclusion

A recently developed 2D THz-MR spectroscopic technique has created new
possibilities for measuring complex molecular magnetic systems. In the present
work, we have illustrated the key features of this technique based on
nonlinear response-function theory and have described a method for simulating
2D THz-MR spectra through the use of the HEOM formalism. Using simulated 1D
and 2D THz-MR spectra for a chiral magnetic material, we demonstrated that the
2D technique allows us to evaluate the pitch and direction of chirality
determined from the strength and sign of the DMI, while only the absolute
amplitude of the DMI can be determined from 1D measurements.

The reason the sign of the DMI was detected by 2D spectroscopy in this study
is that there are eigenstates in which the phase of the wavefunction changes
following the sign of the DMI, which is categorized as (III). Although there
have been studies of eigenstates of spin systems with the DMI, [58] the
existence of eigenstates that change phase with the sign of the DMI, as found
in this study, has not been explored. It is important to know the causes of
these eigenstates because they may lead to the appearance of novel phenomena.

To make a direct comparison between the results of our simulations and those
obtained experimentally, however, we must increase the number of spins in
accordance with those available in experimental systems. In investigating spin
dynamics in the condensed phase, it is also important to consider the non-
perturbative and non-Markovian system–bath interactions. Nevertheless, we
believe that the present results elucidate the key features of 2D THz-MR
spectroscopic methods with regard to probing the fundamental nature of a
magnetic spin system.

For further investigations to monitor the ultrafast dynamical aspects of the
spin system, such as the dynamics of spin waves, it is necessary to conduct a
variety of advanced nonlinear spectroscopic approaches, such as pump-probe and
transient absorption measurements to foster the development of this
spectroscopic method. We leave such extensions to future studies, depending on
progress in experimental and simulation techniques.

## Acknowledgments

The authors are thankful to Professor Jun-ichiro Kishine and Professor Keisuke
Tominaga for helpful discussions. Y.T. was supported by JSPS KAKENHI (Grant
No. B21H01884). J.Z. was supported by JST SPRING (Grant No. JPMJSP2110).

### Conflict of Interest

The authors have no conflicts to disclose.

## Data availability

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

## REFERENCES

## References

  * Schraml and Bellama [1988] J. Schraml and J. M. Bellama,  _Two-dimensional NMR spectroscopy_ (John Wiley & Sons, 1988).
  * Levitt [2013] M. H. Levitt,  _Spin dynamics: basics of nuclear magnetic resonance_ (John Wiley & Sons, 2013).
  * Sato  _et al._ [2016] T. J. Sato, D. Okuyama, T. Hong, A. Kikkawa, Y. Taguchi, T. Arima, and Y. Tokura, “Magnon dispersion shift in the induced ferromagnetic phase of noncentrosymmetric mnsi,” Phys. Rev. B 94, 144420 (2016).
  * Neu and Schmuttenmaer [2018] J. Neu and C. A. Schmuttenmaer, “Tutorial: An introduction to terahertz time domain spectroscopy (thz-tds),” Journal of Applied Physics 124, 231101 (2018), https://doi.org/10.1063/1.5047659 .
  * Zhang  _et al._ [2023] Z. Zhang, F. Sekiguchi, T. Moriyama, S. C. Furuya, M. Sato, T. Satoh, Y. Mukai, K. Tanaka, T. Yamamoto, H. Kageyama, Y. Kanemitsu, and H. Hirori, “Generation of third-harmonic spin oscillation from strong spin precession induced by terahertz magnetic near fields,” Nature Communications 14, 1795 (2023).
  * Mukamel [1995] S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_, Oxford series in optical and imaging sciences (Oxford University Press, 1995).
  * Herrmann [1963] G. Herrmann, “Resonance and high frequency susceptibility in canted antiferromagnetic substances,” Journal of Physics and Chemistry of Solids 24, 597–606 (1963).
  * Nakajima  _et al._ [2010] M. Nakajima, A. Namai, S. Ohkoshi, and T. Suemoto, “Ultrafast time domain demonstration of bulk magnetization precession at zero magnetic field ferromagnetic resonance induced by terahertz magnetic field,” Opt. Express 18, 18260–18268 (2010).
  * Yamaguchi, Nakajima, and Suemoto [2010] K. Yamaguchi, M. Nakajima, and T. Suemoto, “Coherent control of spin precession motion with impulsive magnetic fields of half-cycle terahertz radiation,” Phys. Rev. Lett. 105, 237201 (2010).
  * Kampfrath  _et al._ [2011] T. Kampfrath, A. Sell, G. Klatt, A. Pashkin, S. Mährlein, T. Dekorsy, M. Wolf, M. Fiebig, A. Leitenstorfer, and R. Huber, “Coherent terahertz control of antiferromagnetic spin waves,” Nature Photonics 5, 31 (2011).
  * Satoh  _et al._ [2010] T. Satoh, S.-J. Cho, R. Iida, T. Shimura, K. Kuroda, H. Ueda, Y. Ueda, B. A. Ivanov, F. Nori, and M. Fiebig, “Spin oscillations in antiferromagnetic nio triggered by circularly polarized light,” Phys. Rev. Lett. 105, 077402 (2010).
  * Turgut  _et al._ [2017] E. Turgut, A. Park, K. Nguyen, A. Moehle, D. A. Muller, and G. D. Fuchs, “Chiral magnetic excitations in fege films,” Phys. Rev. B 95, 134416 (2017).
  * Dou  _et al._ [2023] K. Dou, W. Du, Z. He, Y. Dai, B. Huang, and Y. Ma, “Theoretical prediction of antiferromagnetic skyrmion crystal in janus monolayer CrSi2N2As2,” ACS Nano 17, 1144–1152 (2023), pMID: 36448916, https://doi.org/10.1021/acsnano.2c08544 .
  * Soumyanarayanan  _et al._ [2017] A. Soumyanarayanan, M. Raju, A. L. Gonzalez Oyarce, A. K. C. Tan, M.-Y. Im, A. P. Petrovic, P. Ho, K. H. Khoo, M. Tran, C. K. Gan, F. Ernult, and C. Panagopoulos, “Tunable room-temperature magnetic skyrmions in ir/fe/co/pt multilayers,” Nature Materials 16, 898–904 (2017).
  * Ozerov  _et al._ [2014] M. Ozerov, J. Romhányi, M. Belesi, H. Berger, J.-P. Ansermet, J. van den Brink, J. Wosnitza, S. A. Zvyagin, and I. Rousochatzakis, “Establishing the fundamental magnetic interactions in the chiral skyrmionic mott insulator ${\mathrm{cu}}_{2}{\mathrm{oseo}}_{3}$ by terahertz electron spin resonance,” Phys. Rev. Lett. 113, 157205 (2014).
  * Hayashi  _et al._ [2021] Y. Hayashi, Y. Okamura, N. Kanazawa, T. Yu, T. Koretsune, R. Arita, A. Tsukazaki, M. Ichikawa, M. Kawasaki, Y. Tokura, and Y. Takahashi, “Magneto-optical spectroscopy on Weyl nodes for anomalous and topological hall effects in chiral MnGe,” Nature Communications 12, 5974 (2021).
  * Tsuruta  _et al._ [2018] K. Tsuruta, M. Mito, H. Deguchi, J. Kishine, Y. Kousaka, J. Akimitsu, and K. Inoue, “Nonlinear magnetic responses at the phase boundaries around helimagnetic and skyrmion lattice phases in MnSi: Evaluation of robustness of noncollinear spin texture,” Phys. Rev. B 97, 094411 (2018).
  * Biesner  _et al._ [2022] T. Biesner, S. Roh, A. Pustogow, H. Zheng, J. F. Mitchell, and M. Dressel, “Magnetic terahertz resonances above the néel temperature in the frustrated kagome antiferromagnet averievite,” Phys. Rev. B 105, L060410 (2022).
  * Peedu  _et al._ [2022] L. Peedu, V. Kocsis, D. Szaller, B. Forrai, S. Bordács, I. Kézsmárki, J. Viirok, U. Nagel, B. Bernáth, D. L. Kamenskyi, A. Miyata, O. Portugall, Y. Tokunaga, Y. Tokura, Y. Taguchi, and T. Rõ om, “Terahertz spectroscopy of spin excitations in magnetoelectric ${\mathrm{lifepo}}_{4}$ in high magnetic fields,” Phys. Rev. B 106, 134413 (2022).
  * Kozuki, Nagashima, and Hangyo [2011] K. Kozuki, T. Nagashima, and M. Hangyo, “Measurement of electron paramagnetic resonance using terahertz time-domain spectroscopy,” Opt. Express 19, 24950–24956 (2011).
  * Lu  _et al._ [2017] J. Lu, X. Li, H. Y. Hwang, B. K. Ofori-Okai, T. Kurihara, T. Suemoto, and K. A. Nelson, “Coherent two-dimensional terahertz magnetic resonance spectroscopy of collective spin waves,” Phys. Rev. Lett. 118, 207204 (2017).
  * Hamm and Zanni [2011] P. Hamm and M. T. Zanni,  _Concepts and Methods of 2D Infrared Spectroscopy_ (Cambridge University Press, 2011).
  * Cho [2009] M. Cho,  _Two-Dimensional Optical Spectroscopy_ (CRC Press, 2009).
  * Fert, Reyren, and Cros [2017] A. Fert, N. Reyren, and V. Cros, “Magnetic skyrmions: advances in physics and potential applications,” Nature Reviews Materials 2, 17031 (2017).
  * Jungwirth  _et al._ [2016] T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, “Antiferromagnetic spintronics,” Nature Nanotechnology 11, 231–241 (2016).
  * Inoue [2021] K. Inoue, “Chiral magnetism: Coupling static and dynamic chirality,” Chemistry Letters 50, 742–751 (2021), https://doi.org/10.1246/cl.200840 .
  * Togawa  _et al._ [2016] Y. Togawa, Y. Kousaka, K. Inoue, and J.-i. Kishine, “Symmetry, structure, and dynamics of monoaxial chiral magnets,” Journal of the Physical Society of Japan 85, 112001 (2016), https://doi.org/10.7566/JPSJ.85.112001 .
  * Kato, Yamamoto, and Kishine [2022] A. Kato, H. M. Yamamoto, and J.-i. Kishine, “Chirality-induced spin filtering in pseudo jahn-teller molecules,” Phys. Rev. B 105, 195117 (2022).
  * Kishine, Inoue, and Yoshida [2005] J.-i. Kishine, K. Inoue, and Y. Yoshida, “Synthesis, Structure and Magnetic Properties of Chiral Molecule-Based Magnets,” Progress of Theoretical Physics Supplement 159, 82–95 (2005), https://academic.oup.com/ptps/article-pdf/doi/10.1143/PTPS.159.82/5294973/159-82.pdf .
  * Lodahl  _et al._ [2017] P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeutel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller, “Chiral quantum optics,” Nature 541, 473–480 (2017).
  * Mito  _et al._ [2009] M. Mito, K. Iriguchi, H. Deguchi, J.-i. Kishine, K. Kikuchi, H. Ohsumi, Y. Yoshida, and K. Inoue, “Giant nonlinear magnetic response in a molecule-based magnet,” Phys. Rev. B 79, 012406 (2009).
  * Okamura  _et al._ [2015] Y. Okamura, F. Kagawa, S. Seki, M. Kubota, M. Kawasaki, and Y. Tokura, “Microwave magnetochiral dichroism in the chiral-lattice magnet ${\mathrm{cu}}_{2}{\mathrm{oseo}}_{3}$,” Phys. Rev. Lett. 114, 197202 (2015).
  * Clements  _et al._ [2017] E. M. Clements, R. Das, L. Li, P. J. Lampen-Kelley, M.-H. Phan, V. Keppens, D. Mandrus, and H. Srikanth, “Critical behavior and macroscopic phase diagram of the monoaxial chiral helimagnet cr1/3nbs2,” Scientific Reports 7, 6545 (2017).
  * Nomura  _et al._ [2019] T. Nomura, X.-X. Zhang, S. Zherlitsyn, J. Wosnitza, Y. Tokura, N. Nagaosa, and S. Seki, “Phonon magnetochiral effect,” Phys. Rev. Lett. 122, 145901 (2019).
  * Dijkstra and Tanimura [2010] A. G. Dijkstra and Y. Tanimura, “Non-Markovian entanglement dynamics in the presence of system-bath coherence,” Phys. Rev. Lett. 104, 250401 (2010).
  * Weiss [2012] U. Weiss,  _Quantum Dissipative Systems_, 4th ed. (WORLD SCIENTIFIC, 2012) https://www.worldscientific.com/doi/pdf/10.1142/8334 .
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath,” Journal of the Physical Society of Japan 58, 101–114 (1989).
  * Ishizaki and Tanimura [2005] A. Ishizaki and Y. Tanimura, “Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach,” Journal of the Physical Society of Japan 74, 3131–3134 (2005).
  * Tanimura [2006] Y. Tanimura, “Stochastic Liouville, Langevin, Fokker-Planck, and master equation qpproaches to quantum dissipative systems,” Journal of the Physical Society of Japan 75, 082001 (2006).
  * Tanimura [2020] Y. Tanimura, “Numerically ”exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM),” The Journal of Chemical Physics 153, 020901 (2020).
  * Bloch [1946] F. Bloch, “Nuclear induction,” Physical Review 70, 460–474 (1946).
  * Redfield [1965] A. Redfield, “The theory of relaxation processes,”  _Advances in Magnetic Resonance_ ,  Advances in Magnetic and Optical Resonance, 1, 1–32 (1965).
  * Kubo [1969] R. Kubo,  _Advances in Chemical Physics_ (John Wiley and Sons, Ltd, 1969) pp. 101–127, https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470143605.ch6 .
  * Joutsuka and Tanimura [2008] T. Joutsuka and Y. Tanimura, “Detecting the Dzyaloshinskii-Moriya interaction by means of pulsed EPR spectroscopy,” Chemical Physics Letters 457, 237–240 (2008).
  * Takahashi and Tanimura [2020] H. Takahashi and Y. Tanimura, “Open quantum dynamics theory of spin relaxation: Application to $\mu$SR and low-field NMR spectroscopies,” Journal of the Physical Society of Japan 89, 064710 (2020).
  * Wan and Armitage [2019] Y. Wan and N. P. Armitage, “Resolving continua of fractional excitations by spinon echo in thz 2d coherent spectroscopy,” Phys. Rev. Lett. 122, 257401 (2019).
  * Hamm and Shalit [2017] P. Hamm and A. Shalit, “Perspective: Echoes in 2D-Raman-THz spectroscopy,” The Journal of Chemical Physics 146, 130901 (2017), https://doi.org/10.1063/1.4979288 .
  * Ikeda, Ito, and Tanimura [2015] T. Ikeda, H. Ito, and Y. Tanimura, “Analysis of 2D THz-Raman spectroscopy using a non-Markovian Brownian oscillator model with nonlinear system-bath interactions,” The Journal of Chemical Physics 142, 212421 (2015).
  * Perakis  _et al._ [2016] F. Perakis, L. De Marco, A. Shalit, F. Tang, Z. R. Kann, T. D. Kühne, R. Torre, M. Bonn, and Y. Nagata, “Vibrational spectroscopy and dynamics of water,” Chemical Reviews 116, 7590–7607 (2016), https://doi.org/10.1021/acs.chemrev.5b00640 .
  * Takahashi and Tanimura [2023] H. Takahashi and Y. Tanimura, “Discretized hierarchical equations of motion in mixed liouville–wigner space for two-dimensional vibrational spectroscopies of liquid water,” The Journal of Chemical Physics 158, 044115 (2023), https://doi.org/10.1063/5.0135725 .
  * Reimann, Woerner, and Elsaesser [2021] K. Reimann, M. Woerner, and T. Elsaesser, “Two-dimensional terahertz spectroscopy of condensed-phase molecular systems,” The Journal of Chemical Physics 154, 120901 (2021), https://doi.org/10.1063/5.0046664 .
  * Okumura and Tanimura [1998] K. Okumura and Y. Tanimura, “Two-dimensional THz spectroscopy of liquids: non-linear vibrational response to a series of thz laser pulses,” Chemical Physics Letters 295, 298–304 (1998).
  * Lu  _et al._ [2016] J. Lu, Y. Zhang, H. Y. Hwang, B. K. Ofori-Okai, S. Fleischer, and K. A. Nelson, “Nonlinear two-dimensional terahertz photon echo and rotational spectroscopy in the gas phase,” Proceedings of the National Academy of Sciences 113, 11800–11805 (2016), https://www.pnas.org/doi/pdf/10.1073/pnas.1609558113 .
  * Iwamoto and Tanimura [2019] Y. Iwamoto and Y. Tanimura, “Open quantum dynamics of a three-dimensional rotor calculated using a rotationally invariant system-bath Hamiltonian: Linear and two-dimensional rotational spectra,” The Journal of Chemical Physics 151, 044105 (2019).
  * Braun  _et al._ [2005] H.-B. Braun, J. Kulda, B. Roessli, D. Visser, K. W. Krämer, H.-U. Güdel, and P. Böni, “Emergence of soliton chirality in a quantum antiferromagnet,” Nature Physics 1, 159–163 (2005).
  * Thoma  _et al._ [2022] H. Thoma, V. Hutanu, R. Dutta, A. Gukasov, V. Kocsis, Y. Tokunaga, Y. Taguchi, Y. Tokura, I. Kézsmárki, G. Roth, and M. Angst, “Magnetic order and sign of the Dzyaloshinskii-Moriya interaction in 2-D antiferromagnet Ba2CoGe2o7 under applied magnetic field,” IEEE Transactions on Magnetics 58, 1–5 (2022).
  * Sakai, Cépas, and Ziman [2000] T. Sakai, O. Cépas, and T. Ziman, “An electron spin resonance selection rule for spin-gapped systems,” Journal of the Physical Society of Japan 69, 3521–3524 (2000), https://doi.org/10.1143/JPSJ.69.3521 .


# Wall-anchored semiflexible polymer under large amplitude oscillatory shear
flow

Antonio Lamura  Corresponding author<EMAIL_ADDRESS>Istituto
Applicazioni Calcolo, CNR, Via Amendola 122/D, 70126 Bari, Italy  Roland G.
Winkler<EMAIL_ADDRESS>Theoretical Physics of Living Matter,
Institute for Advanced Simulation and Institute of Biological Information
Processing, Forschungszentrum Jülich, 52425 Jülich, Germany  Gerhard Gompper
<EMAIL_ADDRESS>Theoretical Physics of Living Matter, Institute for
Advanced Simulation and Institute of Biological Information Processing,
Forschungszentrum Jülich, 52425 Jülich, Germany

###### Abstract

The properties of semiflexible polymers tethered by one end to an impenetrable
wall and exposed to oscillatory shear flow are investigated by mesoscale
simulations. A polymer, confined in two dimensions, is described by a linear
bead-spring chain, and fluid interactions are incorporated by the Brownian
multiparticle collision dynamics approach. At small strains, the polymers
follow the applied flow field. However, at high strain, we find a strongly
nonlinear response, with major conformational changes. Polymers are stretched
along the flow direction and exhibit U-shaped conformations while following
the flow. As a consequence of confinement in the half-space, a frequency
doubling in the time-dependent polymer properties appears along the direction
normal to the wall.

##  I Introduction

Studies of the nonequilibrium behavior of polymer systems provide a link
between their microscopic molecular characteristics and the emerging
macroscopic dynamical properties Bird, Armstrong, and Hassager (1987); Larson
(1999); Rubinstein and Colby (2003); Larson and Desai (2015); Shaw and
MacKnight (2018); Schroeder (2018). A major experimental breakthrough in
resolving and visualizing the nonequilibrium properties of individual
molecules was achieved by studies on single DNA filaments Perkins  _et al._
(1995), which paved the way to explore a large variety of polymer properties
under flows.

From the theoretical side, several detailed computational molecular models
have been considered, which enhance our understanding of molecular processes.
In particular, single-polymer studies provide the opportunity to directly
observe the microscopic conformational properties of individual polymers close
to equilibrium or under flow conditions, thereby facilitating access to their
non-equilibrium conformational properties, which are ultimately responsible
for the macroscopic rheological behavior. The desire to visualize individual
polymer conformations in flow, both from an experimental and simulation point
of view, is strongly linked with advances in the statistical description of
their properties provided by molecular theoretical models Schroeder (2018);
Winkler (2006a, 2010).

Semiflexible polymer-type molecular structures are omnipresent in biological
systems, e.g., DNA, actin filaments, and microtubules, and their rigidity is
fundamental for their functions. Indeed, actin filaments contribute with their
rigidity to the mechanical properties of the cytoskeleton and the ability of
DNA to pack in the genome or inside a virus capsid is controlled by its
persistence length. As a consequence, the properties of semiflexible polymers
have been intensively investigated Wilhelm and Frey (1996); Götter  _et al._
(1996); Harnau, Winkler, and Reineker (1996); Everaers  _et al._ (1999);
Samuel and Sinha (2002); Le Goff  _et al._ (2002); Winkler (2003); Petrov  _et
al._ (2006). Here, theoretical Bird  _et al._ (1987); Öttinger (1996); Winkler
(2006b); Munk  _et al._ (2006); Winkler (2010) and computational Hur and
Shaqfeh (2000); Jendrejack, de Pablo, and Graham (2002); Hsieh and Larson
(2004); Liu, Ashok, and Muthukumar (2004); Celani, Puliafito, and Turitsyn
(2005); Ryder and Yeomans (2006); Ripoll, Winkler, and Gompper (2006); Sendner
and Netz (2008); He  _et al._ (2009); Zhang  _et al._ (2009); Huang  _et al._
(2010, 2011); Huang, Gompper, and Winkler (2012); Lang, Obermayer, and Frey
(2014); Nikoubashman and Howard (2017); Kong  _et al._ (2019); Romo-Uribe
(2021); Shee  _et al._ (2021); Nikoubashman (2021) studies revealed novel
dynamical, conformational, and rheological properties.

The large majority of these studies focuses on three-dimensional systems, and
predominantly single-polymer theoretical and simulation studies have been
performed under steady shear and extensional flow in the bulk Shaqfeh (2005);
Winkler (2010); Huang  _et al._ (2011). Similarly, in experiments single DNA
filaments have been considered, in particular filaments grafted on a wall by
one end under steady shear, and their mechanical properties Smith, Finzi, and
Bustamante (1992) and relaxation dynamics Perkins  _et al._ (1994) have been
examined. Further experiments Doyle, Ladoux, and Viovy (2000); Lueth and
Shaqfeh (2009) revealed the existence of a cyclic dynamics, which has been
confirmed by simulations Lueth and Shaqfeh (2009); Delgado-Buscalioni (2006).

Polymers in two dimensions are of interest on their on. Compared to the three-
dimensional case, there are two major differences. On the one hand, excluded-
volume interactions are more important and, on the other hand, hydrodynamics
can be neglected, specifically, in the case of strongly adsorbed polymers
Maier, Seifert, and Rädler (2002). Such systems can be realized experimentally
by considering filaments strongly adsorbed onto a surface, a membrane, or at
interfaces separating immiscible fluids Maier, Seifert, and Rädler (2002);
Cherstvy and Petrov (2014); Herold, Schwille, and Petrov (2010). In
particular, end-anchored semiflexible polymers have been considered, where the
central monomer is subject to an oscillatory force, and a transition from a
limit cycle to an aperiodic dynamics with increasing rigidity has been
predicted theoretically Chattopadhyay and Marenduzzo (2007). The dynamics of
semiflexible polymers under shear has been investigated numerically and
substantial stiffness-dependent conformational changes at high flow rates have
been obtained Lamura and Winkler (2012). Specifically, more flexible polymers
are extended by the flow, while stiffer filaments contract. Moreover,
filaments are aligned by the flow with a tumbling motion that, at high shear
rates, resembles the motion of flexible polymers in three dimensions.

Theory, simulations, and experiments typically consider stationary flows in
the linear viscoelastic regime. However, it is important to unravel the
polymer dynamics when time-dependent flows are applied, specifically large
amplitude oscillatory shear (LAOS) flows. Such studies provide additional
insight into the macroscopic polymer properties, specifically their
viscoelastic behavior Hyun  _et al._ (2011); Rogers (2018). Despite the deep
interest and experimental relevance, so far very few theoretical and
simulation studies have been performed under such flow conditions. In Ref.
Chen  _et al._ (2005), the extension and migration of a chain, confined in a
microchannel and subject to oscillatory pressure-driven flow, were observed.
Later, Brownian dynamics studies of macromolecules under oscillatory shear
flow have been performed Thomas, DePuit, and Khomami (2009); Lamura and
Winkler (2019). In contrast, many more LAOS experiments have been conducted
Hyun  _et al._ (2011); Zhou and Schroeder (2016); Rogers (2018). Yet, no study
so far considered a single tethered semiflexible chain subject to oscillatory
shear flow.

In this paper, we investigate the conformational and dynamical properties of a
semiflexible polymer in two dimensions, tethered by one of its ends on a
reflecting wall, in the presence of an oscillatory shear flow by mesoscale
simulations. A polymer is modeled as a self-avoiding worm-like filament.
Neglecting hydrodynamic interactions, the polymer is assumed to be in contact
with a Brownian heat bath implemented via the Brownian version Ripoll,
Winkler, and Gompper (2007); Gompper  _et al._ (2009) of the multiparticle
collision dynamics approach Kapral (2008); Gompper  _et al._ (2009). The
polymer properties are characterized in terms of the strain defined as the
ratio of the shear rate to the shear frequency. For low strain, the polymers
behave roughly as at equilibrium. In the opposite limit of high strain, they
show in each shear half-cycle average properties comparable to those at steady
shear—they fully elongate and flip back and forth following the external flow.
For intermediate values of strain, flipping is disfavored, and polymers may
remain preferentially in a coiled state when shear changes sign without re-
orientating along the instantaneous flow direction. The analysis of the
amplitudes of the undulation modes suggests that polymers exhibit a flexible
polymer-like behavior at larger length scales along the shear direction,
despite their significant stiffness in absence of the external flow. The
periodic motion of the center-of-mass displays the same frequency as the
external shear along the flow direction, while frequency doubling in the
normal direction appears by wall reflection, which “repels” the polymer as
flow drags it from one side to the other in a cycle.

The numerical models for the polymer and the Brownian fluid are introduced in
Sec. II. The results for the conformational and the dynamical behavior are
presented in Sec. III. Finally, Sec. IV summarizes and discusses the main
findings of this study.

##  II Model and Method

We model a filament as a linear bead-spring chain with $N$ beads of mass $M$
separated by bonds of length $r_{0}$, and confined in the positive half-plane
of the two dimensional space. The first bead is tethered at the origin
$(0,0)^{T}$ of the Cartesian reference frame, with no preferred orientation of
the first bond. The beads are subjected to forces by the total potential
$U=U_{bond}+U_{bend}+U_{ex}+U_{w}$. Bonds are described by the harmonic
potential

$U_{bond}=\frac{\kappa_{h}}{2}\sum_{i=1}^{N-1}(|{\bm{r}}_{i+1}-{\bm{r}}_{i}|-r_{0})^{2},$
(1)

where ${\bm{r}}_{i}$ denotes the position vector of the $i-$th bead
($i=1,\ldots,N$) and $\kappa_{h}$ is the force constant. The stiffness of the
polymer is implemented by the bending potential

$U_{bend}=\kappa\sum_{i=1}^{N-2}(1-\cos\varphi_{i}),$ (2)

where $\kappa$ is the bending rigidity and $\varphi_{i}$ is the angle between
two consecutive bond vectors. The filament persistence length $L_{p}$ is
related to $\kappa$ via $L_{p}=2\kappa r_{0}/k_{B}T$, where $k_{B}T$ is the
thermal energy, with $T$ the temperature and $k_{B}$ Boltzmann’s constant. The
shifted and truncated Lennard-Jones potential

$U_{ex}=4\epsilon\Big{[}\Big{(}\frac{\sigma}{r}\Big{)}^{12}-\Big{(}\frac{\sigma}{r}\Big{)}^{6}+\frac{1}{4}\Big{]}\Theta(2^{1/6}\sigma-r)$
(3)

ensures the self-avoidance of non-connected beads. Here, $r$ is the distance
between two beads and $\Theta(r)$ is the Heaviside function ($\Theta(r)=0$ for
$r<0$ and $\Theta(r)=1$ for $r\geq 0$). Confinement in the half-plane $y>0$ is
achieved by a repulsive wall implemented via the potential

$U_{w}=4\epsilon\Big{[}\Big{(}\frac{\sigma_{w}}{y}\Big{)}^{12}-\Big{(}\frac{\sigma_{w}}{y}\Big{)}^{6}+\frac{1}{4}\Big{]}\Theta(2^{1/6}\sigma_{w}-y),$
(4)

where $y$ is the distance of a bead from the wall. The dynamics of the beads
is described by Newton’s equations of motion, which are integrated by the
velocity-Verlet algorithm with time step $\Delta t_{p}$ Swope  _et al._
(1982); Allen and Tildesley (1987).

Shear flow and thermal fluctuations are implemented by the Brownian
multiparticle collision dynamics approach (B-MPC) Ripoll, Winkler, and Gompper
(2007); Gompper  _et al._ (2009); Malevanets and Yeomans (2000), where
hydrodynamic interactions are neglected. In this method, stochastic collisions
between each bead and a number $\rho$ of “fluid” phantom particles of mass $m$
mimic interactions of a fluid volume surrounding a bead. The moment of a
phantom particle is Maxwellian distributed, with variance $mk_{B}T$ and mean
$(m\dot{\gamma}y\sin(2\pi t/T),0)^{T}$ at time $t$ in the presence of the
oscillating shear flow of shear rate $\dot{\gamma}$, period $T$, and
orientation along the $x$-axis. The collision process itself is implemented
via the stochastic rotation dynamics (SRD) variant of the MPC method Ihle and
Kroll (2001); Lamura  _et al._ (2001); Gompper  _et al._ (2009). Here, the
relative velocity of a bead, with respect to the center-of-mass velocity of
the bead and its related phantom particles, is rotated in the $xy-$plane by a
fixed angle $\pm\alpha$ of uniformly distributed sign. Collisions occur at
time intervals $\Delta t$, where $\Delta t>\Delta t_{p}$.

Simulations are performed for the parameters: $\alpha=130^{o}$, $\Delta
t=0.1t_{u}$, with the time unit $t_{u}=\sqrt{mr_{0}^{2}/(k_{B}T)}$, $M=\rho
m$, $\rho=5$, $\kappa_{h}r_{0}^{2}/(k_{B}T)=4\times 10^{3}$,
$\epsilon/(k_{B}T)=1$, $\sigma=r_{0}$, $\sigma_{w}=r_{0}/2$, $N=101$, and
$\Delta t_{p}=10^{-2}\Delta t$. The value of $\kappa_{h}$ ensures the length
$L=100r_{0}$ of the polymer within $1\%$ for all systems and flow conditions.

##  III Results

We study semiflexible polymers with the persistence lengths $L_{p}/L=0.5$ and
$2$. B-MPC simulations of free filaments yield the end-to-end vector
relaxation times $\tau_{r}\simeq(1.9,3.9)\times 10^{6}t_{u}$ Lamura and
Winkler (2012). The tethered polymers are initialized with beads aligned along
the $y$-direction, and are equilibrated up to $5\times 10^{6}t_{u}$.
Subsequently, data are collected up to the longest simulated time
$10^{8}t_{u}$ and averaged over half-periods with positive shear flow. The
flow is characterized in terms of the Weissenberg number
$Wi=\tau_{r}\dot{\gamma}$ and the Deborah number $De=\tau_{r}\omega$, where
$\omega=2\pi/T$ is the shear frequency. Since in oscillatory flow the strain
in a half-cycle is proportional to $\dot{\gamma}/\omega=Wi/De$, averages will
be characterized as function of this ratio in the following. The shear rate is
varied the range $0\leq Wi\lesssim 10^{3}$, while the frequency is altered
such that $0.1\leq Wi/De\leq 100$.

###  III.1 Conformational properties

We characterize the filament conformational properties by the end-to-end
vector ${\bm{R}}_{e}={\bm{r}}_{N}-{\bm{r}}_{1}$. The oscillatory shear implies
a cyclic dynamics of filament collapse, stretching, and alignment along the
flow direction (see movies movie1.mp4, movie2.mp4, movie3.mp4, and movie4.mp4
in supplementary material.).

####  III.1.1 Filament stretching

Stretching along the flow direction is characterized by the maximum extension
$R_{ex,max}$ of the polymer along the $x$-axis in a cycle, or, equivalently,
by the average deficit length-ratio $\varepsilon=1-\left\langle
R_{ex,max}\right\rangle/L$, where the average is performed over periods.

Figure 1: Mean deficit length-ratio $\varepsilon$ as a function of the strain
$Wi/De$ for $L_{p}/L=0.5$ (left), $2$ (right) and $De=10$ ($\circ$, black),
$25$ ($\bullet$, red), $50$ ($\star$, green), $100$ ($\triangle$, blue), and
$200$ ($\square$, magenta). The slope of the full lines is $-1/3$.

Figure 1 displays $\varepsilon$ as a function of the strain $Wi/De$. The data
for the various $De$ nearly collapse when plotted as function of $Wi/De$, a
feature applying to other quantities too. Three regimes can be identified,
which we will denote as low- ($Wi/De\lesssim 1$), intermediate- ($1\lesssim
Wi/De\lesssim 10$), and high-strain regime ($Wi/De\gtrsim 10$). In the low-
strain regime, either the shear rate or/and the period are so small that the
shear is hardly able to deform the polymer compared to the equilibrium value.
With increasing strain beyond unity, the deficit length decreases very rapidly
and for $Wi/De>10$ a power-law regime seems to appears, with an exponent of
$-1/3$ for $De=10$. This exponent is consistent with experimental results on
tethered semiflexible DNA molecules under steady shear Ladoux and Doyle
(2000). The high-strain regime corresponds to high values either of the shear
rate or the shear period, implying a large polymer stretching
($\varepsilon\lesssim 10^{-2}$), and the approximations of Ref. Ladoux and
Doyle (2000) in the derivation of shear-rate dependence apply. For $De>10$,
there seem to be some deviations from the power law $-1/3$, or the power-law
regime has not yet been reached; a detailed analysis is hampered by necessary
very large shear rates. Noteworthy, stiffer polymers are stronger stretched
than flexible ones.

####  III.1.2 Height above wall

Confinement to the positive semi-space breaks spatial symmetry and leads to a
geometry-induced polymer stretching normal to the wall. The dependence of the
height above the wall, defined as the average distance $\left\langle
y_{N}\right\rangle$ of bead $N$ from the wall, on strain is depicted in Fig.
2.

Figure 2: Mean polymer height above wall as a function of the strain $Wi/De$
for $L_{p}/L=0.5$ (left), $2$ (right), and $De=10$ ($\circ$, black), $25$
($\bullet$, red), $50$ ($\star$, green), $100$ ($\triangle$, blue), and $200$
($\square$, magenta). The solid lines indicate the slope $-1/3$. $\left\langle
y_{N,0}\right\rangle$ is the mean polymer height at equilibrium.

The height decreases with increasing strain according to the power-law
$\left\langle y_{N}\right\rangle\sim(Wi/De)^{-1/3}$ for high strain,
independent of the considered stiffness. The range of the scaling regime
roughly coincides with that in Fig. 1 for the deficit length. Moreover, our
simulations yield the dependence $\left\langle
y_{N}^{2}\right\rangle\sim(Wi/De)^{-2/3}\sim\left\langle
y_{N}\right\rangle^{2}$ in the high $Wi/De$ regime. This dependence agrees
with the scaling prediction for the fluctuations normal to a wall of a free-
draining semiflexible polymers under shear flow Ladoux and Doyle (2000).

####  III.1.3 Polymer deformation

More detailed insight onto the polymer deformation is gained by the
distribution function of the end-to-end vector $R_{e}=|\bm{R}_{e}|$ (see Fig.
3). In the absence of shear, the distribution function exhibits a single peak
at $R_{e}/L\simeq 1$ for $L_{p}/L=2$ corresponding to a fully elongated
configuration. With decreasing persistence length, the distribution broadens
and the peak moves to $R_{e}/L\simeq 0.87$. For the stiffer chain at small
strain, $Wi/De=2.5$, the conformations are already substantially affected by
shear, and $R_{e}$ values over a broad range,
$0.3\lesssim|\bm{R}_{e}/L|\lesssim 0.85$, are nearly equally probable. This
range broadens with decreasing $De$ at fixed $Wi$, the probabilities in the
plateau-like regime decrease, and a peak appears in the vicinity of full
stretching. The latter indicates an increasing probability of strongly-
stretched polymers. At the same time, $P(R_{e})$ decreases at smaller end-to-
end vectors. We obtain approximately the same distribution functions for the
more flexible polymer.

Figure 3: Probability distribution function of the polymer end-to-end distance
$|\bm{R}_{e}|$ for $L_{p}/L=0.5$ (left) and $2$ (right), and the Weissenberg
number $Wi=0$ ($\Box$, green) as well as $Wi=500$ for $De=10$ ($\circ$,
black), $50$ ($\bullet$, red), $200$ ($\star$, blue).

The variations of the distribution function with $De$ affect the average end-
to-end vector $R_{e}$. Associated with the overall decrease of the
distribution function by a change of $Wi/De$ from zero to $Wi/De=2.5$, is a
small drop of the average mean end-to-end distance, $\left\langle
R_{e}\right\rangle$. A further increase of $Wi/De$ implies a monotonic
swelling of $\left\langle R_{e}\right\rangle$. At the same time,
$\left\langle(R_{e}-\langle R_{e}\rangle)^{2}\right\rangle$ swells first with
increasing strain, reaches a maximum at $Wi/De\approx 8$, and decreases then
by a power-law with the exponent $-1/3$.

###  III.2 Dynamical properties

Figure 4 displays conformations of polymers during one period.

Figure 4: Polymer conformations at the times $t/T=1,1.25,1.5,1.75$ (black,
red, green, and blue lines, respectively) after equilibration for
$L_{p}/L=0.5$ ((a), (b)) and $L_{p}/L=2$ ((c), (d)), $Wi=500$, and $De=10$
((a), (c)) and $200$ ((b), (d)). (See movies movie1.mp4 for the case (a),
movie2.mp4 for the case (c), movie3.mp4 for the case (b), and movie4.mp4 for
the case (d) in supplementary material.)

At high strain, a polymer flips back and forth following the applied shear.
However, the flipping mechanism strongly depends on the interplay of shear and
bending rigidity. In the case of a more flexible polymer, the polymer recoils
after the flow changes sign, then flips and elongates (see supplementary movie
movie1.mp4). In contrast, the stiffer polymer bends, assumes a U-shaped
conformation during flipping, and finally elongates (see the supplementary
movie movie2.mp4). Smaller strain disfavors chain flip when the flow changes
sign. Here, a polymer remains coiled without flipping (see Fig. 4 (b) and the
supplementary movie movie3.mp4), and, consequently, remains essentially in the
positive half-space $x>0$. This behavior is also observed for higher bending
stiffness. In the latter case, however, the chain flip occurs more frequently
than in the flexible case with the polymer recoiling and reorienting (see Fig.
4 (d) and the supplementary movie movie4.mp4).

The described mechanism of chain flipping is supported by the probability
distribution functions $P(\Phi)$ of the inclination angle $\Phi$, which is
defined as the angle between the polymer center-of-mass position vector and
the flow direction (Fig. 5).

Figure 5: (Top) Probability distribution function of the polymer inclination
angle $\Phi$ based on data from positive half-cycle for $L_{p}/L=0.5$, $De=10$
($\circ$, black), $50$ ($\bullet$, red), $200$ ($\star$, blue), and $Wi=500$.
The distribution functions for negative half-cycle are mirror symmetric with
respect to $\Phi/\pi=1/2$. (Bottom) Average of $\cos\Phi$ as a function of the
strain $Wi/De$ for $L_{p}/L=0.5$ and $De=10$ ($\circ$, black), $25$
($\bullet$, red), $50$ ($\star$, green), $100$ ($\triangle$, blue), and $200$
($\square$, magenta). The inset shows the variance
$\langle(\cos\Phi-\langle\cos\Phi\rangle)^{2}\rangle$ as a function of $Wi/De$
for $L_{p}/L=0.5$ and $De=10$ ($\circ$, black), $25$ ($\bullet$, red), $50$
($\star$, green), $100$ ($\triangle$, blue), and $200$ ($\square$, magenta).
Figure 6: Time dependence of the polymer center-of-mass position along the
flow ($\circ$, blue), $x_{cm}$, and the gradient ($\bullet$, red), $y_{cm}$,
directions for $L_{p}/L=0.5$, $Wi=500$, and $De=10$(a), $50$(b), and $200$(c).
The black lines indicate the externally applied shear flow with arbitrary
amplitude, and the dashed lines are guides for the eye.

For large strain, the polymers spend most of their time aligned with the flow
direction (the main peak is located at $\Phi<\pi/2$) during the positive half-
cycle of the shear oscillation (Fig. 5 (left)), and opposite to the flow
direction (the main peak is at $\Phi>\pi/2$) during the negative half-cycle of
the shear oscillation (the distribution functions for the negative half-cycle
are mirror symmetric with respect to $\Phi=\pi/2$). This suggests that for
high values of the strain the inclination angle and the shear flow are of the
same sign before the shear flow changes sign. Interestingly, for the lowest
considered value $Wi/De=2.5$, the distribution function $P(\Phi)$ for the
positive and negative half-cycles are indistinguishable for the considered
persistence lengths—$P(\Phi)$ shows two broad peaks of comparable height
symmetric with respect to the shear direction ($\Phi=\pi/2$). This indicates
that polymers can either flip or keep their orientation with equal
probability, when the shear rate changes sign.

The average, $\langle\cos\Phi\rangle$ (Fig. 5 (right)), of the inclination
angle depends weakly on shear and increases logarithmically,
$\langle\cos\Phi\rangle\sim\log(Wi/De)$, on strain, while its variance has a
maximum for intermediate values of strain.

####  III.2.1 Center-of-mass motion

Figure 6 illustrates the time dependence of the center-of-mass,
$\bm{r}_{cm}=(x_{cm},y_{cm})^{T}$, during three periods for $L_{p}/L=0.5$. For
the highest considered strain $Wi/De=50$ (Fig. 6(a)), $x_{cm}$ exhibits a
nearly in-phase periodic motion with a small phase shift $\theta$, and the
dynamics is apparently no longer harmonic. The component $y_{cm}$ shows very
narrow peaks whenever $x_{cm}$ changes sign and becomes very small as soon as
the polymer is stretched. Remarkably, $y_{cm}$ exhibits frequency doubling. We
like to stress that the frequency doubling is independent of the considered
quantity used to characterize the center-of-mass motion. Specifically, we
considered the time dependence of the inclination angle $\Phi$ and that of the
distance $R_{cm}=\sqrt{x_{cm}^{2}+y_{cm}^{2}}$. While $\Phi$ exhibits the same
frequency as the applied flow, $R_{cm}$ shows frequency doubling. With
decreasing $Wi/De$, the motion is still periodic and becomes more harmonic
with a larger phase shift (Fig. 6(b)). The $y_{cm}$ peaks are now broader and
their amplitude decreases. Finally, for the smallest value $Wi/De=2.5$ (Fig.
6(c)), $x_{cm}$ partially follows the external flow with a reduced amplitude
which is no longer symmetric with respect to $x=0$ due to entropic effects,
and the magnitude of the angle $\theta$ further increases. The time sequence
of $y_{cm}$ indicates that the polymer is no longer fully elongated along the
wall when the flow reaches its extreme values. The peaks are very broad and no
clear indication of a frequency doubling can be observed, since the amplitude
of $y_{cm}$ is significantly reduced. The comparable data for $L_{p}/L=2$ show
no appreciable differences. In order to evaluate the phase shift $\theta$, the
values of $x_{cm}$ are fitted to the function $A\sin(\omega t+\theta)$, where
$A$ and $\theta$ are fitting parameters. The values of $\theta$ are in the
range $(-\pi/2,0)$ and appear to be independent on the persistence length
(Fig. 7). The extracted dependence on $Wi/De$ can well be described by the
function $\theta=\operatorname{arccot}\left(-aWi/De\right)$, with $a=0.054$.

Figure 7: Dependence of the phase shift $\theta$ on the strain $Wi/De$ for
$L_{p}/L=0.5$ ($\star$, blue) and $2$ ($\bullet$, red), with $Wi=500$. The
full line represents the fitted functional dependence
$2\operatorname{arccot}\left(-0.054Wi/De\right)/\pi$.  Figure 8: Time auto-
correlations of the center-of-mass $x_{cm}$ (left) and $y_{cm}$ (right)
Cartesian coordinates for $L_{p}/L=2$, $Wi=500$, and $De=10$ ($\circ$, black),
$50$ ($\bullet$, red), $200$ ($\star$, blue).

Further evidence of the discussed frequency doubling is obtained by the
autocorrelation functions of the Cartesian center-of-mass coordinates (Fig.
8). Consistent with Fig. 6, the correlation functions display a periodic
motion. The $x_{cm}$-component reveals the same frequency as the external
flow. More importantly, the $y_{cm}$-component clearly shows frequency
doubling for $Wi/De\gtrsim 10$. The wall at $y=0$ reflects the polymer every
time the flow moves the polymer from one side to the other. This process
occurs twice in every flow cycle causing the observed frequency doubling in
the dynamics along the $y$ axis.

###  III.3 Internal polymer dynamics

The snapshots of Fig. 4 reveal strong polymer conformational changes during
oscillation cycles. We characterize the appearing internal polymer dynamics by
the normal mode expansion

${\bm{r}}_{i}=\sum_{n=1}^{N-1}{\bm{A}}_{n}(t)\sin[q_{n}(i-1)]\;,\;\;\;i=1,\ldots,N,$
(5)

with the wave vectors $q_{n}=(n-1/2)\pi/(N-1)$ ($n=1,\ldots,N-1$). The normal
mode amplitudes $\bm{A}_{n}(t)=(A_{x,n}(t),A_{y,n}(t))^{T}$ are given by

$\displaystyle A_{x,n}$ $\displaystyle=$
$\displaystyle\frac{2}{N-1}\sum_{i=1}^{N}r_{x,i}\sin[q_{n}(i-1)],$ (6)
$\displaystyle A_{y,n}$ $\displaystyle=$
$\displaystyle\frac{2}{N-1}\sum_{i=1}^{N}r_{y,i}\sin[q_{n}(i-1)].$ (7)

The dependence of the variance of the mode amplitudes
$\langle\delta\bm{A}_{n}^{2}\rangle=\langle\left(\bm{A}_{n}-\langle\bm{A}_{n}\rangle\right)^{2}\rangle$
on the mode number is illustrated in Fig. 9 for the stiffer polymer.

Figure 9: Variance of the mode amplitudes $A_{x,n}$ (empty symbols connected
by full lines) and $A_{y,n}$ (filled symbols connected by dashed lines) as
functions of mode number $n$ for $L_{p}/L=2$ with $Wi=500$ and $De=10$
($\square$, black), $50$ ($\circ$, red), $200$ ($\star$, blue). The slopes of
full and dashed lines are $-4$ and $-2$, respectively.

Evidently, the variances of $A_{x,n}$ and $A_{y,n}$ are distinctly different
for the considered strains. In the absence of shear, both components of the
mode amplitudes exhibit the mode number dependence $n^{-4}$, in accordance
with the semiflexible character of the considered polymers Harnau, Winkler,
and Reineker (1995); Aragón and Pecora (1985). Note that the plateau for large
mode numbers is a consequence of the discreteness of the polymer with a finite
number of modes. Along the flow direction, the $n^{-4}$ dependence persists
even under shear due to the strong polymer stretching. However, we observe a
pronounced strain effects on $A_{y,n}$ for small mode numbers $3\leq n\leq 5$.
The indicated dependence $n^{-2}$, typical for flexible polymers Harnau,
Winkler, and Reineker (1995); Rouse (1953); Doi and Edwards (1986), suggests
that the polymer acquires a flexible polymer-like behavior on larger length
scales, with a crossover to semiflexibility, $n^{-4}$ dependence, on smaller
scale, along the gradient direction. A similar behavior is found for the more
flexible chain.

The dynamics of the mode amplitudes with $n=1,2$ is analyzed in terms of the
mode-autocorrelation functions $\left\langle
A_{x,n}(t)A_{x,n}(0)\right\rangle$ and $\left\langle
A_{y,n}(t)A_{y,n}(0)\right\rangle$. Results for various strains and the
persistence length $L_{p}/L=2$ are presented in Fig. 10. The simulation data
are fitted by the function

$\displaystyle f(t)=$ $\displaystyle\ A\left[\exp(-\gamma
t/T)-1\right]+B\left[\cos(\Omega\omega t)-1\right]$
$\displaystyle+C\left[\cos(2\Omega\omega t)-1\right]+1,$ (8)

where the dimensionless factors $\gamma$ and $\Omega$ characterizes the
relaxation time and possible variations of the frequency with respect to the
imposed $\omega$, respectively. The correlation functions for $De=50$ and
$De=200$ are well fitted by Eq. (III.3), while data for $De=10$ do not allow
for a suitable fitting. The factor $\gamma$ of the $x$ component decreases
with increasing mode number, $n$, revealing a faster relaxation on larger
length-scales along the flow direction. Thereby, $\gamma$ depends on $De$ and
is larger for $De=200$. It increases with increasing ratio $Wi/De$, reflecting
the strong dependence of the relaxation process on the external field.

At equilibrium, the variance $\langle\delta\bm{A}_{n}^{2}\rangle$ of the
amplitudes is determined by the polymer relaxation times Doi and Edwards
(1986). The mode-number dependence of $\gamma$ is inconsistent with that of
$\langle\delta\bm{A}_{n}^{2}\rangle$ (Fig. 9), which indicates a strong
influence of the external field on the time dependence of the internal
relaxation process. Along the $y$ direction, $\gamma$ is approximately
constant at fixed strain, indicating that the chain relaxation time does not
vary significantly on large scales along the shear direction.

The factor $\Omega$ is close to unity and oscillations occur with the external
frequency. For both, the $x$ and $y$ components, the term with $2\omega$
contributes to the oscillations of the correlation function, i.e., the
observed frequency doubling found for the center-of-mass motion is also
reflected in the internal dynamics.

##  IV Conclusions

We have analyzed the conformational and dynamical properties of semiflexible
polymers tethered at an impenetrable wall under oscillatory shear flow. We
identify three different regimes in terms of polymer deformation as a function
of strain. At small strain, $Wi/De\lesssim 1$, the polymer structures are
close to the equilibrium conformations and they are hardly affected by flow.
For intermediate strain, $1<Wi/De<10$, the polymer is gradually stretched and
its end-to-end distance increases logarithmically. Above $Wi/De\gtrsim 10$,
the deficit length exhibits a power-law decrease with increasing strain, in
agreement with experimental results under steady shear flow Doyle, Ladoux, and
Viovy (2000); Ladoux and Doyle (2000).

Figure 10: Autocorrelation function of the mode amplitudes for the modes $n=1$
along the (a) $x$\- and (b) $y$-direction as well as for $n=2$ along the (c)
$x$\- and (d) $y$-direction as a function of time with $L_{p}/L=2$, $Wi=500$,
and $De=10$ ($\circ$, black), $50$ ($\bullet$, red), $200$ ($\star$, blue).
The black dotted and solid lines are fits of Eq. ((III.3)) to the data for
$De=50$ and $De=200$, respectively.

The polymer conformations are tightly coupled with their dynamics. At high
strain, the polymers follow the external oscillations of the flow with a small
phase shift, and are strongly stretched and aligned with the wall. The polymer
dynamics is far less enslaved by the flow field for intermediate strains,
which is reflected by a larger phase difference between the flow-induced
polymer structures and the external flow. Most importantly, this allows for
more swollen polymer-like conformations and a higher polymer density above the
wall.

As might be expected for the considered geometry, conformational properties
normal to the wall exhibit frequency doubling compared to the external
frequency $\omega$, for example, the autocorrelation function of the
$y$-center-of-mass coordinate. While the properties along the $x$ direction
are typically symmetric in terms of positive and negative values, those along
$y$ are not and only positive values are possible. Hence, positive $y$ values
correspond to positive and negative $x$ values, which determines the frequency
doubling.

In absence of chirality and torsion, we expect semiflexible polymers in three
dimensions to exhibit similar features in terms of conformational changes.
This is reflected in the similarity of our results with experiments of DNA
molecules tethered to a wall under shear flow Ladoux and Doyle (2000).
Hydrodynamic interactions influence the polymer dynamics, however, the extent
depends on their stiffness. As far as semiflexible polymers are concerned,
hydrodynamic interactions are of minor importance, so that the observed
conformational and dynamical features will also be present for polymers
embedded in a fluid even in three dimensions. Hydrodynamic interactions are
expected to more severely affect the dynamics of flexible polymers. They are
again weak for rather stretched conformations Ladoux and Doyle (2000), but can
dominate the dynamics in the coiled state and accelerate the polymer motion.
In particular, this affects the internal dynamics along the $y$ direction,
with a rather Zimm- than Rouse-like mode dependence of the relaxation times on
larger length scales Doi and Edwards (1986).

Our simulations provide valuable insight into the properties of polymers
exposed to oscillatory flows. Specifically, we demonstrate that such
simulations are feasible for experimentally relevant parameters. Our studies
are a first step only toward systematic investigations of polymer properties
under oscillatory shear and we hope that they will prompt further simulations
along that line.

## Supplementary material

The supplementary movies show the animations of the polymer dynamics over two
shear periods for the cases reported in Fig. 4. Movie 1 is relative to the
case with $L_{p}/L=0.5$, $Wi=500$, $De=10$, movie 2 is relative to the case
with $L_{p}/L=2$, $Wi=500$, $De=10$, movie 3 is relative to the case with
$L_{p}/L=0.5$, $Wi=500$, $De=200$, and movie 4 is relative to the case with
$L_{p}/L=2$, $Wi=500$, $De=200$.

###### Acknowledgements.

A. L. thanks G. G., R. G. W., and co-workers for their hospitality during the
stay at the Forschungszentrum Jülich and acknowledges support from the DAAD
through the “Research Stays for University Academics and Scientists” program.

## Data Availability Statement

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

## References

  * Bird, Armstrong, and Hassager (1987) R. B. Bird, R. C. Armstrong, and O. Hassager,  _Dynamics of Polymer Liquids_ , Vol. 1 (John Wiley & Sons, New York, 1987).
  * Larson (1999) R. G. Larson,  _The Structure and Rheology of Complex Fluids_ (Oxford University Press, New York, 1999).
  * Rubinstein and Colby (2003) M. Rubinstein and R. C. Colby,  _Polymer Physics_ (Oxford University Press, Oxford, 2003).
  * Larson and Desai (2015) R. G. Larson and P. S. Desai, “Modeling the rheology of polymer melts and solutions,” Ann. Rev. Fluid Mech. 47, 47 (2015).
  * Shaw and MacKnight (2018) M. T. Shaw and W. J. MacKnight,  _Introduction to polymer viscoelasticity_ (John Wiley & Sons, 2018).
  * Schroeder (2018) C. M. Schroeder, J. Rheol. 62, 371 (2018).
  * Perkins  _et al._ (1995) T. T. Perkins, D. E. Smith, R. G. Larson, and S. Chu, Science 268, 83 (1995).
  * Winkler (2006a) R. G. Winkler, “Semiflexible polymers in shear flow,” Phys. Rev. Lett. 97, 128301 (2006a).
  * Winkler (2010) R. G. Winkler, “Conformational and rheological properties of semiflexible polymers in shear flow,” J. Chem. Phys. 133, 164905 (2010).
  * Wilhelm and Frey (1996) J. Wilhelm and E. Frey, “Radial distribution function of semiflexible polymers,” Phys. Rev. Lett. 77, 2581 (1996).
  * Götter  _et al._ (1996) R. Götter, K. Kroy, E. Frey, M. Bärmann, and E. Sackmann, Macromolecules 29, 30 (1996).
  * Harnau, Winkler, and Reineker (1996) L. Harnau, R. G. Winkler, and P. Reineker, “Dynamic structure factor of semiflexible macromolecules in dilute solution,” J. Chem. Phys. 104, 6355 (1996).
  * Everaers  _et al._ (1999) R. Everaers, F. Jülicher, A. Ajdari, and A. C. Maggs, “Dynamic fluctuations of semiflexible filaments,” Phys. Rev. Lett. 82, 3717 (1999).
  * Samuel and Sinha (2002) J. Samuel and S. Sinha, “Elasticity of semiflexible polymers,” Phys. Rev. E 66, 050801 (2002).
  * Le Goff  _et al._ (2002) L. Le Goff, O. Hallatschek, E. Frey, and F. Amblard, “Tracer studies on f-actin fluctuations,” Phys. Rev. Lett. 89, 258101 (2002).
  * Winkler (2003) R. G. Winkler, “Deformation of semiflexible chains,” J. Chem. Phys. 118, 2919 (2003).
  * Petrov  _et al._ (2006) E. P. Petrov, T. Ohrt, R. G. Winkler, and P. Schwille, “Diffusion and segmental dynamics of double-stranded DNA,” Phys. Rev. Lett. 97, 258101 (2006).
  * Bird  _et al._ (1987) R. B. Bird, C. F. Curtiss, R. C. Armstrong, and O. Hassager,  _Dynamics of Polymer Liquids_ , Vol. 2 (John Wiley & Sons, New York, 1987).
  * Öttinger (1996) H. C. Öttinger,  _Stochastic Processes in Polymeric Fluids_ (Springer, Berlin, 1996).
  * Winkler (2006b) R. G. Winkler, “Semiflexible polymers in shear flow,” Phys. Rev. Lett. 97, 128301 (2006b).
  * Munk  _et al._ (2006) T. Munk, O. Hallatschek, C. H. Wiggins, and E. Frey, “Dynamics of semiflexible polymers in a flow field,” Phys. Rev. E 74, 041911 (2006).
  * Hur and Shaqfeh (2000) J. S. Hur and E. S. G. Shaqfeh, “Brownian dynamics simulations of single DNA molecules in shear flow,” J. Rheol. 44, 713 (2000).
  * Jendrejack, de Pablo, and Graham (2002) R. M. Jendrejack, J. J. de Pablo, and M. D. Graham, “Stochastic simulations of dna in flow: Dynamics and the effects of hydrodynamic interactions,” J. Chem. Phys. 116, 7752 (2002).
  * Hsieh and Larson (2004) C.-C. Hsieh and R. G. Larson, “Modelling hydrodynamic interaction in Brownian dynamics: Simulation of extensional and shear flows of dilute solutions of high molecular weight polystyrene,” J. Rheol. 48, 995 (2004).
  * Liu, Ashok, and Muthukumar (2004) S. Liu, B. Ashok, and M. Muthukumar, “Brownian dynamics simulations of beadrod-chain in simple shear flow and elongational flow,” Polymer 45, 1383 (2004).
  * Celani, Puliafito, and Turitsyn (2005) A. Celani, A. Puliafito, and K. Turitsyn, “Polymers in linear shear flow: A numerical study,” Europhys. Lett. 70, 464 (2005).
  * Ryder and Yeomans (2006) J. F. Ryder and J. M. Yeomans, “Shear thinning in dilute polymer solutions,” J. Chem. Phys. 125, 194906 (2006).
  * Ripoll, Winkler, and Gompper (2006) M. Ripoll, R. G. Winkler, and G. Gompper, “Star polymers in shear flow,” Phys. Rev. Lett. 96, 188302 (2006).
  * Sendner and Netz (2008) C. Sendner and R. R. Netz, EPL 81, 54006 (2008).
  * He  _et al._ (2009) G.-L. He, R. Messina, H. Löwen, A. Kiriy, V. Bocharova, and M. Stamm, “Shear-induced stretching of adsorbed polymer chains,” Soft Matter 5, 3014 (2009).
  * Zhang  _et al._ (2009) Y. Zhang, A. Donev, T. Weisgraber, B. J. Alder, M. G. Graham, and J. J. de Pablo, “Tethered dna dynamics in shear flow,” J. Chem. Phys. 130, 234902 (2009).
  * Huang  _et al._ (2010) C.-C. Huang, R. G. Winkler, G. Sutmann, and G. Gompper, “Semidilute polymer solutions at equilibrium and under shear flow,” Macromolecules 43, 10107 (2010).
  * Huang  _et al._ (2011) C.-C. Huang, G. Sutmann, G. Gompper, and R. G. Winkler, “Tumbling of polymers in semidilute solution under shear flow,” EPL 93, 54004 (2011).
  * Huang, Gompper, and Winkler (2012) C.-C. Huang, G. Gompper, and R. G. Winkler, “Non-equilibrium relaxation and tumbling times of polymers in semidilute solution,” J. Phys.: Condens. Matter 24, 284131 (2012).
  * Lang, Obermayer, and Frey (2014) P. S. Lang, B. Obermayer, and E. Frey, Phys. Rev. E 89, 022606 (2014).
  * Nikoubashman and Howard (2017) A. Nikoubashman and M. P. Howard, Macromolecules 50, 8279 (2017).
  * Kong  _et al._ (2019) X. Kong, Y. Han, W. Chen, F. Cui, and Y. Li, Soft Matter 15, 6353 (2019).
  * Romo-Uribe (2021) A. Romo-Uribe, J. Appl. Polym. Sci. 138, 49712 (2021).
  * Shee  _et al._ (2021) A. Shee, N. Gupta, A. Chaudhuri, and D. Chaudhuri, Soft Matter 17, 2120 (2021).
  * Nikoubashman (2021) A. Nikoubashman, J. Chem. Phys. 154, 090901 (2021).
  * Shaqfeh (2005) E. S. G. Shaqfeh, J. Non-Newtonian Fluid Mech. 130, 1 (2005).
  * Smith, Finzi, and Bustamante (1992) S. B. Smith, L. Finzi, and C. Bustamante, Science 258, 1122 (1992).
  * Perkins  _et al._ (1994) T. T. Perkins, S. R. Quake, D. E. Smith, and S. Chu, Science 264, 822 (1994).
  * Doyle, Ladoux, and Viovy (2000) P. S. Doyle, B. Ladoux, and J.-L. Viovy, “Dynamics of a tethered polymer in shear flow,” Phys. Rev. Lett. 84, 4769 (2000).
  * Lueth and Shaqfeh (2009) C. A. Lueth and E. S. G. Shaqfeh, Macromolecules 42, 9170 (2009).
  * Delgado-Buscalioni (2006) R. Delgado-Buscalioni, “Cyclic motion of a grafted polymer under shear flow,” Phys. Rev. Lett. 96, 088303 (2006).
  * Maier, Seifert, and Rädler (2002) B. Maier, U. Seifert, and J. O. Rädler, Europhys. Lett. 60, 622 (2002).
  * Cherstvy and Petrov (2014) A. G. Cherstvy and E. P. Petrov, “Modeling dna condensation on freestanding cationic lipid membranes,” Phys. Chem. Chem. Phys. 16, 2020 (2014).
  * Herold, Schwille, and Petrov (2010) C. Herold, P. Schwille, and E. P. Petrov, “Dna condensation at freestanding cationic lipid bilayers,” Phys. Rev. Lett. 104, 148102 (2010).
  * Chattopadhyay and Marenduzzo (2007) A. K. Chattopadhyay and D. Marenduzzo, Phys. Rev. Lett. 98, 088101 (2007).
  * Lamura and Winkler (2012) A. Lamura and R. G. Winkler, “Semiflexible polymers under external fields confined to two dimensions,” J. Chem. Phys. 137, 244909 (2012).
  * Hyun  _et al._ (2011) K. Hyun, M. Wilhelm, C. O. Klein, K. S. Cho, J. G. Nam, K. H. Ahn, S. J. Lee, R. H. Ewoldt, and G. H. McKinley, “A review of nonlinear oscillatory shear tests: Analysis and application of large amplitude oscillatory shear (LAOS),” Prog. Poly. Sci. 36, 1697 (2011).
  * Rogers (2018) S. Rogers, “Large amplitude oscillatory shear: Simple to describe, hard to interpret,” Physics Today 71, 34 (2018).
  * Chen  _et al._ (2005) Y.-L. Chen, M. D. Graham, J. J. de Pablo, K. Jo, and D. C. Schwartz, Macromolecules 38, 6680 (2005).
  * Thomas, DePuit, and Khomami (2009) D. G. Thomas, R. J. DePuit, and B. Khomami, J. Rheol. 53, 275 (2009).
  * Lamura and Winkler (2019) A. Lamura and R. G. Winkler, “Tethered semiflexible polymer under large amplitude oscillatory shear,” Polymers 11, 737 (2019).
  * Zhou and Schroeder (2016) Y. Zhou and C. M. Schroeder, Phys. Rev. Fluids 1, 053301 (2016).
  * Ripoll, Winkler, and Gompper (2007) M. Ripoll, R. G. Winkler, and G. Gompper, “Hydrodynamic screening of star polymers in shear flow,” Eur. Phys. J. E 23, 349 (2007).
  * Gompper  _et al._ (2009) G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler, “Multi-particle collision dynamics: A particle-based mesoscale simulation approach to the hydrodynamics of complex fluids,” Adv. Polym. Sci. 221, 1 (2009).
  * Kapral (2008) R. Kapral, “Multiparticle collision dynamics: Simulations of complex systems on mesoscale,” Adv. Chem. Phys. 140, 89 (2008).
  * Swope  _et al._ (1982) W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, “A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters,” J. Chem. Phys. 76, 637 (1982).
  * Allen and Tildesley (1987) M. P. Allen and D. J. Tildesley,  _Computer Simulation of Liquids_ (Clarendon Press, Oxford, 1987).
  * Malevanets and Yeomans (2000) A. Malevanets and J. M. Yeomans, “Dynamics of short polymer chains in solution,” Europhys. Lett. 52, 231–237 (2000).
  * Ihle and Kroll (2001) T. Ihle and D. M. Kroll, “Stochastic rotation dynamics: A Galilean-invariant mesoscopic model for fluid flow,” Phys. Rev. E 63, 020201(R) (2001).
  * Lamura  _et al._ (2001) A. Lamura, G. Gompper, T. Ihle, and D. M. Kroll, “Multiparticle collision dynamics: Flow around a circular and a square cylinder,” Europhys. Lett. 56, 319–325 (2001).
  * Ladoux and Doyle (2000) B. Ladoux and P. S. Doyle, “Stretching tethered dna chains in shear flow,” Europhys. Lett. 52, 511 (2000).
  * Harnau, Winkler, and Reineker (1995) L. Harnau, R. G. Winkler, and P. Reineker, J. Chem. Phys. 102, 7750 (1995).
  * Aragón and Pecora (1985) S. R. Aragón and R. Pecora, Macromolecules 18, 1868 (1985).
  * Rouse (1953) P. E. Rouse, J. Chem. Phys. 21, 1272 (1953).
  * Doi and Edwards (1986) M. Doi and S. F. Edwards,  _The Theory of Polymer Dynamics_ (Clarendon Press, Oxford, 1986).