# Curiosity Based Reinforcement Learning on Robot Manufacturing Cell

Mohammed Sharafath Abdul Hameed, Md Muzahid Khan, Andreas Schwung
Department of Automation Technology
South Westphalia University of Applied Sciences
Soest, Germany.
_sharafath.mohammed, khan.mdmuzahid<EMAIL_ADDRESS>
###### Abstract

This paper introduces a novel combination of scheduling control on a flexible
robot manufacturing cell with curiosity based reinforcement learning.
Reinforcement learning has proved to be highly successful in solving tasks
like robotics and scheduling. But this requires hand tuning of rewards in
problem domains like robotics and scheduling even where the solution is not
obvious. To this end, we apply a curiosity based reinforcement learning, using
intrinsic motivation as a form of reward, on a flexible robot manufacturing
cell to alleviate this problem. Further, the learning agents are embedded into
the transportation robots to enable a generalized learning solution that can
be applied to a variety of environments. In the first approach, the curiosity
based reinforcement learning is applied to a simple structured robot
manufacturing cell. And in the second approach, the same algorithm is applied
to a graph structured robot manufacturing cell. Results from the experiments
show that the agents are able to solve both the environments with the ability
to transfer the curiosity module directly from one environment to another. We
conclude that curiosity based learning on scheduling tasks provide a viable
alternative to the reward shaped reinforcement learning traditionally used.

###### Index Terms:

reinforcement learning, manufacturing cell, curiosity based learning, planning
robots

##  I Introduction

Reinforcement Learning (RL) is becoming a popular algorithm to solve complex
games and tasks like Atari [1], physics simulation [2], Go [3] and Chess [4],
Robotics [5], and optimization problems [6] etc. In these approaches, the RL
agent is given a continuous supply of dense extrinsic reward by the
environment for the action it takes. While providing a reward or defining a
reward function is not a problem in common gaming tasks, it quickly becomes
cumbersome in complex engineering tasks like the optimization problem [7]. And
traditionally when the rewards are not directly available, they need to be
shaped to guide the agent in the direction of an optimal solution. This
approach creates three problems: one, this introduces a huge inductive bias in
the reward scheme that it cripples the generalized learning of the RL agent
from the start. So much so that it might be only able to solve in the known
solution space; two, reward shaping is a notoriously difficult engineering
problem, where un-optimized rewards could break the learning process or
provide only sub-optimal solutions; and finally, it creates solutions which
cannot be transferred to new environments, since the rewards have been heavily
tuned for a particular environment. In this paper, we propose to apply
curiosity based learning to a flexible Robot Manufacturing Cell.

In [8], a Robot Manufacturing Cell (RMC) was proposed as an optimization
problem similar to the scheduling problems encountered in the production
planning. And in [9], the same RMC was solved using central and distributed
advantage actor critic (A2C), through hand engineered reward signals. Similar
reward shaping methods are also followed in [10] [11] [12], where rewards for
each of the scheduling problem is hand engineered to suit the environment they
work in. This obviously creates a problem of non-generalized solutions, where
the learning of an agent from one environment is not transferable to an
another, even though they both solve similar optimization problems.

To overcome this limitation of the reinforcement learning application in
engineering problems, in this paper we apply curiosity [13] to the RMC
environment. Curiosity uses the intrinsic motivation of the agent to create
rewards based on the surprise factor of the agents familiarity with the
environment states. Through the application of curiosity in reinforcement
learning of the RMC we overcome all the problems stated above. The main
contributions of this paper as follows:

  * •

We propose curiosity as an effective tool to eliminate the inductive bias
induced by the reward shaping normally used in reinforcement learning of
scheduling tasks. Further, curiosity eliminates complex reward shaping for
engineering tasks like scheduling.

  * •

We show that a direct transfer of the reward function from one manufacturing
environment to another is possible without any tuning requirements.

  * •

We combine curiosity with gradient monitoring, an effective approach to
stabilize training. We apply the curiosity-based RL to two multi-robot
environments where we particularly show the improvements achieved by the
proposed approach in both environments.

The paper is organized as follows. Section II introduces the learning problem,
RMC, in detail, describing the operations and setup. Section III introduces
the RL algorithms used in this paper along with further optimizations that
were used, while Section IV, Section IV-A, and Section IV-B explain the
curiosity module along with other optimizations used. Section V describes the
results that were achieved on the RMC through the application of the RL along
with its optimizations, while Section VI provides the conclusion with future
scope of work.

##  II Robot Manufacturing Cell

In this paper, RL model is applied to a cooperative self-learning RMC [8]. The
RMC consists of two industrial robots, an Adept Cobra i600 SCARA-robot and an
ABB IRB1400 6-DOF-robot, both of which share a common work platform as shown
in Figure 1. The setup has six stations totally, two input stations (IB1,
IB2), three processing or work stations (M1, M2, M3) and one output station
(OB). There are two different types of work-piece in this system (WP-1, WP-2).
The WP-1 goes through IB1, M1, M2, M3, and then OB, while WP-2 goes through
IB2, M2, M1, M3 and then OB. Both the robots have their own proprietary
control software. Hence a supervisory Siemens Programmable Logic Control (PLC)
system is developed to control and coordinate the robot movement. Two
experiments are conducted to demonstrate the applicability of the curiosity
based approach on scheduling problems. The first experiment is conducted on a
simple structured Robot Manufacturing Cell (sRMC) environment, similar to [9]
and the second experiment is conducted on a graph-network based Robot
Manufacturing Cell (gRMC) environment, which was developed in [14]. gRMC is a
scalable variant of the sRMC environment, where the machines and the buffers
are represented by the nodes of a graph while the connections between them are
expressed as edges. A work-piece is transported between the nodes when the
edge is activated.

Figure 1: Schematic drawing of the robot manufacturing cell

In sRMC no buffers are considered for the processing stations environment,
i.e., there is no queue in front or after the processing stations. The RL
agent’s actions here are the movements of the work-pieces between the
processing stations. While in gRMC three buffer-stations are considered, one
after each of the processing stations (MB1, MB2, MB3) as shown in Figure 2.
Here the agent’s action is the activation of the edges connecting the nodes.
If a node is activated then the work-piece is moved from the ’from’-node to
the ’to’-node. Two different types of edges are defined in the system, to move
the different work-pieces. Both the work-pieces are embedded with the RFID
chip which gives the required information to the robots.

Figure 2: In the gRMC environment, the nodes are indexed from 0-8 and are
colored according to their function, while the different colored lines
represent the two work-pieces to be processed and transported.

##  III Reinforcement Learning

RL is a field of machine learning where an agent is trained to take a sequence
of decisions on an environment, modeled as a Markov Decision Process (MDP), to
maximize the cumulative reward signal. The RL system has two major components
namely, the agent and the environment. The agent interacts with the
environment, takes an action $a_{t}$ based on the agent’s policy
$\pi$($a_{t}|s_{t})$ from its current state $s_{t}$ and moves to a new state
$s_{t+1}$. Based on the action taken by the agent, the environment sends a
reward signal $r_{t}$ to the agent. The goal of the agent is to find an
optimal policy that maximizes the cumulative reward [15]. Hence, the agent
tries to maximize the expected return $\mathbb{E}[R_{t}|\pi]$ with the
discounted reward given by $R$ = $\sum_{t}\gamma^{t}r_{t}$ and discount factor
0 $\leq\gamma\leq$ 1\. The MDP is defined by the tuple
($\mathcal{S,A,P,R},p_{0}$), with the set of states $\mathcal{S}$, the finite
set of actions $\mathcal{A}$, a transition model $\mathcal{P}$ for a given
state $s$ and an action $a$, the reward function $\mathcal{R:S\times A\times
S\rightarrow\mathbb{R}}$, which provides a reward $r$ for each state
transition $s_{t}\rightarrow s_{t+1}$ with an initialization probability
$p_{0}$.

The policy is defined using either value based methods or policy gradient
methods. Value based methods use either a state-value function,
$v(s)=\mathbb{E}[R_{t}|S_{t}=s]$ or a state-action-value function,
$q(s,a)=\mathbb{E}[R_{t}|S_{t}=s,A_{t}=a]$. The policy is then defined by an
$\epsilon-greedy$ strategy where greedy actions are given by
$\pi(a|s)=argmax(q(s,\mathcal{A}))$, where the agent exploits, or a non-greedy
action is given, where the agent explores. In policy gradient methods, a
parameterized policy, $\pi(a|s,\theta)$, is used to take action. The policy
can still be improved using value functions as seen in [1]. The the agent’s
objective is to derive an optimal policy, $\pi^{*}(a|s),$ that maximizes the
reward collection. Correspondingly, the policy is updated in a way to maximize
a cost function $J(\theta_{t})$ given by,

$\theta_{t+1}=\theta_{t}+\rho\nabla J(\theta_{t}),$ (1)

where $\theta$ is the policy parameters of $\pi$ and $\rho$ is the learning
rate.

###  III-A Advantage Actor Critic (A2C):

Advantage Actor Critic is a policy gradient algorithm which combines policy
gradient method with value based approach like SARSA [15]. In A2C, the actor
tries to find the optimal policy and the critic evaluates the actor’s policy.
Here, the actor refers to the policy $\pi(a|s,\theta_{1})$ and the critic
refers to the value function $\nu(s,\theta_{2})$, where $\theta_{1}$ is the
parameter of the actor and the $\theta_{2}$ is the parameter of the critic.
The cost function equation of the A2C is given by Eqn. 2 and 3 respectively.

$J(\theta)=E_{\sim\pi_{\theta}}\bigg{[}\sum\limits_{(s_{t},a_{t})\epsilon}log\pi_{\theta}(a_{t},s_{t})A_{\pi_{\theta}}\bigg{]}$
(2)

$A_{\pi_{\theta}}(s_{t},a_{t})=Q_{\pi_{\theta}}(s_{t},a_{t})-V_{\pi_{\theta}}(s_{t})\\\
$ (3)

##  IV Curiosity Driven Learning

Children often learn through intrinsic motivation such as curiosity.
Psychologists define a behavior as intrinsically motivated when a human works
for his own sake. Similarly, intrinsic motivation in RL keeps an agent engaged
in the exploration of new knowledge and skills when the rewards are rare [16].
There are lots of approaches to develop the intrinsic reward of an agent like
prediction error, prediction uncertainty, or development of the forward model
trained along with the policy of an agent [17]. These approaches drive the
agent to explore the environment in the absence of a dense extrinsic reward.
In this paper, we use prediction error as the intrinsic reward, which is
produced by the agent itself with a model called Intrinsic Curiosity Model
(ICM). ICM is composed of two subsystems: an inverse dynamic model that
produce an intrinsic reward signal and forward dynamic model that outputs a
sequence of actions which guides the agent to maximize the reward signal [13].
Figure 3 represents the schematic architecture of the ICM. The ICM encodes the
current state $s_{t}$ and the next state $s_{t+1}$ into the feature of current
state $\phi(s_{t})$ and feature of next state $\phi(s_{t}+1)$. The feature of
current state $\phi(s_{t})$ and the feature of next state $\phi(s_{t}+1)$ are
trained with inverse dynamic model that predicts action $\hat{a}_{t}$ [13].
Mathematically, this transition can be represented as,

$\hat{a}_{t}=g(s_{t},s_{t+1};\theta_{I})$ (4)

where $g$ is a soft-max distribution function of all possible actions that
could be taken and $\theta_{I}$ is the neural network parameter. The loss
function of the inverse dynamic model can represent by the following,

${\min_{\theta_{I}}}=L_{I}(\hat{a_{t}},a_{t})$ (5)

where $L_{I}$ the loss function that calculates the deviation between the
actual action $a_{t}$ and the predicted action $\hat{a_{t}}$. The forward
dynamic model takes $\phi({s_{t}})$ and action $a_{t}$ as the inputs and
outputs the predicted feature representation $\hat{\phi}(a_{t+1})$ of next
state $s_{t+1}$. The forward dynamic model can be represented by the
following,

$\hat{\phi}(s_{t+1})=f(\phi(s_{t}),a_{t};\theta_{F})$ (6)

where $f$ is the learned function and $\theta_{F}$ is the ICM neural network
parameter. The loss function of the forward model is given below:

$L_{F}(\phi(s_{t}),\hat{\phi}(s_{t+1})=\frac{1}{2}||\phi(s_{t+1})-\phi(s_{t})||^{2}_{2}$
(7)

where $L_{F}$ is the loss function that calculates the variance between the
predicted feature representation of next state $\hat{\phi}(a_{t+1})$ and the
actual predicted feature of next state $\phi(s_{t}+1)$. So, the intrinsic
reward $r_{t}^{i}$ can be calculated as:

$r_{t}^{i}=\frac{\eta}{2}||\phi(s_{t+1})-\phi(s_{t})||^{2}_{2}$ (8)

where $\eta>0$ is the scaling factor of the intrinsic reward signal. The
intrinsic reward is generated by the agent at time $t$. Besides intrinsic
rewards, the agent may also get extrinsic reward $r_{t}^{e}$ from the
environment.

Figure 3: Graphical representation of the Intrinsic Curiosity Module

To generate the intrinsic reward signal and to optimize the ICM, the inverse
dynamic model and forward dynamic model loss should be combined. So, the
overall optimization problem of ICM that helps the agent to learn can be
represented as:

$\min_{\theta_{P},\theta_{I},\theta_{F}}\Bigg{[}{-\lambda
E_{\pi(s_{t};\theta_{P})}}\bigg{[}\sum\nolimits_{t}r_{t}\bigg{]}+(1-\beta)L_{I}+\beta
L_{F}\Bigg{]}$ (9)

where 0$\leq$$\beta$$\leq$1 is a scaling factor that weights the inverse model
against the forward model, and ${\lambda>0}$ is a scaling factor of policy
gradient loss against the intrinsic reward [13]. In this paper, to speed up
the learning process of the RL agents two optimization techniques called
Gradient Monitoring and Curriculum Learning are used. They are explained
further below.

###  IV-A Gradient Monitoring

RL has a non stationary learning problem. This makes ensuring stability in
learning a critical part of RL. In this paper an optimization method called
Gradient Monitoring (GM) is used to stabilize the RL training of the RMC
environments. GM provides for a faster convergence and better generalization
performance of the trained RL agent, for a small computational overhead, by
steering the gradients during the learning process towards the most important
weights. The added benefit of using GM is that it allows for a continuous
adjustment of the neural network’s utilization capacity, i.e., the size of the
neural network need not be tuned with as much care since GM automatically
derives the required size of the neural network required for the training
process [18].

###  IV-B Curriculum Learning

Humans learn things efficiently, when the tasks are organized in such a way
that the complexity levels increase gradually. This strategy is called
Curriculum Learning (CL). CL is a training process where the criterion of the
training process is co-related with different set of weights. It is a form of
re-weighting of the training distribution. At the beginning, the associated
weights represent the simpler task, that can be learned most easily. The next
training criterion involves a small change in the weighting of task that
increases the probability of sampling slightly more complex task. At the end,
the re-weighting of the tasks is uniform, and the final target is trained
[19]. In RL, sometimes it becomes difficult for the RL agent to handle complex
tasks directly and to solve this problem CL technique is used.

##  V Experimental Results

In this paper, two experiments are conducted. First sRMC is solved using ICM
and CL, since the final task of planning of 20 WP-1 and 20 WP-2 was too
difficult for the agent. Then in the same experiment, GM is introduced which
helps the agent to converge without the use of CL. Next in the second
experiment, the ICM used in sRMC is used in the gRMC environment. Here, only
CL is used since the overhead GM did not provide a corresponding gain. The
results thence obtained are discussed below in two parts based on their
environment, sRMC and gRMC. The parts output in each of the case is reported,
along with the rewards collected by in setup.

###  V-A sRMC:

In this section, first the results of curriculum learning without gradient
monitoring approach in sRMC environment is discussed. The result in parts
output, WP-1 and WP-2, with respect to iterations is shown in Figure 4. Here,
the target is 20 work-pieces for each agent. Initially, without the additional
optimizations the agents were able to solve only seven work-pieces each as
opposed to the target of 20 each. Hence, CL was introduced to divide the final
target into smaller targets of increasing difficulty, i.e. target of 5, 10,
15, and 20 work-pieces of WP-1 and WP-2. Both the agents are trained for
$1e^{5}$ episodes and the convergence for 20 target is achieved after
$0.78e^{5}$ episodes.

Figure 4: Parts Output (sRMC without gradient monitoring)

To increase the convergence speed of the agents, GM was introduced into the
agent learning the sRMC environment. The improvement posed by GM was good
enough to remove the CL completely from the training regime of sRMC
environment. This quickened the learning progress as shown in Figure 4. A
target of 20 work-pieces for each agent was set directly rather than dividing
the targets into smaller tasks. The agents solved the target directly in
$5e^{3}$ steps compared to the $0.78e^{5}$ without GM. With this, the
objective of training the agents on the sRMC environment was completed. The
results graph shows an interesting trend. Figure 6 shows the combined,
intrinsic and extrinsic, rewards collected by the agents. The rewards graph
can be used to describe the agent’s learning behavior. Normally, in extrinsic
reward-based RL algorithm, the slope of the reward curves becomes zero after
reaching a solution. But the intrinsic reward-based system works differently.
When the agent remains curious about completing the target, it generates the
intrinsic reward by itself. But once the agent finishes exploring its
environment and finds the an optimal policy, the agent remains no more curious
and loses its motivation, hence receives less reward for the same task. That
means, intrinsic rewards obviously reduce once the agent has converged. This
is illustrated in the Figure 6, where after the iteration number $1e5$, the
rewards collected reduces. The next objective is to transfer the ICM
incorporated here in to another environment without any architectural changes.

Figure 5: Parts Output (sRMC with gradient monitoring) Figure 6: Combined
Rewards (sRMC with gradient monitoring)

###  V-B Graph network based curriculum learning:

In the next experiment, the intrinsic curiosity module is plugged to train the
agents on gRMC. The learning process was a lot more difficult in this setup.
This difficulty was also noted in the reward shaping method, where the agent
did not converge in regular training. Hence, first CL was introduced in to the
training regime. The agent’s target in curriculum learning was designed
keeping in mind the difficulty i.e., targets of 1, 2, 3, 4, etc., were used.
The agent solved the targets until a target of seven work-pieces each as shown
in Figure 7, although the training times were extremely high. Another
optimization trial using GM produced similar results, but instead increased
the training times marginally. This could be because of the complex nature of
the gRMC interaction and the higher size of the action space. When looking at
Figure 8 though, the agent starts exploring the environment from its early
stage and the agent reaches to the peak of exploration faster compared to
sRMC. When the agent learns to solve its target, the curiosity of the agent
reduces, and this leads to a reduction of further exploration of the
environment. In addition to that, the combined rewards collection by the agent
in gRMC is much more higher compared to the sRMC. The objective of solving the
a different environment like gRMC using the same architecture of ICM as used
in sRMC is achieved.

Figure 7: Parts Output (gRMC without gradient monitoring) Figure 8: Combined
Rewards (gRMC without gradient monitoring)

##  VI Conclusion

RL on industrial scheduling problems are difficult to train because of the
inherent bias related to hand tuned rewards. A novel combination of the
flexible RMC and a curiosity based RL has been proposed in this study. We
found that the application of the curiosity based RL provides good solutions
even with very scarce external reward. Additionally, we found that the reward
function architecture used in one environment could be easily transferred to
another, which was is impossible in hand tuned rewards. In future work, we
plan to improve the performance of the agent by using a hybrid reward function
providing rewards with little inductive bias for a quicker initial learning
but with enough generalized performance. This combined with a gRMC with a
custom message passing network could provide highly generalized solutions that
can be directly transferable from one robot to another.

## References

  * [1] Volodymyr Mnih et al.  “Human-level control through deep reinforcement learning”  In _nature_ 518.7540 Nature Publishing Group, 2015, pp. 529–533
  * [2] Emanuel Todorov, Tom Erez and Yuval Tassa  “MuJoCo: A physics engine for model-based control”  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 5026–5033
  * [3] David Silver et al.  “A general reinforcement learning algorithm that masters chess, shogi, and Go through self-play”  In _Science_ 362, 2018, pp. 1140–1144
  * [4] Julian Schrittwieser et al.  “Mastering atari, go, chess and shogi by planning with a learned model”  In _arXiv preprint arXiv:1911.08265_ , 2019
  * [5] S. Gu, E. Holly, T. Lillicrap and S. Levine  “Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates”  In _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , 2017, pp. 3389–3396
  * [6] Bruno Cunha, Ana M. Madureira, Benjamim Fonseca and Duarte Coelho  “Deep Reinforcement Learning as a Job Shop Scheduling Solver: A Literature Review” 923, Advances in intelligent systems and computing, 2194-5357  Springer, 2020, pp. 350–359
  * [7] Kaishu Xia et al.  “A digital twin to train deep reinforcement learning agent for smart manufacturing plants: Environment, interfaces and intelligence”  In _Journal of Manufacturing Systems_ , 2020
  * [8] D. Schwung, F. Csaplar, A. Schwung and S.. Ding  “An application of reinforcement learning algorithms to industrial multi-robot stations for cooperative handling operation”  In _2017 IEEE 15th International Conference on Industrial Informatics (INDIN)_ , 2017, pp. 194–199
  * [9] Andreas Schwung, Dorothea Schwung and Mohammed Sharafath Abdul Hameed  “Cooperative Robot Control in Flexible Manufacturing Cells: Centralized vs. Distributed Approaches”  In _2019 IEEE 17th International Conference on Industrial Informatics (INDIN)_ , 2019, pp. 233–238
  * [10] In-Beom Park, Jaeseok Huh, Joongkyun Kim and Jonghun Park  “A Reinforcement Learning Approach to Robust Scheduling of Semiconductor Manufacturing Facilities: IEEE Transactions on Automation Science and Engineering, 1-12”  In _IEEE Transactions on Automation Science and Engineering_ , 2020, pp. 1–12
  * [11] Thomas Gabel and Martin Riedmiller  “Adaptive reactive job-shop scheduling with reinforcement learning agents”  In _International Journal of Information Technology and Intelligent Computing_ 24.4, 2008
  * [12] Bernd Waschneck et al.  “Optimization of global production scheduling with deep reinforcement learning”  In _Procedia CIRP_ 72.1, 2018, pp. 1264–1269
  * [13] Deepak Pathak, Pulkit Agrawal, Alexei A. Efros and Trevor Darrell  “Curiosity-Driven Exploration by Self-Supervised Prediction”  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) Workshops_ , 2017
  * [14] Mohammed Sharafath Abdul Hameed and Andreas Schwung  “Reinforcement Learning on Job Shop Scheduling Problems Using Graph Networks”, 2020  arXiv:2009.03836
  * [15] Richard S Sutton and Andrew G Barto  “Reinforcement learning: An introduction”  MIT press, 2018
  * [16] Nuttapong Chentanez, Andrew G Barto and Satinder P Singh  “Intrinsically motivated reinforcement learning”  In _Advances in neural information processing systems_ , 2005, pp. 1281–1288
  * [17] Yuri Burda et al.  “Large-Scale Study of Curiosity-Driven Learning”  In _International Conference on Learning Representations_ , 2018
  * [18] Mohammed Sharafath Abdul Hameed, Gavneet Singh Chadha, Andreas Schwung and Steven X. Ding  “Gradient Monitored Reinforcement Learning”, 2020  arXiv:2005.12108
  * [19] Yoshua Bengio, Jérôme Louradour, Ronan Collobert and Jason Weston  “Curriculum learning”  In _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 41–48


# Dark discrete breather modes in monoaxial chiral helimagnet with easy-plane
anisotropy

I.G. Bostrem1, E.G. Ekomasov2,3, J. Kishine4,5, A.S. Ovchinnikov1,6, and V.E.
Sinitsyn1 1 Institute of Natural Science and Mathematics, Ural Federal
University, Ekaterinburg 620002, Russia  2 Bashkir State University, Institute
of Physics and Technology, Ufa 450076, Russia  3 Tyumen State University,
Institute of Physics and Technology, Tyumen 625003, Russia  4 Division of
Natural and Environmental Sciences, The Open University of Japan, Chiba
261-8586, Japan  5 Institute for Molecular Science, 38 Nishigo-Naka, Myodaiji,
Okazaki, 444-8585, Japan  6 Institute of Metal Physics, Ural Division, Russian
Academy of Sciences, Ekaterinburg 620219, Russia

###### Abstract

Nonlinearity and discreteness are two pivotal factors for an emergence of
discrete breather excitations in various media. We argue that these
requirements are met in the forced ferromagnetic phase of the monoaxial chiral
helimagnet CrNb3S6 due to specific domain structure of the compound. The
stationary, time-periodic breather modes appear as the discrete breather
lattice solutions whose period mismatches with a system size. Thanks to easy-
plane single-ion anisotropy intrinsic to CrNb3S6, these modes are of the dark
type with frequencies lying within the linear spin-wave band, close to its
bottom edge. They represent cnoidal states of magnetization, similar to the
well-known soliton lattice ground state, with differing but limited number of
embedded $2\pi$-kinks. Linear stability of these dark breather modes is
verified by means of the Floquet analysis. Their energy controlled by two
parameters, namely the breather lattice period and amplitude, falls off
linearly with a growth of the kink number. These results may pave a new path
to design spintronic resonators on the base of chiral helimagnets.

###### pacs:

Valid PACS appear here

##  I Introduction

Nonlinear magnetic phenomena, where soliton-like waves are involved, have been
extensively discussed for many decades, and a wealth of information has been
accumulated on the subject since that time Kosevich1990 . Among various
nonlinear coherent structures, such as kinks, vortices, or monopoles,
breathers may appear as generic solutions of nonlinear dynamics in a wide
range of physical systems Infeld1990 . Unlike other soliton-like solutions,
which preserve their shape as they propagate, the breathers are characterized
by spatial localization and periodic oscillations in time. This localization
is realized in perfectly regular systems as opposed to disorder-induced
Anderson localization.

When considering the discrete aspects of nonlinear lattice problems, the
modified concept of discrete breathers (DB) is introduced Flach1998 ;
Flach2008 ; Dmitriev2016 . Discreteness is essential to preventing resonances
between the breather modes and the system characteristic frequencies. Contrary
to continuous solitons, DBs do not require integrability for their existence
and stability. They are not confined to a specific spatial dimension and have
an amplitude dependent frequency. The DB basic properties must be complemented
by their division into the stationary and moving modes Bickham1997 . These
nonlinear excitations have been observed in a variety of physical systems,
where discreteness arises either as an intrinsic property of the underlying
lattice, such as lattice dynamics of halide-bridged transition metal complexes
Swanson1999  and PbSe crystal Manley2019 , layered crystal insulator
Russel2007 , or in artificially designed arrays of coupled Josephson junctions
Binder2000 ; Trias2000 , coupled pendula Cuevas2009 , micromechanical
cantilevers Sato2003 , optical waveguides Eisenberg1998 , in compressed
diatomic granular one-dimensional (1D) crystal Boechler2010 , and in
electrical circuit with series-connected tunnel diodes Narahara2020 .

As for DB excitations in spin lattices, they have been predicted for
antiferromagnetic chains Lai1999  with subsequent direct experimental
verification Schwarz1999 ; Sato2004 . Currently, considerable numerical and
theoretical progress has been achieved in understanding the nature of DB in
Heisenberg spin chains with additional interactions, for instance, onsite
anisotropy or next-nearest-neighbor exchange interactionsWallis1995 ;
Rakhmanova1996 ; Zolotaryuk2001 ; Khalack2003 ; Lakshmann2014 ; Lakshmanan2018
; Kamburova2019 . Recently, formation of discrete breathers has been analyzed
for weak ferromagnetic chains where the presence of the Dzyaloshinskii-Moriya
(DM) interaction leads to a small canting between the interacting moments
Djoufack2016 ; Kavitha2016 . The DB solutions have been examined in dynamics
of the 1D array of magnetic particles (dots) with the easy-plane anisotropy
and interparticle magnetic dipole interaction Pylypchuk2015 .

Most of realistic discrete systems may be reformulated for appropriate
continuous media that often provide an adequate description of nonlinear
properties and, in some cases, analytical expressions in closed form may be
derived Kosevich1990 ; Kosevich1974 ; Ivanov1979 ; Barjaktar1994 ;
Papanicolaou1997 . However, it is quite possible that some nonlinear phenomena
in continuous models can not appear in their discrete counterparts, and more
comprehensive treatment is required to gain a deeper insight into the problem
Gorbach2005 ; Quan2009 . One such issue is whether the discrete system
supports breather lattice (BL) solutions, which represent a regularly ordered
array of single breathers. These periodic solutions may be expressed in terms
of Jacobi elliptic functions Byrd  and occur naturally for the continuous
nonlinear equations, such as the sine-Gordon equation McLachlan1994 ;
Kevrekidis2001 ; Fu2007a , sinh-Gordon equation Fu2007c , Korteweg–de Vries
(KdV) equation Kevrekidis2003 , modified KdV equation Kevrekidis2004 ; Yan2008
; Zhao2009 , nonlinear Schrödinger equation Fu2007b . Applications of the
periodic breather modes to various continuous media have been discussed in
Refs.Tankeyev2009 ; Tankeyev2010 ; Smagin2009 ; Yin2018 . The BL solution has
been predicted for continuous $\beta$-Fermi-Pasta-Ulam (FPU) chainQuan2008 ,
finding periodic DB configurations in a triangular $\beta$-FPU lattice has
been recently reported Babicheva2021 . As for discrete spin lattices,
spatially periodic breather type solutions have been traced numerically in
pioneering work Rakhmanova1998 , where they were named intrinsic localized
spin modes of dark types by analogy with the optic solitonsWeiner1988 ;
Bao2018 . However, to our knowledge, their analytical expressions are still
lacking.

In this paper, we argue that the phase of forced ferromagnetism in the
monoaxial chiral helimagnet CrNb3S6 assumes the existence of BL excitations.
These modes are intrinsically stationary due to specific domain structure of
CrNb3S6 samples composed of 1 $\mu$m-wide grains with different, left- or
right-hand, crystallographic structural chirality Togawa2015 . We demonstrate
analytically and numerically that these spatially periodic standing solutions
with amplitude varying in time may be categorized into four types, depending
on the position their frequencies with respect to the spectrum of linear spin
waves. This position is set by a value and sign of single-ion anisotropy, and,
as a consequence, the only type of BL solutions, which is allowed by fairly
small easy-plane anisotropy in CrNb3S6, is the so-called dark breather modes
emerging within spin-wave band, close to its bottom edge. They represent
cnoidal states of magnetization, similar to the soliton lattice ground state,
with different number of embedded kinks. The salient features of the ”bottom
dark” breather lattice are its period incommensurate with a system size and
its energy falling off linearly with a growth of a number of kinks. The linear
stability of these modes is verified by means of Floquet analysis.

It should be emphasized at the outset that previous theoretical investigations
of DB modes in ferromagnetic chains with on-site anisotropies Wallis1995 ;
Rakhmanova1996 ; Rakhmanova1998  were inspired by experimental studies of spin
dynamics in the linear chain compound CsFeCl3 Schmid1994  for which intrinsic
single-ion anisotropy exceeds intersite exchange coupling. A similar situation
was dealt in Ref. Pylypchuk2015  where magnetic dots with a strong easy-plane
anisotropy are coupled by weak magnetic dipole interaction. In contrast, our
treatment is targeted at search of DB solutions for the case when exchange
coupling between the nearest moments is much larger than on-site magnetic
anisotropy.

The paper is organized as follows. In Sec. II, we describe a model and give
classification of periodic breather solutions. In Sec. III, we present results
of numerical simulation of breather lattices based on discrete equations of
spin motion. In Sec. IV, we present a detailed analysis of dark breather modes
emerging near the bottom edge of the spin-wave spectrum. The conclusions and
discussions are given in Sec. V.

##  II Model

The model Hamiltonian of the chiral monoaxial helimagnet is of the form

$\mathcal{H}=-2J\sum_{n}\textbf{S}_{n}\cdot\textbf{S}_{n+1}+A\sum_{n}\left(S^{z}_{n}\right)^{2}$
$+D\sum_{n}\left[\textbf{S}_{n}\times\textbf{S}_{n+1}\right]_{z}-H\sum_{n}S^{z}_{n},$
(1)

where the first term is the exchange coupling along the chiral axis ($z$-axis)
with $J>0$. The second describes the single-ion anisotropy with the constant
$A$, while the third does Dzyaloshinskii-Moryia interaction of the strength
$D$. The last term denotes Zeeman coupling with an external magnetic field $H$
directed along the $z$ axis. The magnetic field is assumed to be strong
enough, so that in the ground state all spins are ordered along the external
field direction. This forced ferromagnetic arrangement requires
$H>2S\left(\sqrt{4J^{2}+D^{2}}-J+A\right)$ Kishine2015 .

The appropriate spin variables $s^{\pm}_{n}=\left(S^{x}_{n}\pm
iS^{y}_{n}\right)/S$ and $s^{z}_{n}=S^{z}_{n}/S$ describe spin deviations from
the ground state, where $S$ is the magnitude of spin.

The equations of motion for these variables become

$\frac{i\hbar}{2JS}\frac{d}{dt}s^{+}_{n}=s^{+}_{n}\left(s^{z}_{n+1}+s^{z}_{n-1}\right)-s^{z}_{n}\left(s^{+}_{n-1}+s^{+}_{n+1}\right)$
$-2Bs^{+}_{n}s^{z}_{n}+i\frac{D}{2J}s^{z}_{n}\left(s^{+}_{n-1}-s^{+}_{n+1}\right)+\frac{H}{2JS}s^{+}_{n},$
(2)

where $s^{z}_{n}=\sqrt{1-\left|s^{+}_{n}\right|^{2}}$ and $B=A/2J$.

Time-dependent solutions $s^{+}_{n}=s_{n}(t)\exp\left(ikna-i\omega t\right)$
with the wave number $k$ and the frequency $\omega$ are to be found, where the
amplitude $s_{n}$ is called the transverse spin accumulation and $a$ being the
lattice constant. Substituting this form into Eq. (2) and collecting
separately real and imaginary parts we obtain the system

$\Omega
s_{n}=s_{n}\left(\sqrt{1-s^{2}_{n+1}}+\sqrt{1-s^{2}_{n-1}}\right)-2Bs_{n}\sqrt{1-s^{2}_{n}}$
$-\sqrt{1-s^{2}_{n}}\left(s_{n-1}+s_{n+1}\right)\sqrt{1+\frac{D^{2}}{4J^{2}}}\cos(ka+\delta),$
(3)
$\frac{ds_{n}}{d\tau}=\sqrt{1-s^{2}_{n}}\left(s_{n-1}-s_{n+1}\right)\sqrt{1+\frac{D^{2}}{4J^{2}}}\sin(ka+\delta),$
(4)

where $\tau=t/t_{0}$ is the dimensionless time with $t_{0}=\hbar/(2JS)$, the
effective frequency $\Omega=\left(\hbar\omega-H\right)/2JS$, and the phase
shift is determined from $\delta=\textrm{tan}^{-1}\left(D/2J\right)$.

Breather spin modes are expected to exist in perfect discrete magnetic chains
due to salient nonlinearity in the exchange and anisotropy interactions. These
modes may lie inside the linear wave band and splitting off from its upper or
lower edges, so-called DB solutions of the dark type Rakhmanova1998 . By
contrast, the bright type modes can appear at frequencies either just above
the upper band edge or just below the lower edge of the linear spin wave
spectrum.

In the small spin deviation limit, the linear dispersion curve for spin wave
is obtained from Eq.(3) as

$\Omega(k)=2-2B-2\sqrt{1+\frac{D^{2}}{4J^{2}}}\cos ka.$ (5)

The top and bottom zone boundaries of the plane wave spectrum take the values
$\Omega(\pi)=2-2B+2\sqrt{1+\frac{D^{2}}{4J^{2}}}$, and
$\Omega(0)=2-2B-2\sqrt{1+\frac{D^{2}}{4J^{2}}}$, respectively.

Concerning the DB wave number, it coincides with that of the simple or conical
spirals in the static case when $ds_{n}/d\tau=0$, i.e.
$ka=-\textrm{tan}^{-1}\left(D/2J\right)$. It immediately follows from Eq. (3).

To start with, we focus on the upper band edge and consider both bright and
dark excitations. As orientational order of transversal spin components in
these breather modes is close to that of the corresponding linear spin waves,
i.e. to antiferromagnetic one, it is convenient to introduce the envelope
function $\psi(z)=(-1)^{n}s_{n}$, where $z=na$, which varies slowly in space.

As will be seen in the next Section, the DB lattice modes may be compared with
periodic solutions of the nonlinear Duffing equationRakhmanova1996 ;
Tankeyev2010

$\frac{d^{2}\psi}{d\tilde{z}^{2}}-\alpha\psi+\beta\psi^{3}=0,$ (6)

which can be deduced from Eq. (3) if to keep only the leading cubic
nonlinearity (for details, see Appendix A). Here, $\tilde{z}=z/a$ is the
dimensionless coordinate expressed in lattice units, and

$\alpha=\frac{\Omega-\Omega(\pi)}{\sqrt{1+\frac{D^{2}}{4J^{2}}}},\quad\beta=-\frac{\Omega(\pi)}{2\sqrt{1+\frac{D^{2}}{4J^{2}}}}.$
(7)

For the case $\alpha>0$, $\beta>0$ the bright BL solution centered at
$\tilde{z}_{0}$ is obtained as

$\psi^{(+)}_{b}(\tilde{z})=\left[\sqrt{\frac{2c}{\beta}+\frac{\alpha^{2}}{\beta^{2}}}+\frac{\alpha}{\beta}\right]^{\frac{1}{2}}\textrm{cn}\left[\frac{4K}{L^{(+)}_{b}}\left(\tilde{z}-\tilde{z}_{0}\right),\kappa^{2}\right],$
(8)

where $\textrm{cn}(\ldots)$ is Jacobi elliptic function with the modulus

$\kappa^{2}=\frac{1}{2}\frac{\alpha+\sqrt{\alpha^{2}+2\beta
c}}{\sqrt{\alpha^{2}+2\beta c}}.$ (9)

The constant $c$ can take only positive values for these solutions that the
superscript ”+” refers to. Together with the requirement of smallness of
$\psi^{(+)}_{b}$, it results in the condition

$0<c<\frac{\beta}{2}-\alpha.$ (10)

The constant $c$ specifies the amplitude and the bright BL period
$L^{(+)}_{b}=4K\left(\alpha^{2}+2\beta c\right)^{-\frac{1}{4}}$, where $K$ is
the elliptic integral of the first kind. For a finite size system $c$ may be
derived from the requirement that dynamics of the edge and interior spins must
be consistent (see Sec. III for details).

In the opposite case $c<0$, the solution is given by

$\psi^{(-)}_{b}(\tilde{z})=\left[\frac{\alpha}{\beta}-\sqrt{\frac{2c}{\beta}+\frac{\alpha^{2}}{\beta^{2}}}\right]^{\frac{1}{2}}\textrm{nd}\left[\frac{2K}{L^{(-)}_{b}}\left(\tilde{z}-\tilde{z}_{0}\right),\kappa^{2}\right],$
(11)

where $\textrm{nd}(\ldots)$ is the Jacobi elliptic function
$\textrm{dn}^{-1}(\ldots)$ with the modulus

$\kappa^{2}=\frac{2\sqrt{\frac{\alpha^{2}}{\beta^{2}}+\frac{2c}{\beta}}}{\frac{\alpha}{\beta}+\sqrt{\frac{\alpha^{2}}{\beta^{2}}+\frac{2c}{\beta}}},$
(12)

and the period $L^{(-)}_{b}=\sqrt{8}K\left[\alpha+\sqrt{\alpha^{2}+2\beta
c}\right]^{-\frac{1}{2}}$. In contrast to (10), the $c$-values are bounded
from below $c>-\alpha^{2}/(2\beta^{2})$.

As expected, in the limit $\kappa^{2}=1$, when $c=0$, the bright breather
lattice (8) reduces to the single breather solution

$\psi_{b}(\tilde{z})=\sqrt{\frac{2\alpha}{\beta}}\frac{1}{\cosh\left[\sqrt{\alpha}\left(\tilde{z}-\tilde{z}_{0}\right)\right]}$
(13)

with a peak shape of width $1/\sqrt{\alpha}$ centered around $\tilde{z}_{0}$,
which can be loosely be thought as a ”bound state” of a kink and antikink. In
contrast, the solution (11) approaches zero when $\kappa^{2}$ goes to 1, but
for all other $\kappa$ values the classification based on the number of kink-
antikink bound pairs may be retained.

Figure 1:  The top bright breather mode with one kink-antikink bound pair in
the envelope function built from numerical simulations (upper panel).
Transversal spin components (relative to the chain axis) are shown with a
factor 50 expansion. Lower panel: the spin-wave spectrum at $B=4.0$ (black
solid), the breather frequency $\Omega=-3.95$ (red solid) is limited from
above by
$\Omega_{u}=\frac{3}{2}\left[1-B+\sqrt{1+\left(D/2J\right)^{2}}\right]$ (blue
dashed).

The numerical solution corresponding to the bright BL mode with the envelope
function (8) is presented in Fig. 1 (for details, see Sec. III). It is known
that the existence of such intrinsic localized modes (ILSM) for a
ferromagnetic chain with nearest-neighbor interactions requires that the
strength of the single-ion anisotropy exceeds a certain critical value so that
the resulting ILSM frequencies can appear above the linear spin wave band
Lai1999 . In our case, the requirement $\beta>0$ imposes the restriction on
the allowed anisotropy strength, $B>1+\sqrt{1+D^{2}/4J^{2}}$, whereby it must
exceed the Heisenberg exchange coupling. The condition is not suitable for the
chiral helimagnet CrNb3S6, where $B\sim 0.15$ Miyadai1983  is close to the
value for the antisymmetric exchange $D/2J\sim 0.16$. Shinozaki2016
Notwithstanding, there are non-chiral magnetic materials, such as CsFeCl3 or
some quasi-1D metal-organic compounds, where this somewhat exotic situation
takes place. Anders1989

Figure 2:  The top dark breather mode with 4 successive kink/ antikinks in the
envelope function built from numerical simulations (upper panel). Transversal
spin components (relative to the chain axis) are shown with a factor 20
expansion. Lower panel: Dispersion of the linear spin waves at the easy-plane
anisotropy $B=1.0$ (solid black). The breather frequency $\Omega=1.8$ is
limited by the top edge of the spin-wave spectrum and
$\Omega_{d}=1-B+\sqrt{1+\left(D/2J\right)^{2}}$ (blue dashed).

The opposite situation, $\alpha<0$ and $\beta<0$, leads to the lattice
breather solution of the dark type

$\psi^{(+)}_{d}(\tilde{z})=\pm\sqrt{\frac{\alpha}{\beta}-\sqrt{\frac{\alpha^{2}}{\beta^{2}}+\frac{2c}{\beta}}}\,\textrm{sn}\left[\frac{4K}{L_{d}}(\tilde{z}-\tilde{z}_{0}),\kappa^{2}\right]$
(14)

with the elliptic modulus

$\kappa^{2}=\frac{\alpha+\sqrt{\alpha^{2}+2\beta
c}}{\alpha-\sqrt{\alpha^{2}+2\beta c}}$ (15)

and the period $L_{d}=4K\left(-\alpha/2+\sqrt{\alpha^{2}/4+\beta
c/2}\right)^{-\frac{1}{2}}$. Condition on the parameter $c$, which is
positive, repeats that of the bright modes (10). Note that at $c<0$ only
unbounded solutions appear what makes them inappropriate as finite-amplitude
excitations.

The corresponding single breather solution is retrieved when $\kappa^{2}\to
1$, or $c\to-\alpha^{2}/(2\beta)$,

$\psi_{d}(\tilde{z})=\pm\sqrt{\frac{\alpha}{\beta}}\tanh\left[\sqrt{-\frac{\alpha}{2}}\left(\tilde{z}-\tilde{z}_{0}\right)\right].$
(16)

It describes a kink whose inverse width is $\sqrt{|\alpha|/2}$. Based on the
form (16), breather lattices of the dark type may be categorized according to
the number of consecutively embedded S-shaped kinks/antikinks.

An example of the dark type DB excitations found numerically is shown in Fig.
2. Its envelope function fits well with the expression (14) (see Sec. III for
details). Here we just note that the derivation of the dark breather solution
is based on the assumption $\alpha<0$ thereby meaning that the effective
frequency of the localized mode cannot be upper than the top of the spin-wave
spectrum, $\Omega<\Omega(\pi)$ [see Eq. (7)]. Similarly to the case of bright
modes, this fact results from the restriction on the single-ion anisotropy
constant, $B<1+\sqrt{1+D^{2}/4J^{2}}$ originated formally from $\beta<0$.
Although this estimate allows the top dark BL solutions (14) for fairly small
easy-plane anisotropy valid for CrNb3S6, the requirement $B>0$ turns out to be
incompatible with smooth spin arrangement at boundaries (see Sec. III).

We now turn to the lower edge of the spin wave band, where bright and dark DB
modes may be also found. Similar to the above classification, the bright
localized excitations should lie in the gap below the bottom of the spin-wave
spectrum, whereas their dark counterparts may occur at frequencies above the
bottom spectrum edge.

Figure 3:  The bottom bright breather mode with one kink-antikink bound pair
in the envelope function built from numerical simulations (upper panel).
Transversal spin components (relative to the chain axis) are shown with a
factor 50 expansion. Lower panel: Dispersion curve $\Omega(k)$ for the linear
spin wave band (black solid) at the easy-axis anisotropy $B=-1.0$ (solid
black). Frequency of the breather excitations $\Omega=1.95$ (red solid) with
the lowest boundary
$\Omega_{d}=\frac{3}{2}\left[1-B-\sqrt{1+\left(D/2J\right)^{2}}\right]$ (blue
dashed).

By invoking the continuum approximation the same Eq. (6) is recovered for the
slowly varying envelope $\psi_{n}=s_{n}$, but with different coefficients

$\alpha=\frac{\Omega(0)-\Omega}{\sqrt{1+\frac{D^{2}}{4J^{2}}}},\quad\beta=\frac{\Omega(0)}{2\sqrt{1+\frac{D^{2}}{4J^{2}}}}.$
(17)

Repeating the steps that led from Eq.(6) to DB modes, the previous solutions
(8,11) and (14) are recovered as breather lattices of the bright and dark
types, respectively. However, there is a significant discrepancy in conditions
for $\Omega$ and $B$ owing to the different definition for $\alpha$ and
$\beta$. As expected, one gets $\Omega<\Omega(0)$ for the bright mode and
$\Omega>\Omega(0)$ for the dark mode what originates formally from the
requirements $\alpha>0$ or $\alpha<0$, respectively. From a magnetic
viewpoint, the emergence of these excitations is related to the restriction
imposed on the constant of the single-ion anisotropy $B$. The bright modes
($\beta>0$) appear at easy-axis anisotropy, $B<1-\sqrt{1+D^{2}/4J^{2}}$, that
rules out automatically excitations of this type for observation in CrNb3S6.
Spatial distribution of magnetic moments of this type of DB mode, found
numerically and corresponding to the model solution (11), is illustrated in
Fig. 3. In addition to this, one may find a family of solutions (8) for the
same set of parameters. In contrast, the dark modes ($\beta<0$) comply with
the opposite restriction $B>1-\sqrt{1+D^{2}/4J^{2}}$, thereby admitting both
fairly small easy-axis anisotropy ($B<0$) and that of the easy-plane type
($B>0$). Visualization of the spin arrangement for $B=0.15$ is given in Fig.
4. It is precisely this type of DB excitations, which can be approximated by
the bottom dark solution (14), is relevant to the chiral helimagnet CrNb3S6.

##  III Numerical simulations of DB modes

Figure 4:  The bottom dark breather mode with one kink in the envelope
function built from numerical simulations (upper panel). Transversal spin
components (relative to the chain axis) are shown with a factor 20 expansion.
Lower panel: The spin-wave spectrum at the easy-plane anisotropy $B=0.15$
(solid black). The breather frequency $\Omega=-0.28$ (red solid) is upper
bounded by $\Omega_{u}=1-B-\sqrt{1+\left(D/2J\right)^{2}}$ (blue dashed).

When spin deviations of the breather modes grow to be large enough, the basic
equation (6) of the theoretical treatment becomes invalid. This warrants
additional numerical simulations to verify analytical results.

Our numerical procedure follows a scheme that has been suggested in Ref.
Rakhmanova1996 . Giving an initial guess for $s_{1}$, one can let $s_{2}$ be
determined by the truncated version of Eq. (3)

$\Omega
s_{1}=-2Bs_{1}\sqrt{1-s^{2}_{1}}+s_{1}\sqrt{1-s^{2}_{2}}-s_{2}\sqrt{1-s^{2}_{1}}\sqrt{1+q^{2}_{0}}$

with $q^{2}_{0}=\frac{D^{2}}{4J^{2}}$.

After the initialization steps the general recursion takes the form

$s_{n+1}=\frac{-A_{n}\sqrt{\left(1-s^{2}_{n}\right)\left(1+q^{2}_{0}\right)}+\sqrt{1+\left(1-s^{2}_{n}\right)q^{2}_{0}-A^{2}_{n}}}{1+\left(1-s^{2}_{n}\right)q^{2}_{0}},$
(18)

where

$\displaystyle A_{n}$ $\displaystyle=$ $\displaystyle
s_{n-1}\sqrt{1-s^{2}_{n}}\sqrt{1+q^{2}_{0}}$
$\displaystyle+s_{n}\left(\Omega\right.\left.-\sqrt{1-s^{2}_{n-1}+2B\sqrt{1-s^{2}_{n}}}\right).$

The node index $n$ runs through the lattice, until it eventually reaches the
right end. The truncated equation for the right edge spins

$\displaystyle\Omega s_{L}$ $\displaystyle=$ $\displaystyle-
s_{L-1}\sqrt{1-s^{2}_{L}}\sqrt{1+q^{2}_{0}}$
$\displaystyle+s_{L}\sqrt{1-s^{2}_{L-1}}-2Bs_{L}\sqrt{1-s^{2}_{L}}.$ (19)

may be used as an inexpensive way to check convergence of the iterative
procedure. Solutions are generated by scanning $2\cdot 10^{11}$ trial runs for
$s_{1}$ from the range $[s_{1},s_{1}+h]$ with the step $h=10^{-6}$.

Examples of the four DB types, which are determined by their frequency
positions with respect to the spin wave band, are illustrated in Figs. 1-4
(their temporal oscillations may be seen in Supplemental Materials). The chain
length $L=101$ and the constant DM strength $D/2J=0.16$ have been taken in
each case, only values for $\Omega$, $B$ and $s_{1}$ varied. In each of the
plots, an additional line, bounding possible frequency range of the breather
excitations, is shown either inside or outside the spin-wave band. The
boundary is established from the requirement that the amplitude of the
appropriate single breather solution is limited to be less than 1. It can be
easily seen from Eqs.(13,16) that this is equivalent to the conditions
$2\alpha<\beta$ or $\beta<\alpha$ for bright or dark modes, respectively.

As such, only the bright breather modes may be associated with self-localized
spin-wave excitations (or, intrinsic localized modes), where a scale of the
localization may be on a scale comparable to the lattice spacing. In the case
of the bottom bright breathers, the reason for this is that both the magnetic
anisotropy and the magnetic field tend to align the moments all in the
direction of the chain axis and only the DM interaction tries to prevent them.
For the top bright breathers, the single-ion anisotropy acts already in
alliance with the antisymmetric exchange coupling, but the interior
antiferromagnetic order of moments restricts the spatial extent of these
excitations.

The latter plays a key role in a difference between the dark breather
solutions located at opposite ends of the spin-wave band. Although fairly
small easy-plane anisotropy is permissible for both of them there is an
additional restriction on $B$ originated from the requirement that the
frequency $\Omega$ of the interior spins must be the same as that of the end
moments.

Figure 5:  The envelope function of the top bright breather lattice for
$L=101$ as predicted by Eq.(8). The solid red (dotted blue) line corresponds
to analytical (numerical) calculation for $D/2J=0.16$ with $B=4.0$,
$\Omega=-3.95$. The edge spin $s_{-L/2}=0.00048823$.

Indeed, the top dark breathers may arise at

$\Omega>\Omega_{d}=1-B+\sqrt{1+q^{2}_{0}}.$ (20)

On the other hand, Eq.(19) results in

$\Omega\approx\sqrt{1-s^{2}_{L}}\left(1-2B+\sqrt{1+q^{2}_{0}}\right)$

provided the antiparallel alignment of the edge spins, $s_{L}\approx-s_{L-1}$.
Then, as an immediate consequence,

$\sqrt{1-s^{2}_{L}}>\frac{1-B+\sqrt{1+q^{2}_{0}}}{1-2B+\sqrt{1+q^{2}_{0}}}$
(21)

which leads to the obvious contradiction $\sqrt{1-s^{2}_{L}}>1$ for the easy-
plane anisotropy $B>0$.

Frequency of the bottom dark breathers is bounded above by
$\Omega_{u}=1-B-\sqrt{1+q^{2}_{0}}$. Given ferromagnetic arrangement of
moments at the edge, $s_{L}\approx s_{L-1}$, we have from (19)

$\Omega\approx\sqrt{1-s^{2}_{L}}\left(1-2B-\sqrt{1+q^{2}_{0}}\right).$

Figure 6:  The envelope function of the bottom bright breather lattice for
$L=100$ as predicted by Eq.(8). The solid red (dotted blue) line corresponds
to analytical (numerical) calculation for $D/2J=0.16$ with $B=-1.0$,
$\Omega=1.95$. The edge spin $s_{-L/2}=0.24708308$.

Therefore,

$\sqrt{1-s^{2}_{L}}<\frac{B+\sqrt{1+q^{2}_{0}}-1}{2B+\sqrt{1+q^{2}_{0}}-1},$

which amounts to the desired result $\sqrt{1-s^{2}_{L}}<1$ for any positive
$B$.

In the foregoing analysis it was tacitly assumed that the spatial arrangement
of the spin variables $s_{n}$ is smooth in the vicinity of the edges and takes
zero values only inside the chain. When, however, $s_{L}$ is small enough, the
relationship $|s_{L}|\approx|s_{L-1}|$ is not longer applicable and dark
breather solutions become possible only under a suitable choice of $B$,
ensuring matching of a spin configuration scale to the chain length. That is
the situation shown in Fig. 2. Moreover, solutions obtained thereby vary
fairly rapidly with distance that makes the continuum approximation
ineffective for them.

Figure 7:  The envelope function of the bottom bright breather lattice for
$L=101$ as predicted by Eq. (11). The solid red (dotted blue) line corresponds
to analytical (numerical) calculation for $D/2J=0.16$ with $B=-1.0$,
$\Omega=1.95$. The edge spin $s_{-L/2}=0.016876628$.

We shall now turn to fitting the analytical expressions of Sec. II to the
envelope function $\psi(z)$ obtained numerically to confirm the validity of
the continuum approximation.

Details of such a procedure for the top bright breather lattices have been
discussed by us in Ref. Bostrem2021 . The main purpose is to find the $c$
parameter and then deduce the elliptic modulus $\kappa$ from Eq. (9) with
defined $\alpha$ and $\beta$. Note that in the case of the soliton lattice
ground state, the constant $c$ is excluded from minimization of energy per
unit length, since $c$ is related with the elliptic modulus $\kappa$.
Dzyaloshinskii1964

By definition, the constant $c$ arises as the first integral of Eq. (6)

$c=\left(\frac{d\psi}{d\tilde{z}}\right)^{2}-\alpha\psi^{2}+\frac{\beta}{2}\psi^{4}.$
(22)

On the other hand, the derivative of the envelope function at the chain edge
is given by

$\displaystyle\left.\left(\frac{d\psi}{d\tilde{z}}\right)\right|_{L/2}\approx$
$\displaystyle\frac{1}{\sqrt{1+q^{2}_{0}}}\left(1+\sqrt{1+q^{2}_{0}}-\Omega-2B\right)\psi_{L/2}$
$\displaystyle+\frac{1}{2\sqrt{1+q^{2}_{0}}}\left(2B-1-\sqrt{1+q^{2}_{0}}\right)\psi^{3}_{L/2},$

which directly follows from Eq. (3), if to account the open boundary
conditions and neglect nonlinear terms including derivatives of the function
$\psi$.

It can be envisaged that the envelope function is small at edges of the chain
due to intrinsic localization of the top bright breathers. Then, Eqs. (22) and
(III) result in the transcendental equation

$c\approx\left(\alpha^{2}+\alpha\frac{2+\sqrt{1+q^{2}_{0}}}{\sqrt{1+q^{2}_{0}}}+2\frac{1+\sqrt{1+q^{2}_{0}}}{\sqrt{1+q^{2}_{0}}}\right)\psi^{2}_{L/2}.$

Comparison of the envelope functions computed both numerically and
analytically from Eq. (8) is presented in Fig. 5. Note that the top bright
breather modes may be both on-site centered Takeno1988  and intersite centered
Page1990 . This type of solutions may be found below the spin-wave band as
well (see Fig. 6), however, it requires the chain consisting of an even number
of sites, since these modes are being intersite centered only. A comparing
procedure between the numerical data and the model solution is similar to that
of the dark BL modes as explained below.

Fig. 7 shows the numerical data superimposed on the analytical curve for the
bottom bright breather, as given by Eq. (11). To render the fit with known
$\alpha$ and $\beta$, it is convenient to make use the property of the
solution, $\psi_{\textrm{min}}/\psi_{\textrm{max}}=\textrm{dn}K=1-\kappa^{2}$,
that yields $\kappa^{2}\approx 0.994298$ in this case. Then, the argument
coefficient in the Jacobi nd-function is easily recovered as

$\frac{1}{\sqrt{2}}\left(\alpha+\sqrt{\alpha^{2}+2\beta
c}\right)^{\frac{1}{2}}=\sqrt{\frac{\alpha}{2-\kappa^{2}}}\approx 0.155293.$

Figure 8:  The envelope function of the top dark breather lattice for $L=101$
as predicted by Eq. (14). The solid red (dotted blue) line corresponds to
analytical (numerical) calculation for $D/2J=0.16$ with $B=1.0$, $\Omega=1.8$.
The edge spin $s_{-L/2}=-0.079802639$. Figure 9:  The envelope function of the
bottom dark breather lattice for $L=101$ as predicted by Eq. (14). The solid
red (dotted blue) line corresponds to analytical (numerical) calculation for
$D/2J=0.16$ with $B=0.15$, $\Omega=-0.28$. The edge spin
$s_{-L/2}=0.4858544687$.

It should be noted that although the top bright BL solutions (11) are admitted
by the continuum theory, they are not reproduced numerically. This is probably
related to antiparallel alignment of magnetic moments in these excitations.

Fitting for the dark BL solutions plotted in Figs. 8,9 follows the same
template. At first, one finds the constant

$c=\psi^{2}_{\textrm{max}}\left(|\alpha|-\frac{|\beta|}{2}\psi^{2}_{\textrm{max}}\right),$

where $\psi_{\textrm{max}}$ is maximum value of the envelope function adopted
from numerical data. Next, the elliptic modulus (15) and the argument in the
Jacobi sn-function (14) can be specified.

It is evident that predictions of the treatment based on the continuum
approximation are in good agreement with their numerical counterparts. From
now on, solutions requiring either easy-axis or strong easy-plane anisotropy
are beyond the scope of our study, and we narrow focus on the bottom dark
breathers only as the most likely candidate to be detected in CrNb3S6.

Figure 10:  The envelope function of the bottom dark breather lattice
containing $N$ kinks: $N=1$ (purple) with $s_{-L/2}=0.495373$, $N=2$ (teal)
with $s_{-L/2}=0.494843$, $N=3$ (brown) with $s_{-L/2}=0.485854$. The other
parameters as in Fig. 9.

##  IV Bottom dark breathers

To gain a deeper insight into the properties of the dark bottom breathers we
start off with their classification. Numerical simulations of these
excitations for specific choices of the edge spin are depicted in Fig. 10. The
main feature of the breather lattices is that their spatial period does not
match the system size what resembles, in some respects, standing waves with
soft pinning at boundaries Pincus1960 ; Kishine2019 . Obviously, these
solutions differ in a number $N$ of embedded kinks/antikinks, or,
equivalently, by the number of nodes, that may be adopted as a criterion to
categorize them. Finite spin deviations $s_{-L/2}$ needed to generate spin
configurations of a given $N$ are slightly different from each other, but
extremely high precision typical for top dark breather solutions is not
required.

Next, we use the Floquet theory to examine linear stability of the identified
breather modes Marin1998 ; Archilla2003 ; Khalack2003 . For the purpose, one
has to study the evolution of a perturbation $\varepsilon_{n}(t)$ added to the
DB periodic solution $s^{(0)}_{n}(t)$, i.e.
$s_{n}(t)=s^{(0)}_{n}(t)+\varepsilon_{n}(t)$, for the lattice of $L$ sites
with $n=1,2,\ldots,L$. Assuming that the size of perturbation is suitably
small one may linearize the resulting equations for $\varepsilon_{n}(t)$
deduced directly from the equations of motion (2). Their explicit form looks a
rather cumbersome, but more importantly, the linearized system defines the
monodromy matrix $\hat{M}$, which maps $s_{n}(t)$ into $s_{n}(t+T)$,

$\left(\begin{array}[]{c}\varepsilon^{{}^{\prime}}_{n}(T)\\\
\varepsilon^{{}^{\prime\prime}}_{n}(T)\end{array}\right)=\hat{\mathcal{M}}\left(\begin{array}[]{c}\varepsilon^{{}^{\prime}}_{n}(0)\\\
\varepsilon^{{}^{\prime\prime}}_{n}(0)\end{array}\right).$ (24)

Here, the period $T=\pi/\omega$ is only half the size of that of the periodic
solutions $s^{(0)}_{n}(t)$. In addition, the complex nature of
$\varepsilon_{n}$ has to be taken into account to separate it into the real
$\varepsilon^{{}^{\prime}}_{n}$ and imaginary
$\varepsilon^{{}^{\prime\prime}}_{n}$ parts.

Figure 11:  Eigenvalues $\lambda$ of the monodromy matrix for the bottom dark
breather lattices of the size 101 with a different number of embedded kinks
$N=1,2,3$ in dependence on the magnetic field $\beta$. The parameters are
$B=0.15$, $\Omega=-0.28$, and $D/2J=0.16$. The critical magnetic field
$\beta_{\textrm{cr}}=1.325$.

Any solution can be determined by the column matrix of the initial conditions
$\Lambda(0)=\left[\varepsilon^{{}^{\prime}}_{1}(0),\varepsilon^{{}^{\prime\prime}}_{1}(0),\ldots,\varepsilon^{{}^{\prime}}_{L}(0),\varepsilon^{{}^{\prime\prime}}_{L}(0)\right]^{T}$.
The monodromy matrix can be easily constructed numerically by integrating
differential equations for $\varepsilon_{n}$ $2L$ times from $t=0$ to $T$ with
the initial conditions $\Lambda^{\nu}(0)$, $\nu=1,2,\ldots,2L$, with the
elements $\Lambda^{\nu}_{\mu}(0)=\delta_{\nu\mu}$. If the eigenvalues
$\lambda$ of the monodromy matrix $\hat{\mathcal{M}}$ lie on the unit circle
of the complex plane, then according to the Floquet theorem, the periodic
orbit is stable, otherwise it is unstable.

Performing these calculations we obtain that the dark breather solutions are
stable for any $\beta>\beta_{\textrm{cr}}=2\sqrt{1+q^{2}_{0}}+2B-1$
(hereinafter, $\beta=H/2JS$ is the magnetic field strength measured in the
$2JS$ units). Close to the threshold point the Floquet eigenvalues are located
almost uniform on the unit circle, while far from $\beta_{\textrm{cr}}$ they
are redistributed being mostly concentrated on the left half of the unit
circle, but stability nonetheless retains (Fig. 11).

Figure 12:  Dependence of the BL energy per bond $E_{b}$ on the number of
embedded kinks $N$ shown by dots for $\Omega=-0.18$ (blue), $\Omega=-0.22$
(brown), $\Omega=-0.28$ (green) and $\Omega=-0.32$ (red). Linear fitting to
the numerical data is shown by the dashed lines. The parameters are $B=0.15$,
$D/2J=0.16$ and $\beta=1.5$. The frequencies are restricted to the range
$\Omega(0)<\Omega<\Omega_{u}$ with $\Omega(0)=-0.325438$ and
$\Omega_{u}=-0.162719$ (see the text) and corresponds to the physical
frequencies $\omega$ as 1.485 THz ($\Omega=-0.18$), 1.440 THz
($\Omega=-0.22$), 1.373 THz ($\Omega=-0.28$) and 1.328 THz ($\Omega=-0.32$).
The chain length $L=101$.

An interesting feature of the bottom dark breather solutions is their energy
$E$ dependence on the number $N$ of embedded kinks. By rewriting the
Hamiltonian (1) in terms of the spin variables $s_{n}$ we get the form

$\frac{E}{2JS^{2}}=-\sqrt{1+q^{2}_{0}}\sum_{n=1}^{L-1}s_{n}s_{n+1}+B\sum_{n=1}^{L-1}\left(1-s^{2}_{n}\right)$
$-\sum_{n=1}^{L-1}\sqrt{1-s^{2}_{n}}\sqrt{1-s^{2}_{n+1}}-\beta\sum_{n=1}^{L-1}\sqrt{1-s^{2}_{n}}.$
(25)

With the help of the expression, the function $E_{\textrm{b}}$, the energy per
a lattice bond measured from the ground state value $E_{0}=-1+B-\beta$, is
plotted against $N$ in Fig. 12 for different values of the frequency $\Omega$
with $B$ and $\beta$ fixed.

To understand this result recall that the orientation of the moments
perpendicular to the chain axis is favorable from the viewpoint of easy-plane
anisotropy, however, there is a significant increase in the Zeeman energy.
With an addition of a kink/antikink the increase in the latter can be reduced
in the vicinity of a kink center, however, this would give rise to an increase
in the anisotropy energy. The competition of the two contributions, whereof
the Zeeman term dominates since the ratio $B/\beta$ is on the order of 0.1,
determines the details of dependence of $E_{\textrm{b}}$ on $N$, which is
projected on a straight line with negative slope (Fig. 12). The exchange
interaction plays no role owing to its full rotational invariance and the
almost parallel alignment of adjacent magnetic moments.

Being a continuous function of the frequency $\Omega$, the energy $E_{b}$
decreases gradually to zero when $\Omega$ varies from the upper boundary
$\Omega_{u}=1-B-\sqrt{1+q^{2}_{0}}$ to the bottom edge of the spin-wave
spectrum, $\Omega(0)=2-2B-2\sqrt{1+q^{2}_{0}}$. The requirement that $E_{b}$
should be nonnegative imposes a restriction on a number of kinks in the
breather lattice whereby the less breather frequency
$\omega=2JS(\Omega+\beta)/\hbar$, the less maximal number of kinks. It should
be mentioned at this point that a similar trend has been experimentally
observed in electrical lattices, where discrete breathers are produced via the
well-known mechanism of modulation instability English2010 ; Palmero2011 .
Numerical estimations for the physical frequency $\omega$ are given in Fig. 12
for $J=18$ K and $S=3/2$ relevant for CrNb3S6 Shinozaki2016 . Furthermore,
according to our calculations, the amplitude of the breather lattice modes
becomes larger for the modes of higher frequency but with the same number of
kinks. By definition, it corresponds to the regime of hard nonlinearity.

Finally, we point out that each of the discrete BL solutions possesses the
topological charge $\mathcal{Q}={\Delta\varphi}/{2\pi}={kL_{0}^{(d)}}/{2\pi}$,
where $\Delta\varphi$ is the total angle of spin moment rotation around the
$z$-axis per the BL period $L_{d}=4K\sqrt{(1+\kappa^{2})/|\alpha|}$. This
quantity includes the wave vector $k=-a^{-1}\tan^{-1}q_{0}$ resulting from the
DM interaction. The same factor appears in the expression for the linear
momentum density $\mathcal{P}_{n}=\hbar kS\left(1-\sqrt{1-s^{2}_{n}}\right)$ ,
which takes nonzero values whenever the transverse spin accumulation $s_{n}$
emerges. Thus, the important consequences of the DM interaction are
topological protection ensured by the topological charge and a possibility of
sliding motion, which requires no expense of energy, provided by the nonzero
linear momentum.

##  V Results and Discussion.

Noncollinear spin textures in chiral helimagnets and their concomitant
fundamental excitations have aroused research interest as technologically
relevant objects for spintronic applications Gobel2021 ; Inoue2021 ; Yang2020
. Unfortunately, nonlinear breather excitations lie outside the mainstream of
these investigations and the number of relevant works is very limited
Kiselev2012 ; Kiselev2013 . In our study we argue that this kind excitations
may arise in the form of magnetic discrete breathers in the monoaxial chiral
helimagnet CrNb3S6 owing to fairly small easy-plane anisotropy and the
specific domain structure of the compound. The latter has strong impacts on
standing spin waves in the thin films of CrNb3S6 Kishine2019 . Namely the
smallness of microdomains of definite crystallographic chirality necessitates
discrete approach to studying the breather modes.

We treat the phase of the forced ferromagnetic order, which is established
above the threshold value of the external magnetic field directed along the
chiral axis, and classify all types of possible solutions permitted by the
model Hamiltonian of the monoaxial chiral helimagnet. Using the data on
magnetic anisotropy in CrNb3S6, we estimate the relevance of these excitations
and conclude that only the dark breather modes, having the form of the
periodic breather lattice, are supported just above the bottom edge of the
spin wave spectrum. This kind of nonlinear solutions, which is proved to be
stable, supplements diversity of nontrivial spin structures in chiral
helimagnets and may contribute to spin response in ESR experiments.

The proposed classification scheme of DB solutions within the continuum
approach offers a natural interpretation of findings of the early numerical
studies on intrinsic localized spin modes in ferromagnetic Heisenberg chains
Rakhmanova1996 ; Rakhmanova1998 . It should be pointed out, however, that the
inclusion of the Dzyaloshinskii-Moryia interaction is not simply a
consideration for an additional form of magnetic anisotropy, the antisymmetric
exchange provides topological charge and linear momentum for emerging breather
modes.

In accordance with a generally accepted concept, DB modes owe their existence
to two pivotal components, discreteness and nonlinearity. In the case of
CrNb3S6, the discrete nature becomes evident by comparing the characteristic
domain size $L\sim$ 1 $\mu\textrm{m}$ with the BL period $L_{0}$. Our
calculations show (see, for instance, Fig. 10) that the period encompasses
around 100 sites that yields $L_{0}\sim$ 0.1 $\mu\textrm{m}$, i.e. only one
tenth of $L$, if the distance 6.847 $\AA$ between the nearest Cr ions along
the chiral axis Mandrus2013  is accounted for. Nonlinearity of our problem is
guaranteed by the single ion anisotropy $B$ and the constraint that the spin
length is conserved in dynamical processes.

The dark breather lattices may be classified by the number $N$ of embedded
kinks/antikinks, but modes with the same $N$ may appear at different
frequencies going by the rule whereby the higher frequency, the bigger
amplitude. This fact imposes severe restriction upon a maximum value of
kinks/antikins that can be accommodated in a finite length segment. A striking
feature of the dark breather lattice is its energy which decreases linearly
with increasing $N$ in contrast to linear growth for bright modes Bostrem2021
. This strange behavior reflects a competition in the phase of forced
ferromagnetism between the single-ion anisotropy, the antisymmetric exchange
and the Zeeman coupling whereof the last interaction dominates and what
specifically makes the phase stable.

These findings mean that the micron-sized domain of CrNb3S6, in which the DB
lattice is excited, is a system where areas of stored energy are linked to
regions of transverse spin accumulation. As such, there arises a strictly 1D
periodic array of resonators whose stored energy can take only discrete values
and is controlled by two degrees of freedom, namely, the kink/antikink number
$N$ (or equivalently, by the BL period) and the BL frequency $\omega$ (or
equivalently, by the BL amplitude). According to our order-of-magnitude
estimate, $\omega$ lies in the THz range. Because of nonlinearity, an amount
of stored energy in these ”breather capacitors” will be far beyond the
capacity of standing spin waves. This functionality may be used to design
spintronic resonators on the base of chiral helimagnets. Care should be taken
to exclude resonances with linear spin waves, since frequencies of the dark
breather modes lie inside the SW band. Whilst discreteness of the system
ensures that these resonances may be avoided if the breather mode frequency is
appropriately chosen, this issue requires more detailed examination in the
framework of the special theory Kopidakis1999 ; Kopidakis2000 . Another open
issue relates to a way to excite in controllable manner these breather modes
by external sources. The creation of special microfabricated microwave
antennae, like those used in a noncentrosymmetrical ferromagnet LiFe5O8 to
induce magnons with large momentum, appears to be a promising way Iguchi2015 .
In this case, microwave wavelength emitted by the antenna must be matched with
a period of the excited BL mode. These challenges posed by the practical
application together with the basic problem of finding similar discrete
breather modes in the another ordered phases of the monoaxial chiral
helimagnet, i.e. conical and soliton lattice ones, should be addressed in
future works.

##  VI Conclusions.

In summary, we investigate the possibility of experimental observation of
discrete breather excitations in the forced ferromagnetic phase of the
monoaxial chiral helimagnet CrNb3S6. We found that the presence of the DM
interaction does not prevent their emergence but, conversely, ensures their
topological protection. We demonstrate that the model Hamiltonian of the
monoaxial chiral helimagnet supports several kinds of possible spatially
periodic discrete breather solutions, but a specific choice depends on subtle
interplay between interactions accounted for in the Hamiltonian. Using
experimental data for strengths of the DM coupling and the easy-plane single-
ion anisotropy in CrNb3S6 we predict for this magnetic material the appearance
of the dark breather lattice modes with frequencies lying within the linear
spin-wave spectrum near to its bottom edge. Results of numerical simulations
of the DB excitations on a finite-length chain are in very good agreement with
the continuum approximation, which relates the DB lattice modes and periodic
solutions of the nonlinear Duffing equation. The Floquet analysis confirms
stability of the discrete breather excitations relevant for CrNb3S6. Their
classification based on a number of embedded kinks/antikinks is suggested, and
dependence of their energy inversely proportional to this number is
established. This unusual feature of the dark breather lattice modes
positioned just above the spin wave band bottom paves a new way to design
spintronic capacitors based on chiral helimagnets and operating on the
principle energy pumping from the DB modes to the magnon sector.

Figure 13:  Phase portraits of the Duffing equation (6). (a) The bright
solutions at $\alpha=2$ and $\beta=1$; (b) the dark solutions at $\alpha=-2$
and $\beta=-1$.

###### Acknowledgements.

This work was supported by the Act of the Government of the Russian Federation
(contract No. 02.A03.21.0006). I.G.B., E.G.E. and V.E.S. acknowledge financial
support by the Russian Foundation for Basic Research (RFBR), Grant No.
20-02-00213. A.S.O. thanks the Russian Foundation for Basic Research (RFBR),
Grant No. 20-52-50005, and the Ministry of Science and Higher Education of the
Russian Federation, project No. FEUZ-2020-0054. J.K. acknowledges financial
support by JSPS KAKENHI Grant No. 17H02923.

##  Appendix A Duffing equation

In the longwave limit, when a characteristic spatial scale of DB solutions is
much larger than the lattice unit $a$, one may proceed with the continuous
coordinate $z=na$ instead of the discrete parameter $n$ in Eq. (3).

Approximating the spin variables by their continuous counterparts

$s_{n}=(-1)^{n}\psi(z),$ (26) $s_{n\pm 1}=(-1)^{n\pm 1}\left[\psi(z)\pm
a\frac{d\psi}{dz}+\frac{a^{2}}{2}\frac{d^{2}\psi}{dz^{2}}\right],$ (27)
$\sqrt{1-s^{2}_{n}}\approx 1-\frac{1}{2}\psi^{2}(z),$ (28)
$\sqrt{1-s^{2}_{n\pm 1}}\approx 1-\frac{1}{2}\psi^{2}(z)\mp
a\psi(z)\frac{d\psi(z)}{dz},$ (29)

and substituting them into Eq. (3) we obtain eventually the classical Duffing
equation (6) if nonlinear terms involving $d\psi/dz$ and $d^{2}\psi/dz^{2}$
are neglected.

Behavior of restricted solutions of the Duffing equation (6) is clearly
illustrated by their trajectories on the phase plane $(\psi,d\psi/dz)$. A type
of the solutions depends on the signs of the coefficients $\alpha$ and $\beta$
(see, for example, Ref. Ostrovsky2002 ).

If $\alpha>0$, $\beta>0$ (Fig. 13a), the phase trajectories 1,2 and 3 conform
with the breather solutions (8), (13) and (11), respectively. This is the case
of the bright modes.

If $\alpha<0$, $\beta<0$ (Fig. 13b), the phase trajectories 1 and 2 correspond
to (14) and (16), respectively. This is the case of the dark modes.

## References

  * (1) A.M. Kosevich, B.A. Ivanov, and A.S. Kovalev, Phys. Rep. 194, 117 (1990).
  * (2) E. Infeld, G. Rowlands, Nonlinear waves, solitons and chaos (Cambridge University Press, Cambridge, 1990).
  * (3) S. Flach and C.R. Willis, Phys. Rep. 295, 181 (1998).
  * (4) S. Flach and A.V. Gorbach, Phys. Rep. 467, 1 (2008).
  * (5) S.V. Dmitriev, E.A. Korznikova, Yu.A. Baimova, M.G. Velarde, Phys. Usp. 59, 446 (2016).
  * (6) S.R. Bickham, S.A. Kiselev, A.J. Sievers, in Spectroscopy and Dynamics of Collective Excitations in Solids, NATO ASI Series B (Springer, Boston, MA, 1997), Vol. 356, p. 417.
  * (7) B.I. Swanson, J.A. Brozik, S.P. Love, G.F. Strouse, and A.P. Shreve, A. R. Bishop, W.-Z. Wang, and M.I. Salkola, Phys. Rev. Lett. 82, 3288 (1999).
  * (8) M.E. Manley, O. Hellman, N. Shulumba, A.F. May, P.J. Stonaha, J.W. Lynn, V.O. Garlea, A. Alatas, R.P. Hermann, J.D. Budai, H. Wang, B.C. Sales, and A.J. Minnich, Nat. Commun. 10, 1928 (2019).
  * (9) F. M. Russell and J. C. Eilbeck, Europhys. Lett. 78, 10004 (2007).
  * (10) P. Binder, D. Abraimov, A.V. Ustinov, S. Flach, and Y. Zolotaryuk, Phys. Rev. Lett. 84, 745 (2000).
  * (11) E. Trías, J. J. Mazo, and T. P. Orlando, Phys. Rev. Lett. 84, 741 (2000).
  * (12) J. Cuevas, L. Q. English, P. G. Kevrekidis, and M. Anderson, Phys. Rev. Lett. 102, 224101 (2009).
  * (13) M. Sato, B. E. Hubbard, L. Q. English, A. J. Sievers, B. Ilic, D. A. Czaplewski and H. G. Craighead, Chaos 13, 702 (2003).
  * (14) H. S. Eisenberg, Y. Silberberg, R. Morandotti, A. R. Boyd, and J. S. Aitchison, Phys. Rev. Lett. 81 3383 (1998).
  * (15) N. Boechler, G. Theocharis, S. Job, P. G. Kevrekidis, Mason A. Porter, and C. Daraio, Phys. Rev. Lett. 104, 244302 (2010).
  * (16) K. Narahara, J. Phys. Soc. Jpn. 89, 074005 (2020).
  * (17) R. Lai and A. J. Sievers, Phys. Rep. 314, 147 (1999).
  * (18) U. T. Schwarz, L. Q. English, and A. J. Sievers, Phys. Rev. Lett. 83, 223 (1999).
  * (19) M. Sato and A. J. Sievers, Nature 432, 486 (2004).
  * (20) R. F. Wallis, D. L. Mills and A. D. Boardman, Phys. Rev. B 52, R3828 (1995).
  * (21) S. Rakhmanova and D. L. Mills, Phys. Rev. B 54, 9225 (1996).
  * (22) Y.Zolotaryuk, S. Flach and V. Fleurov, Phys. Rev. B 63, 214422 (2001).
  * (23) J. M. Khalack, Y. Zolotaryuk and P. L. Christiansen, Chaos 13, 683 (2003).
  * (24) M. Lakshmanan, B. Subash and A. Saxena, Phys. Lett. A 378, 1119 (2014).
  * (25) M. Lakshmanan and A. Saxena, Phys. Lett. A 382, 1890 (2018).
  * (26) R. S. Kamburova, S. K. Varbev and M. T. Primatarowa, Phys. Lett. A 383, 471 (2019).
  * (27) Z. I. Djoufack, E. Tala-Tebue, J. P. Nguenang and A. Kenfack-Jiotsa, Chaos 26, 103110 (2016).
  * (28) L. Kavitha, E. Parasuraman, D. Gopi, A. Prabhu and R. A. Vicencio, J. Magn. Magn. Mat. 401, 394 (2016).
  * (29) R.L. Pylypchuk, and Y. Zolotaryuk, Low Temp. Phys. 41, 733 (2015).
  * (30) A. M. Kosevich and A. S. Kovalev, Sov. Phys. JETP 67, 1793 (1974).
  * (31) B. A. Ivanov, A. M. Kosevich, and I. M. Babich, JETP Lett. 29, 714 (1979).
  * (32) V. G. Bar’yakhtar, M. V. Chetkin, B. A. Ivanov, and S. N. Gadetskii, Dynamics of Topological Magnetic Solitons: Experiment and Theory (Springer Tracts in Modern Physics) (Springer, New York, 1994).
  * (33) N. Papanicolaou, Phys. Rev. B 55, 12290 (1997).
  * (34) A. V. Gorbach and S. Flach, Phys. Rev. E 72, 056607 (2005).
  * (35) X. Quan, and T. Qiang, Commun. Theor. Phys. 51, 153 (2009).
  * (36) P.F. Byrd, M.D. Friedman, Handbook of Elliptic Integrals for Engineers and Physics (Springer-Verlag, Berlin, 1954).
  * (37) R. McLachlan, Math. Intelligencer 16, 31 (1994).
  * (38) P.G. Kevrekidis, A. Saxena, and A.R. Bishop, Phys. Rev. E 64, 026613 (2001).
  * (39) Z. Fu, S. Liu, and S. Liu, Phys. Scr. 76, 15 (2007).
  * (40) Z. Fu, S. Liu, Z. Naturforsch. 62a, 555 (2007).
  * (41) P.G. Kevrekidis, A. Khare, and A. Saxena, Phys. Rev. E 68, 047701 (2003).
  * (42) P.G. Kevrekidis, A. Khare, A. Saxena, and G. Herring, J. Phys. A: Math. Gen. 37, 10959 (2004).
  * (43) Z. Yan, Phys. Lett. A 372, 969 (2008).
  * (44) X. Zhao, Z.T. Fu, J.Y. Mao, S.K. Liu, Commun. Theor. Phys. 52, 23 (2009).
  * (45) Z. Fu, S. Liu, and S. Liu, Phys. Lett. A 368, 238 (2007).
  * (46) A.P. Tankeyev, V.V. Smagin, M.A. Borich, and A.S. Zhuravlev, Phys. Met. Metallogr. 107, 229 (2009).
  * (47) A.P. Tankeyev, M.A. Borich, and V.V. Smagin, Phys. Met. Metallogr. 110, 689 (2010).
  * (48) V.V. Smagin, A.P. Tankeyev, M.A. Borich, Phys. Met. Metallogr. 108, 425 (2009).
  * (49) F. Yin and B. Tang, Int. J. Theor. Phys. 57, 2843 (2018).
  * (50) X. Quan, and T. Qiang, Chin. Phys. Lett. 25, 2400 (2008).
  * (51) R. I. Babicheva, A. S. Semenov, E. G. Soboleva, A. A. Kudreyko, K. Zhou, and S. V. Dmitriev, Phys. Rev. E 103, 052202 (2021).
  * (52) S.V. Rakhmanova, and A.V. Shchegrov, Phys. Rev. B 57, R14012 (1998).
  * (53) A. M. Weiner, J. P. Heritage, R. J. Hawkins, R. N. Thurston, E. M. Kirschner, D. E. Leaird, and W. J. Tomlinson Phys. Rev. Lett. 61, 2445 (1988).
  * (54) C. Bao, Y. Xuan, C. Wang, A. Fülöp, D. E. Leaird, V. Torres-Company, M. Qi, and A. M. Weiner Phys. Rev. Lett. 121, 257401 (2018).
  * (55) Y. Togawa, T. Koyama, Y. Nishimori, Y. Matsumoto, S. McVitie, D. McGrouther, R. L. Stamps, Y. Kousaka, J. Akimitsu, S. Nishihara, K. Inoue, I. G. Bostrem, Vl. E. Sinitsyn, A. S. Ovchinnikov, and J. Kishine, Phys. Rev. B 92, 220412(R) (2015).
  * (56) B. Schmid, B. Dorner, D. Petitgrand, L.P. Regnault, and M. Steiner, Z. Phys. B 95, 13 (1994);
  * (57) J. Kishine and A. S. Ovchinnikov, Theory of Monoaxial Chiral Helimagnet in Solid State Physics, vol. 66, R. E. Camley and R. L. Stamps (Academic Press, New York, 2015).
  * (58) J. L. Marín, S. Aubry and L. M. Floría, Physica D 113, 283 (1998).
  * (59) J. F. R. Archilla, J. Cuevas, B. Sánchez-Rey, and A. Alvarez, Physica D 180, 235 (2003).
  * (60) L.Q. English, F. Palmero, A.J. Sievers, P.G. Kevrekidis, and D H. Barnak, Phys. Rev. E 81, 046605 (2010).
  * (61) F. Palmero, L.Q. English, J. Cuevas, R. Carretero-González, and P. G. Kevrekidis, Phys. Rev. E 84, 026605 (2011).
  * (62) M. Shinozaki, S. Hoshino, Y. Masaki, J. Kishine, and Y. Kato, J. Phys. Soc. Jpn. 85, 074710 (2016).
  * (63) T. Miyadai, K. Kikuchi, H. Kondo, S. Sakka, M. Arai and Y. Ishikawa, J. Phys. Soc. Jpn. 52, 1394 (1983).
  * (64) A.A. Anders, G.V. Borisenko, and S. V. Volotskii, Sov. J. Low Temp. Phys. 15, 21 (1989).
  * (65) I.G. Bostrem, Vl.E. Sinitsyn, A.S. Ovchinnikov, E.G. Ekomasov, and J. Kishine, AIP Advances 11, 015208 (2021).
  * (66) I. E. Dzyaloshinskii, Sov. Phys. JETP 19, 960 (1964); 20, 665 (1965).
  * (67) A. Sievers, S. Takeno, Phys. Rev. Lett. 61, 970 (1988).
  * (68) J. B. Page, Phys. Rev. B 41, 7835 (1990).
  * (69) P. Pincus, Phys. Rev. 118, 658 (1960).
  * (70) J. Kishine, Vl.E. Sinitsyn, I.G. Bostrem, Igor Proskurin, F.J.T. Goncalves, Y. Togawa, and A.S. Ovchinnikov, Phys. Rev. B 100, 024411 (2019).
  * (71) B. Göbel, I. Mertig, and O.A. Tretyakov, Phys. Rep. 895, 1 (2021).
  * (72) K. Inoue, Chem. Lett. 50, 742 (2021).
  * (73) S.-H. Yang, Appl. Phys. Lett. 116, 120502 (2020).
  * (74) C. Zhang, J. Zhang, C. Liu, S. Zhang, Y. Yuan, P. Li, Y. Wen, Z. Jiang, B. Zhou, Y. Lei, D. Zheng, C. Song, Z. Hou, W. Mi, U. Schwingenschlögl, A. Manchon, Z. Q. Qiu, H. N. Alshareef, Y. Peng, X.-X. Zhang, Adv. Mater. 2101131 (2021).
  * (75) V.V. Kiselev, and A.A. Raskovalov, Theor. Math. Phys. 173, 1565 (2012).
  * (76) V.V. Kiselev, and A.A. Raskovalov, J. Exp. Theor. Phys. 116, 272 (2013).
  * (77) N. J. Ghimire, M. A. McGuire, D. S. Parker, B. Sipos, S. Tang, J.-Q. Yan, B. C. Sales, and D. Mandrus, Phys. Rev. B 87, 104403 (2013).
  * (78) G. Kopidakis, and S. Aubry, Physica D 130, 155 (1999).
  * (79) G. Kopidakis, and S. Aubry, Phys. Rev. Lett. 84, 3236 (2000).
  * (80) Y. Iguchi, S. Uemura, K. Ueno, and Y. Onose, Phys. Rev. B 92, 184419 (2015).
  * (81) L.A. Ostrovsky, A.I. Potapov, Modulated Waves Theory and Application (Johns Hopkins University Press, Baltimore, 2002).