# Unbiased and affordable atomic forces in _ab initio_ Variational Monte Carlo

Kousuke Nakano<EMAIL_ADDRESS>Center for Basic Research on
Materials, National Institute for Materials Science (NIMS), Tsukuba, Ibaraki
305-0047, Japan  Michele Casula  Institut de Minéralogie, de Physique des
Matériaux et de Cosmochimie (IMPMC), Sorbonne Université, CNRS UMR 7590, IRD
UMR 206, MNHN, 4 Place Jussieu, 75252 Paris, France  Giacomo Tenti
International School for Advanced Studies (SISSA), Via Bonomea 265, 34136,
Trieste, Italy

###### Abstract

_Ab initio_ quantum Monte Carlo (QMC) is a state-of-the-art numerical approach
for evaluating accurate expectation values of many-body wavefunctions.
However, one of the major drawbacks that still hinders widespread QMC
applications is the lack of an affordable scheme to compute unbiased atomic
forces. In this study, we propose a very efficient method to obtain unbiased
atomic forces and pressures in the Variational Monte Carlo (VMC) framework
with the Jastrow-correlated Slater determinant ansatz, exploiting the gauge-
invariant and locality properties of its geminal representation. We
demonstrate the effectiveness of our method for H2 and Cl2 molecules and for
the cubic boron nitride crystal. Our framework has a better algorithmic
scaling with the system size than the traditional finite-difference method,
and, in practical applications, is as efficient as single-point VMC
calculations. Thus, it paves the way to study dynamical properties of
materials, such as phonons, and is beneficial for pursuing more reliable
machine-learning interatomic potentials based on unbiased VMC forces.

##  Introduction $-$

_Ab initio_ quantum Monte Carlo (QMC) [1] is a state-of-the-art numerical
approach for evaluating the expectation values of many-body wavefunctions. It
usually provides extremely accurate energies. To date, QMC has been
successfully applied to various materials for which other electronic structure
methods, such as the Density Functional Theory (DFT), lose predictive power.
Examples are molecular crystals [2], two-dimensional materials [3, 4, 5],
superconductors [6], and materials at extreme pressures [7, 8, 9, 10, 11].
Despite several successful applications done so far and the recent development
of sophisticated QMC packages [12, 13, 14, 15, 16], this technique is not as
widely used as other established electronic structure methods. One of the main
QMC drawbacks is the lack of an efficient and affordable scheme to compute
unbiased atomic forces consistent with the derivatives of the total energy
with respect to atomic positions. In contrast, within DFT, atomic forces are
readily evaluated by exploiting the Hellmann-Feynman theorem [17] and can be
used to study dynamical properties of materials, such as phonons [18, 19].
This problem is also relevant in the construction of machine learning
potentials (MLPs), which need large datasets, where energy and forces are
computed with the method of choice. Recently, some QMC-driven MLPs have been
reported [20, 21, 22], where the availability of unbiased forces and pressures
has been a major concern.

There are two main real-space QMC frameworks, the variational Monte Carlo
(VMC) and the fixed-node diffusion Monte Carlo (FN-DMC) methods [1]. In this
study, we focus on VMC because the forces computation within the FN-DMC
framework is much more difficult and it is still a highly debated topic [23,
24, 25, 26, 27, 28, 29, 30]. Let ${\bf R}_{\alpha}$ be the atomic position of
the nucleus $\alpha$. The atomic force acting on $\alpha$ is defined as the
negative gradient of the energy with respect to ${\bf R}_{\alpha}$:

$\displaystyle{\bf F}_{\alpha}=-\frac{dE}{d{\bf R}_{\alpha}}=$
$\displaystyle-\braket{\frac{\partial}{\partial{\bf R}_{\alpha}}E_{\rm L}}$
(1a) $\displaystyle-2\braket{(E_{\rm L}-E)\frac{\partial\log\Psi_{\rm
T}}{\partial{\bf R}_{\alpha}}}$ (1b)
$\displaystyle-\braket{\sum_{i=1}^{N_{p}}{\frac{\partial E}{\partial
p_{i}}}{\frac{dp_{i}}{d{\bf R}_{\alpha}}}},$ (1c)

where $\Psi_{\rm T}$ is the variational wavefunction, $\braket{A}$ indicates
the quantum average of the local operator $A$ over the VMC sampling of
$|\Psi_{\rm T}|^{2}$, $E_{\rm L}$ is the so-called local energy ($E_{\rm
L}\equiv\hat{H}\Psi_{\rm T}/\Psi_{\rm T}$), with $E\equiv\braket{E_{L}}$, and
$\\{p_{1},\cdots,p_{N_{p}}\\}$ is the set of $N_{p}$ variational parameters
included in the $\Psi_{\rm T}$ ansatz. Eqs. (1a), (1b), and (1c) are called
the Hellmann–Feynman (HF), Pulay, and variational terms, respectively. One
usually ignores Eq. (1c) when evaluating atomic VMC forces, resulting in

${\bf F}_{\alpha}^{\rm{VMC}}=-\braket{\frac{\partial}{\partial{\bf
R}_{\alpha}}E_{\rm L}}-2\braket{(E_{\rm L}-E)\frac{\partial\log\Psi_{\rm
T}}{\partial{\bf R}_{\alpha}}}.$

The long-standing problem of obtaining a statistically meaningful ${\bf
F}_{\alpha}^{\rm{VMC}}$ value with a finite variance and at the same cost as
the VMC energy evaluation has been solved by the zero-variance zero-bias
principle [31] together with the so-called reweighting technique [32, 24, 31,
33, 34, 35, 36]. Hereafter, we will denote ${\bf F}_{\alpha}^{\rm{VMC}}$ as
regular VMC force.

Neglecting Eq. (1c) is justified only when the system is at its variational
minimum for all parameters (i.e., $\partial E/\partial p_{i}=0$, $\forall i$)
or when the variational parameters, which are implicitly dependent on the
atomic positions, accidentally or by construction become position-independent
(i.e., $dp_{i}/d{\bf R}_{\alpha}=0$, $\forall i$); otherwise, ${\bf
F}_{\alpha}^{\rm{VMC}}$ can be biased. This bias is referred to as self-
consistency error [37, 36].

We note that the traditional finite-difference method (FDM), where the atomic
forces are evaluated as ${\bf F}_{\alpha}^{\rm{FDM}}\sim-\left(E^{{{\bf
R}_{\alpha}}+\Delta{{\bf R}_{\alpha}}}-E^{{{{\bf R}_{\alpha}}-\Delta{{\bf
R}_{\alpha}}}}\right)/2\Delta{{\bf R}_{\alpha}}$, includes in principle all
terms. Thus, it can provide unbiased forces for small enough $\Delta{{\bf
R}_{\alpha}}$, but it becomes soon unfeasible as it requires an $O(N)$ larger
computational cost than standard VMC calculations, with $N$ the number of
atoms in the system.

In this Letter, we propose a method to obtain unbiased atomic forces and
pressures that does not increase the computational complexity of the VMC
energy calculation, by supplementing the regular VMC force with a suitable
variational term, computed by exploiting the gauge-invariant and locality
properties of the Antisymmetrized Geminal Power (AGP) ansatz [38]. For
assessment, we demonstrate that the potential energy surfaces (PESs) of the H2
and Cl2 molecules, and the equation of state (EOS) of the cubic Boron Nitride
(cBN) are consistent with the forces and pressure obtained by our proposed
method.

##  Illustrating the problem $-$

For the sake of clarity, we present the case of the PES of a dimer expressed
as a function of the interatomic distance $R$, while the present discussion
can be applied for any other system. Fig. 1 shows a schematic picture of
several PESs. Let $E^{\rm true}$(a) be the true PES of the dimer. $E^{\rm
true}$ is the ultimate goal of any electronic structure calculation, but it is
unknown except for hydrogen-like dimers, or more generally, for nodeless
ground states. The best possible PES within a given $\Psi_{\rm T}$ ansatz is
$E^{\rm fullopt}$(b), yielded by a VMC calculation with the fully optimized
$\Psi_{\rm T}$. This is achievable for rather small systems by optimization
methods suitable for noisy data [39, 40]. However, the wavefunction
optimization becomes impractical for a large number of variational parameters.
Therefore, a good compromise between accuracy and computational efficiency is
$E_{\rm JSD}^{\rm true}$(c), obtained by the Jastrow correlated Slater
determinant (JSD) ansatz, with one-body molecular orbitals (MOs) computed by
DFT for each interatomic distance $R$. The JSD is the most common VMC ansatz:
only the Jastrow factor is optimized at the VMC level, while the DFT MOs are
kept frozen in the Slater determinant (SD). However, in this case, the VMC
force $F^{\rm{VMC}}$ is not consistent with the slope of $E_{\rm JSD}^{\rm
true}$, because the variational parameters included in the SD are not at their
VMC minima. Instead, $F^{\rm{VMC}}$ corresponds to the slope of $E_{\rm
JSD}^{\rm biased}$(d), where the DFT MOs obtained at $R=R^{\prime}$ are used
artificially for all $R$, such that $dp_{i}/dR=0$, $\forall i$. In this work,
we propose an efficient method to obtain atomic forces and pressures that are
unbiased, namely consistent with the slope of $E_{\rm JSD}^{\rm true}$(c).

Figure 1: Schematic picture of PESs as a function of the dimer bond length
$R$. (a) The true PES, not accessible in practice. (b) The best possible PES
obtained in the VMC framework by minimizing all variational parameters of
$\Psi_{\rm T}$. (c) The PES obtained with optimized Jastrow factor and Slater
MOs yielded by DFT at each point $R$, whose slope at $R^{\prime}$ is exactly
given by the VMC force $F^{\rm{VMC}}$ supplemented by the variational term
$F^{c}$, as proposed in this study. (d) The PES obtained with frozen DFT
orbitals computed at its minimum $R^{\prime}$, whose slope corresponds to
$F^{\rm{VMC}}$ without the additional term $F^{c}$. Figure 2: (a) and (b): H2
PESs (solid green curves), their numerical derivatives (dashed green curves),
regular VMC forces (dashed red curves), and corrected forces (dashed purple
curves) obtained with (a) small [1s] and (b) large [4s2p1d] Jastrow basis
sets. The PESs are computed from 0.30 Å to 2.00 Å with 18 equally spaced
datapoints plus 5 additional datapoints (0.55 Å, 0.65 Å, 0.75 Å, 0.85Å, and
0.95 Å and fitted by a polynomial of 11th (10th) order for the energies
(forces). (c): Cl2 PES (solid green curve), its numerical derivative (dashed
green curve), regular VMC force (dashed red curve), and corrected force
(dashed purple curve). The PES is computed from 1.50 to 2.80 Å with 14 equally
spaced datapoints and fitted by a polynomial of 6th (5th) order for energies
(forces). In all panels, only the region in the vicinity of the equilibrium
geometry is drawn.

As long as the Jastrow factor is at its variational minimum, the contribution
to the bias comes only from the SD part. This suggests that a straightforward
solution for correcting the bias is to compute the variational term ${\bf
F}_{\alpha}^{c}\equiv-\braket{\sum_{i=1}^{N_{p}^{\rm SD}}{\cfrac{\partial
E}{\partial p_{i}^{\rm SD}}}{\cfrac{dp_{i}^{\rm SD}}{d{\bf R}_{\alpha}}}}$,
where $p_{i}^{\rm SD}$ are SD variational parameters. In the following, we
explain the method to evaluate these terms by combining DFT with VMC gradients
calculations.

##  Method to obtain unbiased atomic forces $-$

We begin by introducing the AGP representation [38] of the SD ansatz made of
MOs. The general AGP ansatz for a system of $N_{\rm{e}}$ electrons is written
as $\Psi_{\rm AGP}=\hat{\cal A}[g({\bf x}_{1},{\bf x}_{2})g({\bf x}_{3},{\bf
x}_{4})\ldots g({\bf x}_{N_{\rm{e}}-1},{\bf x}_{N_{\rm{e}}})],$ where
$\hat{\cal A}$ is the antisymmetrization operator and $g$ is the so-called
geminal function $g({\bf x}_{l},{\bf x}_{m})=f({\bf r}_{l},{\bf
r}_{m})(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow})/\sqrt{2}$. The
spatial part $f({\bf r},\mathbf{r}^{\prime})$ can be written in terms of MOs,
such that $f({\bf r},\mathbf{r}^{\prime})=\sum_{k}^{M}\Phi_{k}\,({\bf
r})\Phi_{k}({\bf r}^{\prime})$, where $\Phi_{k}(x)$ is the $k$-th MO expressed
as $\Phi_{k}({\bf r})=\sum_{i}^{L}c_{i,k}\psi_{i}({\bf r})$, $\psi_{i}({\bf
r})$ is the $i$-th atomic orbital (AO), and, $c_{i,k}$ are the AO coefficients
obtained by a DFT calculation. If the Hilbert space is restricted to the
occupied states, i.e. $M=N_{e}/2$ for spin-unpolarized systems 111The AGP
ansatz has been generalized also to spin-polarized systems [38]., the
resultant AGP is equivalent to the SD ansatz. In other words, the SD ansatz
can be treated as a special case of the more general AGP wavefunction. We
assume that $\Psi_{\rm T}$ is real for the sake of conciseness; thus, the
variational parameters are also real. However, our method can be readily
generalized to complex $\Psi_{\rm T}$ [42]. In this work, the geminal function
is constructed from the MOs obtained from a DFT calculation, and then
converted to the AO representation, namely $f({\bf r}_{l},{\bf
r}_{m})=\sum_{i,j}^{L,L}\lambda_{i,j}\psi_{i}({\bf r}_{l})\psi_{j}({\bf
r}_{m}),$ with $\lambda_{i,j}=\sum_{k}c_{i,k}c_{j,k}$
222$\lambda_{i,j}=\sum_{k}c_{i,k}^{*}c_{j,k}$ in complex cases. Thus, the
variational term needed for correcting the self-consistency error reads

${\bf F}_{\alpha}^{\rm c}=-\braket{\sum_{i,j}^{L,L}{\cfrac{\partial
E}{\partial\lambda_{i,j}}}{\cfrac{d\lambda_{i,j}}{d{\bf R}_{\alpha}}}},$ (2)

which is what one should compute to get unbiased atomic forces in the JSD
ansatz.

The first factor in the brakets of Eq. 2, $\partial E/\partial\lambda_{i,j}$,
used for optimizing $\Psi_{\rm T}$ and often dubbed as generalized force, can
be efficiently computed by VMC [42]. The second factor, the total derivative
$d\lambda_{i,j}/d{{\bf R}_{\alpha}}$, can be numerically evaluated using the
FDM, i.e. $d\lambda_{i,j}/d{{\bf R}_{\alpha}}\sim\left(\lambda_{i,j}^{{{\bf
R}_{\alpha}}+\Delta{{\bf R}_{\alpha}}}-\lambda_{i,j}^{{{{\bf
R}_{\alpha}}-\Delta{{\bf R}_{\alpha}}}}\right)/2\Delta{{\bf R}_{\alpha}}$. We
notice that the latter term is $N$ times more time-consuming than the single-
point DFT calculation, and this is regardless of the number of variational
parameters. Indeed, to compute the correction terms for the geometry
$\mathbf{R}\equiv\\{{\bf R}_{1},\ldots,{\bf R}_{N}\\}$, one only needs
6$N$-times DFT calculations, where 6 is twice the number of Cartesian
components.

The gauge-invariant property of the AGP, inherited from its close relation
with the reduced one-body density matrix [44], allows one to construct a
robust workflow to compute the second term using finite differences. Indeed,
thanks to the gauge invariance of the AGP, one does not suffer from (i) the
global phase (or sign) indetermination of MOs, nor from (ii) their possible
degeneracy. As for (i), a global phase $\theta$ rotating the $k$-th MO
($\Phi_{k}\rightarrow e^{i\theta}\Phi_{k}$), which is reduced to a global sign
$e^{i\theta}=\pm 1$ in case of a real $\Psi_{\rm T}$, does not affect the
total energy, but prevents the calculation of orbital derivatives,
$dc_{i,k}/d{{\bf R}_{\alpha}}$, based on finite differences. Indeed, the
global phase (or sign) is sometimes inconsistent between DFT outcomes with
different atomic displacements. Instead, the sign flip is not problematic in
the AGP representation, because the relation
$\lambda_{i,j}=\sum_{k}c_{i,k}c_{j,k}$ [41] implies that $\lambda_{i,j}$ is
invariant under an MO sign change. As for (ii), when two MOs are degenerate,
e.g., $\Phi_{k^{\prime}}$ and $\Phi_{k^{\prime\prime}}$,
$c_{i,k^{\prime}}^{{{\bf R}_{\alpha}}+\Delta{{\bf R}_{\alpha}}}$ and
$c_{i,k^{\prime}}^{{{\bf R}_{\alpha}}-\Delta{{\bf R}_{\alpha}}}$ might have
very different values due to the presence of the other degenerate MO,
$\Phi_{k^{\prime\prime}}$. Nevertheless, it is straightforward to show that an
MOs degeneracy does not affect the uniqueness of $\lambda_{i,j}$, making
$\lambda_{i,j}$ independent of the choice of the particular DFT implementation
for degenerate MOs. Thus, by exploiting the AO representation of the AGP
wavefunction, one can always devise a well-defined method to compute
$d\lambda_{i,j}/d{{\bf R}_{\alpha}}$, which will be superior to the
calculation of $dc_{i,k}/d{{\bf R}_{\alpha}}$.

By combining the first and second factor in the brackets of Eq. 2, the
variational term can be cast in a form suitable for a VMC estimate, as follows
333Only its real part is taken in complex cases [42]:

$\displaystyle{\bf{F}}_{\alpha}^{\rm c}=-2\left\langle{E_{\rm
L}(x)\sum_{i,j}^{L,L}\left[(O_{i,j}(x)-\bar{O}_{i,j}){\cfrac{d\lambda_{i,j}}{d{\bf
R}_{\alpha}}}\right]}\right\rangle,$ (3)

where we made apparent the dependence of the local operators on the total
electronic coordinate $x$, sampled by VMC, to distinguish them from constant
values. In Eq. 3, $O_{i,j}(x)=\partial\ln\Psi_{T}(x)/\partial\lambda_{i,j}$,
and $\bar{O}_{i,j}\sim\braket{O_{i,j}(x)}$. We remark that ${O}_{i,j}(x)$ can
be efficiently computed in a VMC calculation using the adjoint algorithmic
differentiation [34], and the divergences of the generalized forces can be
cured by the reweighting method proposed by Attaccalite and Sorella [33] or
the one proposed by Pathak and Wagner [46]. It is extremely important that the
variational force is evaluated in a covariance form of random variables to
reduce its fluctuations [39, 47]. In addition, the expression in Eq. 3 implies
that if the variational wavefunction is an exact eigenstate of the
Hamiltonian, ${\bf{F}}_{\alpha}^{\rm c}$ vanishes regardless of the VMC
sample, because the local energy coincides with the corresponding eigenvalue
$E$. Indeed, the zero-variance property holds in this expression, which is
another way to recover the Hellmann-Feynman theorem.

##  Applications to H2 and Cl2 molecules $-$

We determine the interatomic force of the H2 and Cl2 molecules, taken as first
examples to assess the accuracy of our method. The correlation-consistent
effective core potentials (ccECPs) [48, 49, 50, 51] accompanied with the cc-
pVDZ basis sets were employed for both molecules. For Cl2, the He-core ccECP
was employed. DFT-MOs were prepared by PySCF v2.0.1 [52, 53] with the LDA-PZ
exchange-correlation functional [54], and then converted to the TurboRVB
wavefunction format [55] using the TurboGenius package [56] via TREX-IO files
[57]. The inhomogeneous one-body, the two-body, and the three-body Jastrow
factors [55] were added to the SD with frozen DFT MOs and optimized using the
linear method implemented in TurboRVB [55]. The second term,
$d\lambda_{i,j}/d{{\bf R}_{\alpha}}$, was numerically evaluated using the
displacements $\Delta R=\pm 0.001$ Å along the molecular bond direction.

The simple H2 molecule highlights the importance of removing the self-
consistency error in the forces calculation by adding the variational force
term to the regular VMC expression. In H2, the JSD ansatz with DFT MOs is in
principle exact if the Jastrow factor is converged in the basis set. Indeed,
the wavefunction is nodeless, so the difficulty of finding the optimal
variational state can be fully transferred to the Jastrow factor
determination. Thus, H2 allows one to study different situations, from a poor
to a refined Jastrow factor. In this study, we examined a small [1s] and a
large [4s2p1d] basis set expansion, as a poor and refined Jastrow factor,
respectively. Fig. 2(a) shows the H2 PES yielded by the poor Jastrow factor.
It indicates that the DFT parameters are not optimal at the VMC level; thus,
the self-consistency error is present. The equilibrium distance obtained from
the PES (0.7344(2) Å) and the one from regular VMC forces (0.7392(1) Å) are
reported in Tab. 1. Fig. 2 and Tab. 1 demonstrate that the self-consistency
error is mitigated by the proposed force correction $F^{c}$, which gives a
bond distance of 0.7341(1) Å, compatible with the one derived from the PES.
Fig. 2(b) shows, instead, the H2 PES obtained with the refined Jastrow factor.
It indicates that the DFT MOs with the large [4s2p1d] Jastrow expansion are
optimal at the VMC level, i.e., the refined Jastrow factor compensates for the
DFT determinant. In this case, the self-consistency error is negligible, as
expected, since all variational parameters included in the ansatz are optimal.
Thus, the regular VMC force is already consistent with the derivative of the
PES, and the corresponding force correction eventually vanishes, as reported
in Tab. 1. The H2 example is illustrative of the capability of the variational
term in Eq. 2 to correct the force bias due not only to the frozen DFT MOs,
but also to an underconverged Jastrow factor.

Fig. 2 (c) shows the PES of the Cl2 molecule, as yielded by a [3s1p] Jastrow
basis set. Tab. 1 reports the equilibrium geometries obtained from the PES,
regular VMC force, and corrected force. The Figure and Table show that the
self-consistency error is more significant for Cl2 ($Z_{\rm eff}=15$) than H2
($Z_{\rm eff}=1$). This is consistent with the seminal work by Tiihonen et al.
[37], reporting that the self-consistency error increases with the effective
nuclear charge. Fig. 2 and Tab. 1 clearly illustrate that the proposed force
correction works also for heavier molecules, resulting in corrected forces and
equilibrium geometries fully consistent with the VMC PES.

Table 1:  The equilibrium bond distances $r_{\rm eq}$ (Å) of the H2 and Cl2 molecules obtained from the PESs, the regular VMC force, and the corrected force. The corresponding PESs are shown in Fig. 2.  Dimers | Source |  $r_{{\rm eq}}$ (Å)
---|---|---
H2 (Jas. [1s]) | PES | 0.7344(2)
VMC force | 0.7392(1)
Corrected force | 0.7341(1)
Experiment | 0.741 111These values are taken from Ref. 58.
H2 (Jas.[4s2p1d]) | PES | 0.7418(3)
VMC force | 0.7408(6)
Corrected force | 0.7408(6)
Experiment | 0.741 111These values are taken from Ref. 58.
Cl2 (Jas.[3s1p]) | PES | 1.987(1)
VMC force | 1.9979(1)
Corrected force | 1.9864(1)
Experiment | 1.987 111These values are taken from Ref. 58.
|  |

##  Application to cubic Boron Nitride $-$

Not only atomic forces, but also pressures can be corrected in solids using
the same method, just by replacing $d\lambda_{i,j}/d{{\bf R}_{\alpha}}$ with
$d\lambda_{i,j}/dV$. To demonstrate it, we computed the cBN EOS. The
correlation-consistent effective core potentials (ccECPs) [48, 49, 50, 51]
with accompanied cc-pVDZ basis sets were used for the cBN calculation. The
linear dependency of the basis sets is solved at the DFT level by cutting
basis set elements with exponents smaller than 0.20 a.u. This is crucial to
suppress the statistical errors on atomic forces and pressures for periodic
systems in QMC calculations [59]. The 2$\times$2$\times$2 conventional
supercell (256 valence electrons in the simulation cell) with $\mathbf{k}$ =
$\Gamma$ was employed. DFT-MOs were prepared by the built-in DFT module
implemented in TurboRVB [55] with the LDA-PZ exchange-correlation functional
[54]. Then, the inhomogeneous one-body, the two-body, and the three-body
Jastrow factors were added to the SD with frozen DFT molecular orbitals. A
[3s1p] Jastrow basis set was employed for both B and N atoms. The Jastrow
factor was optimized using the linear method implemented in TurboRVB [55] for
each volume. The second term, $d\lambda_{i,j}/dV$, was numerically evaluated
using the built-in DFT module with volume variations $\Delta V=\pm 0.1$ %.

Fig. 3 shows the cBN EOS, its volume derivative, the regular VMC pressures,
and the corrected pressures. The EOS and pressures were fitted by the Vinet
form [60]. The obtained equilibrium lattice parameters and volumes are
reported in Tab. 2. It is apparent that the self-consistency error in pressure
is $\sim 5$ GPa, constant over the whole volume range. Our method gives
corrections that bring the estimated pressures very close to the exact values
for all volumes, as shown in Fig. 3. This result illustrates the possibility
to successfully correct not only atomic forces but also pressures in large
systems with the explicit evaluation of the variational force term.

Figure 3: cBN EOS (solid green curve), its volume derivative (dashed green curve), the regular VMC pressure (dashed red curve), and the corrected pressure (dashed purple curve). The 11 datapoints are fitted by the Vinet EOS [60].  Table 2:  Equilibrium lattice parameters and volumes per atom obtained by fitting the EOS, and from the regular VMC pressure and the corrected one. Zero point energy and temperature effects are not included. Source | Lattice (Å) | Volume (Bohr3)
---|---|---
EOS | 3.5962(3) | 39.232(9)
VMC pressure | 3.5800(1) | 38.704(5)
Corrected pressure | 3.5943(2) | 39.169(5)
Experiment | 3.594 111These values are taken from Ref. 58. | 39.160 111These values are taken from Ref. 58.
|  |

##  Discussion about the computational cost $-$

We first compare our method with the FDM, which is the traditional way to
obtain unbiased atomic forces in the VMC framework. The main drawback of the
FDM is that it requires at least 6$N$ independent QMC runs to compute all 3$N$
force components, preventing its use in routine QMC calculations. Instead, our
proposed method requires just a single VMC run to compute all 3$N$ regular VMC
forces, together with all $O_{i,j}$ terms that appear in the expression for
${\bf{F}}_{\alpha}^{\rm c}$ of Eq. 3. This is thanks to the algorithmic
differentiation [34]. As we mentioned before, the other terms in Eq. 3, namely
$d\lambda_{i,j}/d{\bf R}_{\alpha}$, are computed by FDM using DFT as the
driver, thus leading to DFT calculations $N$-times more time-consuming than a
single DFT run. However, since the DFT cost is negligible compared to QMC, and
it is mainly fast Fourier transform bound, with a favorable $O(N^{2}\log N)$
scaling for a single run, the resulting algorithmic cost of our method is
superior to the FDM evaluation of atomic VMC forces.

Finally, we discuss the scaling of the variance of the variational term
${\bf{F}}_{\alpha}^{\rm c}$ with respect to $N$. The variance of the local
energy $E_{\rm L}$ scales with $N_{\rm e}$, since $E_{\rm L}$ is an extensive
property [42], while the variance of the logarithmic derivatives $O_{i,j}$ is
practically $O(1)$ [34]. Thus, the variance of ${\bf{F}}_{\alpha}^{\rm c}$ is
bound by $O(L^{2}{N_{\rm e}})$, where the factor of $L^{2}$ comes from the
double summation over the extended basis set elements. Since all terms in the
summation are evaluated by a single VMC run, covariance effects can then be
present. They are, however, neglected in the above variance estimation.
Nevertheless, in the actual calculations the statistical error of the sum of
the correlated random values was estimated by the Jackknife method [42]. The
variance can, in practice, be significantly reduced by neglecting
$\left|d\lambda_{i,j}/d{\bf R}_{\alpha}\right|$ with small absolute values,
which is justified because they are obtained deterministically by DFT
calculations. Thus, by exploiting the locality of the $\lambda_{i,j}$ matrix,
the summation in Eq. 3 can be reduced from $L^{2}$ to $L$ terms contributing
with non-negligible values in the $N\rightarrow\infty$ limit. In fact, the
percentages of disregarded elements according to the criterion
$\left|d\lambda_{i,j}/d{\bf R}_{\alpha}\right|/\max\left|d\lambda_{i,j}/d{\bf
R}_{\alpha}\right|\leq 0.01\%$ is 38.5 %, 45.0 %, and 66.0 % for 1 $\times$ 1
$\times$ 1 (8 atoms), 2 $\times$ 2 $\times$ 2 (64 atoms), and 3 $\times$ 3
$\times$ 3 (216 atoms) cBN supercells, respectively, demonstrating that the
larger the system, the more terms can be neglected. Since $L$ and $N_{\rm e}$
are proportional to $N$, the scaling of the variance of our method with
respect to $N$ is bound by $O(N^{2})$ in the $N\rightarrow\infty$ limit, which
is just $O(N)$ larger than the variance of the regular VMC force calculation,
$O(N)$ [34]. However, in our H2, Cl2, and cBN calculations, we got the same
error bars on the regular VMC forces and on the corrected ones with the same
statistics. This points to a very small prefactor in the $O(N^{2})$ variance
terms, such that in the total variance ${\rm Var}({\bf F}_{\alpha})={\rm
Var}({\bf F}_{\alpha}^{\rm{VMC}})+{\rm Var}({\bf F}_{\alpha}^{\rm{c}})\approx
O(N)+\varepsilon O(N^{2})$, and the $O(N^{2})$ contribution can be neglected
given a small $\varepsilon$ and an affordable $N$ for VMC calculations.

##  Concluding remarks $-$

In this work, we analyzed the bias seriously affecting the regular VMC
expression ${\bf F}_{\alpha}^{\rm{VMC}}$. We then proposed a method to
efficiently and robustly compute the missing contribution, i.e. the
variational term ${\bf F}_{\alpha}^{\rm{c}}$, to completely remove that bias
for a JSD ansatz with DFT one-body orbitals, the most common wavefunction in
_ab initio_ VMC calculations, usually the best compromise between accuracy and
computational cost. We demonstrated that the correction works very well for
the systems that have been tested here, namely the equilibrium geometry of H2
and Cl2 molecules, and the EOS evaluation of the cubic Boron Nitride. Unbiased
atomic forces within a JSD ansatz, which is in general much cheaper to
optimize than more refined $\Psi_{\rm T}$s, will be particularly useful to
generate VMC datasets for MLPs construction, which would otherwise be affected
by the self-consistency error. Thus, our approach has the potential to open up
new horizons for VMC applications, also in the context of machine learning.
Our proposed scheme can also be used to evaluate FN-DMC forces and pressures,
even if further work will be necessary to fully adapt the present scheme to
the DMC framework. Finally, the same scheme can be extended to more elaborated
wavefunctions, once a suitable auxiliary method is used to generate their
antisymmetric part.

##  Acknowledgements $-$

K.N. is grateful for computational resources from the Numerical Materials
Simulator at National Institute for Materials Science (NIMS). The authors are
grateful for computational resources of the supercomputer Fugaku provided by
RIKEN through the HPCI System Research Projects (Project IDs: hp230030). K.N.
acknowledges financial support from Grant-in-Aid for Early Career Scientists
(Grant No. JP21K17752), from Grant-in-Aid for Scientific Research (Grant No.
JP21K03400), and from MEXT Leading Initiative for Excellent Young Researchers
(Grant No. JPMXS0320220025). K.N. is grateful for fruitful discussion with A.
Togo (NIMS) about the illustration of the self-consistency problem. M.C.
acknowledges GENCI for granting access to French computational resources
through the DARI application 906493. The molecular and crystal structures were
depicted using VESTA [61]. This work is supported by the European Centre of
Excellence in Exascale Computing TREX - Targeting Real Chemical Accuracy at
the Exascale. This project has received funding from the European Union’s
Horizon 2020 - Research and Innovation program - under grant agreement no.
952165.

##  Code availability $-$

The QMC kernel used in this work, TurboRVB, is available from its GitHub
repository [https://github.com/sissaschool/turborvb].

## References

  * Foulkes  _et al._ [2001] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Zen  _et al._ [2018] A. Zen, J. G. Brandenburg, J. Klimeš, A. Tkatchenko, D. Alfè, and A. Michaelides, Proc. Natl. Acad. Sci. U.S.A. 115, 1724 (2018).
  * Mostaani  _et al._ [2015] E. Mostaani, N. D. Drummond, and V. I. Fal’ko, Phys. Rev. Lett. 115, 115501 (2015).
  * Frank  _et al._ [2019] T. Frank, R. Derian, K. Tokár, L. Mitas, J. Fabian, and I. Štich, Phys. Rev. X 9, 011018 (2019).
  * Nikaido  _et al._ [2022] Y. Nikaido, T. Ichibha, K. Hongo, F. A. Reboredo, K. H. Kumar, P. Mahadevan, R. Maezono, and K. Nakano, J. Phys. Chem. C 126, 6000 (2022).
  * Casula and Sorella [2013] M. Casula and S. Sorella, Phys. Rev. B 88, 155125 (2013).
  * Clay  _et al._ [2014] R. C. Clay, J. Mcminis, J. M. McMahon, C. Pierleoni, D. M. Ceperley, and M. A. Morales, Phys. Rev. B 89, 184106 (2014).
  * Clay  _et al._ [2016] R. C. Clay, M. Holzmann, D. M. Ceperley, and M. A. Morales, Phys. Rev. B 93, 035121 (2016).
  * Drummond  _et al._ [2015] N. D. Drummond, B. Monserrat, J. H. Lloyd-Williams, P. L. Ríos, C. J. Pickard, and R. J. Needs, Nat. Commun. 6, 1 (2015).
  * Mazzola  _et al._ [2018] G. Mazzola, R. Helled, and S. Sorella, Phys. Rev. Lett. 120, 025701 (2018).
  * Monacelli  _et al._ [2023] L. Monacelli, M. Casula, K. Nakano, S. Sorella, and F. Mauri, Nat. Phys. 19, 845 (2023).
  * Wagner  _et al._ [2009] L. K. Wagner, M. Bajdich, and L. Mitas, J. Comput. Phys. 228, 3390 (2009).
  * Scemama  _et al._ [2013] A. Scemama, M. Caffarel, E. Oseret, and W. Jalby, J. Comput. Chem. 34, 938 (2013).
  * Needs  _et al._ [2020] R. Needs, M. Towler, N. Drummond, P. Lopez Rios, and J. Trail, J. Chem. Phys. 152, 154106 (2020).
  * Kent  _et al._ [2020] P. R. C. Kent, A. Annaberdiyev, A. Benali, M. C. Bennett, E. J. Landinez Borda, P. Doak, H. Hao, K. D. Jordan, J. T. Krogel, I. Kylänpää, J. Lee, Y. Luo, F. D. Malone, C. A. Melton, L. Mitas, M. A. Morales, E. Neuscamman, F. A. Reboredo, B. Rubenstein, K. Saritas, S. Upadhyay, G. Wang, S. Zhang, and L. Zhao, J. Chem. Phys. 152 (2020).
  * Wheeler  _et al._ [2023] W. A. Wheeler, S. Pathak, K. G. Kleiner, S. Yuan, J. N. B. Rodrigues, C. Lorsung, K. Krongchon, Y. Chang, Y. Zhou, B. Busemeyer, K. T. Williams, A. Muñoz, C. Y. Chow, and L. K. Wagner, J. Chem. Phys. 158, 114801 (2023).
  * Martin [2004] R. M. Martin,  _Electronic structure : basic theory and practical methods_ (Cambridge University Press, 2004).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Togo [2023] A. Togo, J. Phys. Soc. Jpn. 92, 012001 (2023).
  * Tirelli  _et al._ [2022] A. Tirelli, G. Tenti, K. Nakano, and S. Sorella, Phys. Rev. B 106, L041105 (2022).
  * Niu  _et al._ [2023] H. Niu, Y. Yang, S. Jensen, M. Holzmann, C. Pierleoni, and D. M. Ceperley, Phys. Rev. Lett. 130, 076102 (2023).
  * Ceperley  _et al._ [2023] D. M. Ceperley, S. Jensen, Y. Yang, H. Niu, C. Pierleoni, and M. Holzmann, “Training models using forces computed by stochastic electronic structure methods,” (2023), arXiv:2310.15994 [cond-mat.mtrl-sci] .
  * Reynolds  _et al._ [1986] P. Reynolds, R. Barnett, B. Hammond, R. Grimes, and W. Lester Jr, Int. J. Quantum Chem. 29, 589 (1986).
  * Assaraf and Caffarel [2000] R. Assaraf and M. Caffarel, J. Chem. Phys. 113, 4028 (2000).
  * Filippi and Umrigar [2000] C. Filippi and C. J. Umrigar, Phys. Rev. B 61, R16291 (2000).
  * Chiesa  _et al._ [2005] S. Chiesa, D. M. Ceperley, and S. Zhang, Phys. Rev. Lett. 94, 036404 (2005).
  * Badinski and Needs [2008] A. Badinski and R. J. Needs, Phys. Rev. B 78, 035134 (2008).
  * Assaraf  _et al._ [2011] R. Assaraf, M. Caffarel, and A. C. Kollias, Phys. Rev. Lett. 106, 150601 (2011).
  * Moroni  _et al._ [2014] S. Moroni, S. Saccani, and C. Filippi, J. Chem. Theory Comput. 10, 4823 (2014).
  * Van Rhijn  _et al._ [2021] J. Van Rhijn, C. Filippi, S. De Palo, and S. Moroni, J. Chem. Theory Comput. 18, 118 (2021).
  * Assaraf and Caffarel [2003] R. Assaraf and M. Caffarel, J. Chem. Phys. 119, 10536 (2003).
  * Umrigar [1989] C. J. Umrigar, Int. J. Quantum Chem. 36, 217 (1989).
  * Attaccalite and Sorella [2008] C. Attaccalite and S. Sorella, Phys. Rev. Lett. 100, 114501 (2008).
  * Sorella and Capriotti [2010] S. Sorella and L. Capriotti, J. Chem. Phys. 133, 234111 (2010).
  * Filippi  _et al._ [2016] C. Filippi, R. Assaraf, and S. Moroni, J. Chem. Phys. 144, 194105 (2016).
  * Nakano  _et al._ [2022] K. Nakano, A. Raghav, and S. Sorella, J. Chem. Phys. 156, 034101 (2022).
  * Tiihonen  _et al._ [2021] J. Tiihonen, R. C. Clay III, and J. T. Krogel, J. Chem. Phys. 154, 204111 (2021).
  * Casula and Sorella [2003] M. Casula and S. Sorella, J. Chem. Phys. 119, 6500 (2003).
  * Sorella [2005] S. Sorella, Phys. Rev. B 71, 241103 (2005).
  * Umrigar  _et al._ [2007] C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett. 98, 110201 (2007).
  * Note [1] The AGP ansatz has been generalized also to spin-polarized systems [38].
  * Becca and Sorella [2017] F. Becca and S. Sorella,  _Quantum Monte Carlo approaches for correlated systems_ (Cambridge University Press, 2017).
  * Note [2] $\lambda_{i,j}=\sum@\slimits@_{k}c_{i,k}^{*}c_{j,k}$ in complex cases.
  * Goscinski [1982] O. Goscinski, Int. J. Quantum Chem. 22, 591 (1982).
  * Note [3] Only its real part is taken in complex cases [42].
  * Pathak and Wagner [2020] S. Pathak and L. K. Wagner, AIP Adv. 10, 085213 (2020).
  * Umrigar and Filippi [2005] C. J. Umrigar and C. Filippi, Phys. Rev. Lett. 94, 150201 (2005).
  * Bennett  _et al._ [2017] M. C. Bennett, C. A. Melton, A. Annaberdiyev, G. Wang, L. Shulenburger, and L. Mitas, J. Chem. Phys. 147, 224106 (2017).
  * Bennett  _et al._ [2018] M. C. Bennett, G. Wang, A. Annaberdiyev, C. A. Melton, L. Shulenburger, and L. Mitas, J. Chem. Phys. 149, 104108 (2018).
  * Annaberdiyev  _et al._ [2018] A. Annaberdiyev, G. Wang, C. A. Melton, M. Chandler Bennett, L. Shulenburger, and L. Mitas, J. Chem. Phys. 149, 134108 (2018).
  * Wang  _et al._ [2019] G. Wang, A. Annaberdiyev, C. A. Melton, M. C. Bennett, L. Shulenburger, and L. Mitas, J. Chem. Phys. 151, 144110 (2019).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, Wiley Interdiscip. Rev. Comput. Mol. Sci. 8, e1340 (2018).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, _et al._ , J. Chem. Phys. 153, 024109 (2020).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * Nakano  _et al._ [2020] K. Nakano, C. Attaccalite, M. Barborini, L. Capriotti, M. Casula, E. Coccia, M. Dagrada, C. Genovese, Y. Luo, G. Mazzola, A. Zen, and S. Sorella, J. Chem. Phys. 152, 204121 (2020).
  * Nakano  _et al._ [2023] K. Nakano, O. Kohulák, A. Raghav, M. Casula, and S. Sorella, J. Chem. Phys. 159, 224801 (2023).
  * Posenitskiy  _et al._ [2023] E. Posenitskiy, V. G. Chilkuri, A. Ammar, M. Hapka, K. Pernal, R. Shinde, E. J. Landinez Borda, C. Filippi, K. Nakano, O. Kohulák, S. Sorella, P. de Oliveira Castro, W. Jalby, P. L. Ríos, A. Alavi, and A. Scemama, J. Chem. Phys. 158, 174801 (2023).
  * Huber [2013] K.-P. Huber,  _Molecular spectra and molecular structure: IV. Constants of diatomic molecules_ (Springer Science & Business Media, 2013).
  * Nakano  _et al._ [2021] K. Nakano, T. Morresi, M. Casula, R. Maezono, and S. Sorella, Phys. Rev. B 103, L121110 (2021).
  * Vinet  _et al._ [1987] P. Vinet, J. R. Smith, J. Ferrante, and J. H. Rose, Phys. Rev. B 35, 1945 (1987).
  * Momma and Izumi [2011] K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).


# Switching of the information backflow between a helical spin system and non-
Markovian bath

M. Kaczor  I. Tralle  P.Jakubczyk  S. Stagraczyński  L. Chotorlishvili
<EMAIL_ADDRESS>Institute of Physics, College of Natural
Sciences, University of Rzeszów, Pigonia 1, 35-310 Rzeszów, Poland  Department
of Physics and Medical Engineering, Rzeszów University of Technology, 35-959
Rzeszów, Poland  Institute of Spintronics and Quantum Information, Faculty of
Physics, Adam Mickiewicz University, 61-614 Poznań, Poland

###### Abstract

The dissipative dynamics of the spin chain coupled to the non-Markovian
magnonic reservoir was studied. The chirality of the chain is formed due to
the magnetoelectric coupling. We explored the sign of the trace distance
derivative and found the alternating positive/negative periods in system’s
time evolution. The negative sign is associated with the flow of information
from the system to the bath and decrease in states distinguishability, while
the positive sign is related to the flow of the information in the opposite
direction and increase in distinguishability. We found the distinct effect of
the applied electric and magnetic fields. While the Dzyaloshinskii-Moriya
interaction and external electric field lead to reshuffling of the periods,
the applied magnetic field leads to the swift positive-negative transitions.
Thus, in the helical quantum rings coupled to the non-Markovian magnonic
baths, it is possible to control the directions of information flow through
the external fields.

###### keywords:

Open Quantum Systems , Helical Spin Chains , Non-Markovianity , Quantum
Information

††journal: Annals of Physics

##  1 Introduction

The quantum object in question can be formally split into two parts. Such
procedure allows for describing the smaller piece termed as ”system”
$\hat{H_{0}}$ through the finite set of canonical quantum operators. At the
same time, we postulate the statistical features on the major part termed as a
”bath” $\hat{H_{b}}$. Usually, the expectation values of the bath operators
are described by Bose distribution function $n_{b}=\langle\langle
b^{+}_{k}b_{k}\rangle\rangle$. This concept lies at the heart of the open
quantum systems theory [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] and in the present
work we follow this standard scheme.

The paradigmatic example of an open quantum system is a two-level system
coupled to the phononic (photonic) bath, i.e. the spin-boson model. The
coupling to the bath induces decoherence and dissipation effects in the
system. The stochastic thermostat comprises a set of quasiparticle excitations
and, in many cases, can possess different physical features. For example, in
the case of spin-phonon interaction, the phononic thermostat relaxes swiftly
compared to the spin system. Therefore, the phononic bath recovers instantly
from the spin system’s impact under the continuous flow of information from
the spin system.

The different relaxation time scales between the system and bath permit to
exploit the Markov approximation. However, in several cases, the reservoir has
the same relaxation time scale. For example, the NV centers, are surrounded by
a nuclear spin bath [12, 13, 14, 15]. In this case, the stochastic thermostat
is comprised not of phononic modes, but magnonic modes of the surrounding
nuclear spins, which have relaxation rates similar to the system itself. Due
to the absence of significant difference between system and environment time
scales, Markovian approximation fails. The backflow of information from the
environment towards the system leads to the non-Markovian features [16, 17].

The non-Markovian system-environment information exchange was studied in the
pioneering work [18], within the framework of the stochastic non-Markovian
quantum jump method. Markovian processes reduce the distinguishability between
two states, while non-Markovian processes have the opposite effect. Therefore
non-Markovianity can be quantified through the trace distance between two
density matrices[19, 20, 21].

In the present work, we aim to explore non-Markovian processes in the helical
quantum systems. Multiferroics (MF) constitute a class of materials
simultaneously possessing ferroelectric and magnetic properties[22, 23, 24,
25, 26, 27, 28, 29]. The single-phase MF couples to the external electric
field through the magnetoelectric (ME) term. The mathematical structure of the
ME term is formally identical to the Dzyaloshinskii-Moriya (DM) interaction.
Therefore MF inherits features of helical systems [30]. Our central interest
concerns experimentally feasible helical 1D chains that can be manufactured
and embed on a substrate. Such a chain can be used to transfer quantum
information[26]. Through the exchange interaction, the helical 1D spin chain
interacts with the 2D spin lattices of the substrate, which is the essence of
the non-Markovian spin bath. For the spins of the non-Markovian bath, we
exploit the Holstein-Primakoff transformation and map the spin operators of
the reservoir into the bosonic magnon operators. Spins of the helical ring
coupled to these magnon modes describe the excitation spectrum of the bath.
Before proceeding further, we set the Hamiltonian of the model in question.

Figure 1:  The schematics of the system. The model in question mimics the
atomic Fe chains on the $(5\times 1)-Ir(001)$ surface [26]. In what follows,
we consider the finite length of the chain (up to 50 spins) playing the
helical spin system’s role. The 2D spin substrate contains many spins and
embodies the concept of a non-Markovian spin bath.

Let us consider the Hamiltonian of the 1D chiral ring coupled to the non-
Markovian spin bath. The total Hamiltonian of the system comprises the
Hamiltonian of the chain $\hat{H}_{0}$, the Hamiltonian of the bath
$H_{b}=\sum\limits_{k=k_{1}}^{k_{max}}\omega_{k}b_{k}^{+}b_{k}$, where $k$
innumerates different magnon modes and $\hat{V}$ is the system-bath
interaction term. The total Hamiltonian of the problem reads:

$\displaystyle\hat{H}=\hat{H}_{0}+\hat{H}_{b}+\hat{V}.$ (1)

The Hamiltonian of the system Fig.(1) has the form:

$\displaystyle\hat{H}_{0}$ $\displaystyle=$
$\displaystyle-B\sum\limits_{n=1}^{N}\hat{\sigma}_{n}^{z}+J_{1}\sum\limits_{n}\hat{\sigma}_{n}\hat{\sigma}_{n+1}+J_{2}\sum\limits_{n}\hat{\sigma}_{n}\hat{\sigma}_{n+2}+$
(2)
$\displaystyle+D\sum\limits_{n}\left(\hat{\sigma}_{n}\times\hat{\sigma}_{n+1}\right)_{z}.$

The first two terms in Eq. (2) describe the exchange interaction between
nearest and next-nearest spins, the second term represents the DM interaction,
which is characterized by the ME coupling constant $D=c_{ME}E$.
$\textbf{E}=(0,E,0)$ is the applied external electric field. By changing the
direction of the electric field, we can invert the sign of DM constant. The
Hamiltonian $\hat{H}_{0}$ commutes with the total $z$ projection of the spin
$\hat{\sigma}^{z}=\sum\limits_{n=1}^{N}\hat{\sigma}^{z}_{n}$. Therefore
Hamiltonian has a block diagonal structure. Each sector corresponding to the
particular $\hat{\sigma}^{z}=N,N-1,...,-N$ can be diagonalized separately. In
what follows, we study the low-temperature case and limit the discussion to
the single excitation sector. For further convenience, we diagonalize
Hamiltonian of the ring and present it in spectral decomposition form:
$\hat{H_{0}}=\sum\limits_{n=1}^{N}E_{n}|n\rangle\langle n|$, where
$E_{n},~{}|n\rangle$ are the eigenvalues and eigenfunctions respectively. The
interaction term we assume in the following form:
$V=\sum\limits_{n}^{N}\left(|n\rangle\langle g|B+|g\rangle\langle
n|B^{\dagger}\right)$, where $B=\sum\limits_{k}g_{k}b_{k}$. We explore the
impact of the electric and magnetic fields on the backflow of information from
the non-Markovian bath to the system. The paper is organized as follows: in
Section 2, we describe the analytical solution of Schrödinger equation with
total Hamiltonian of the problem. In Section 3, we present results of
numerical calculation describing the information flow and discuss the
possibility of its control through the electric and magnetic field and
conclude the work.

##  2 Analytical solution

In what follows, we are interested in the low excitation properties of the
system. We define the zero excitation state of the system $\hat{H}_{0}$ as
$|g\rangle=|0_{1}0_{2}...0_{N}\rangle$, with energy $E_{g}=-BN$, for
$J_{1}=-J<0,~{}J_{2}=-J_{1}=J$ and following the standard procedure tackle
non-Markovian problem in the single excitation basis [31]:
$|n\rangle=\frac{1}{\sqrt{N}}\sum\limits_{l=1}^{N}\exp(-inl)|1_{l}\rangle$,
where $|1_{l}\rangle=|0...1_{l}...0_{N}\rangle$, and $E_{n}=J_{1}\cos(2\pi
n/N)+J_{2}\cos(4\pi n/N)+D\sin(2\pi n/N)-B(N-1)$. The general state of the
coupled system we assumed as follows:

$|\phi(t)\rangle=\left(c_{0}|g\rangle+\sum\limits_{n}^{N}c_{n}(t)|n\rangle_{S}\right)|0\rangle_{R}+\sum_{k}f_{k}(t)|g\rangle_{S}|1_{k}\rangle_{R},$
(3)

where
$|1_{k}\rangle_{R}=|1_{{k_{1}}}\rangle|1_{{k_{2}}}\rangle...|1_{k_{max}}\rangle$
mean different magnon modes $\omega_{k}$. In what follows we consider up to
$k_{max}=200$ magnons in the reservoir. The Schrödinger equation reads:

$\displaystyle i\frac{d|\phi(t)\rangle}{dt}=\hat{H}(t)|\phi(t)\rangle,$ (4)

or in the explicit form:

$\displaystyle
i\frac{dc_{n}(t)}{dt}=\omega_{n}c_{n}(t)+\sum_{k}g_{k}f_{k}(t),$
$\displaystyle
i\frac{df_{k}(t)}{dt}=\omega_{k}f_{k}(t)+\sum_{m}g_{k}c_{m}(t),$
$\displaystyle n=(1,...N),~{}~{}n\neq g.$ (5)

Here we introduce the eigenfrequency of the system
$\omega_{n}=(E_{n}-BN)/\hbar$. After performing gauge transformations
$c_{n}\equiv e^{-i\omega_{n}t}c_{n}$, $f_{k}\equiv e^{-i\omega_{k}t}f_{k}$ we
deduce:

$\displaystyle
i\frac{dc_{n}(t)}{dt}=\sum_{k}g_{k}e^{i(\omega_{n}-\omega_{k})t}f_{k}(t),$
$\displaystyle
i\frac{df_{k}(t)}{dt}=\sum_{m}g_{k}e^{-i(\omega_{m}-\omega_{k})t}c_{m}(t),$
$\displaystyle n=(1,...N),~{}~{}n\neq g.$ (6)

The exact set of equations Eq.(2) can be solved numerically, and we use such a
solution in the describtion of a information backflow in a system (see Sec.
3). While the adapted formalism implies a single magnon exchange between the
system and bath [32, 12, 33], the number of magnons in the bath itself is
infinite. Here we present approximated analytical solution, following recipes
from [32, 12, 34, 33] to solve Eq.(2). The formal solution of the second
equation reads:

$\displaystyle
f_{k}(t)=-i\sum_{m}g_{k}\int_{0}^{t}e^{-i(\omega_{m}-\omega_{k})t_{1}}c_{m}(t_{1})dt_{1}.$
(7)

We insert this expression into the first equation in Eq.(2) to obtain the
closed equation

$\displaystyle\frac{dc_{n}(t)}{dt}=-\sum_{m}\int_{0}^{t}c_{m}(t_{1})\sum_{k}g_{k}^{2}e^{-i(\omega_{g}+\omega_{m}-\omega_{k})(t-t_{1})}dt_{1}.$
(8)

Here for the sake of simplicity we took into account that
$|\omega_{m}-\omega_{n}|<|\omega_{m}|,~{}|\omega_{n}|$ and high frequency
osculating terms average to zero. For further simplification we postulate
spectral properties of the non-Markovian bath and rewrite Eq.(8) in the form

$\displaystyle i\frac{dc_{n}(t)}{dt}=-\sum_{m\neq
n}\int\limits_{0}^{t}dt_{1}f_{m}(t-t_{1})c_{m}(t_{1}),$ $\displaystyle
f_{m}(t-t_{1})=\int\limits d\omega
J(\omega)\exp\left[-i(\omega_{m}+\omega)(t-t_{1})\right],$ $\displaystyle
J(\omega)=\frac{1}{2\pi}\frac{\gamma_{0}\lambda^{2}}{(\omega-\omega_{c})^{2}+\lambda^{2}}.$
(9)

Here $\omega_{c}$ is the center frequency of the bath, $\lambda$ defines the
spectral width of the reservoir, $\gamma_{0}\sim 1/\tau$ is related to the
relaxation rate. To solve Eq.(8) we utilize convolution property of Laplace
transform:

$\displaystyle pc_{n}(p)-c_{n}(0)=\sum_{m\neq n}c_{m}(p)f_{m}(p),$
$\displaystyle
c_{n}(t)=\frac{1}{i2\pi}\oint\limits_{\gamma}\exp(pt)c_{n}(p)dp,$ (10)

where contour integration $\gamma$ includes all poles. The explicit solution
is involved and presented in the appendix.

Taking into account the solution we construct the density matrix of the
system. The reduced density matrix reads:

$\displaystyle\hat{\varrho}_{R}^{s}=\left(1-\sum\limits_{n}^{N}|c_{n}(t)|^{2}\right)|g\rangle\langle
g|+\sum_{n=1}^{N}c_{0}c_{n}^{*}|g\rangle\langle n|+$
$\displaystyle\sum\limits_{n,m}^{N}c_{n}(t)c_{m}^{*}(t)|n\rangle\langle
m|+\sum_{n=1}^{N}c_{n}c_{0}^{*}|n\rangle\langle g|.$ (11)

For more details, refer to the appendix.

##  3 Control of the backflow of information

The evolution of the open quantum system can be either Markovian or non-
Markovian. Our interest here concerns the case when entire evolution consists
of an alternation of Markovian and non-Markovian periods. During the Markovian
period, information from the system flows to the bath, and in the non-
Markovian case, information flows back to the system. We aim to control the
process and switch forward and backflows of information. We associate the
direction of the information flow with the sign of the derivative of the trace
distance. The negative sign is associated with the flow of information from
the system to the bath. The positive sign is related to the flow of
information in the opposite direction. The negative sign of the derivative of
trace distance leads to the reduction of distinguishability between the
states, while the positive sign enhances it.

We define the time-derivative of the trace distance as follows:

$\displaystyle\mathcal{\hat{R}}\left(t,\hat{\rho}_{R,1,2}^{s}(0)\right)=\frac{d}{dt}\mathcal{D}\left(\hat{\rho}_{R,1}^{s}(t),\hat{\rho}_{R,2}^{s}(t)\right).$
(12)

where the trace distance is given by

$\displaystyle\mathcal{D}\left(\hat{\rho}_{R,1}^{s}(t),\hat{\rho}_{R,2}^{s}(t)\right)=\frac{1}{2}Tr\left(|\hat{\rho}_{R,1}^{s}(t)-\hat{\rho}_{R,2}^{s}(t)|\right).$
(13)

Two density matrices $\hat{\rho}_{R,1}^{s}(t)$ and $\hat{\rho}_{R,2}^{s}(t)$
differ in the initial conditions.

Figure 2:  Time dependence of the derivative of trace distance $R$. Negative
parts $R<0$ correspond to the flow of information from the system to the bath,
and the positive parts $R>0$ correspond to flow of information in the reverse
direction. One can see that the magnetic field $B$ decreases the modulation
amplitude and leads to the fast switching between positive and negative
regimes. The values of exchange constants $J=-1$, $J_{2}=1$, the value of the
DM constant $D=0.5$. The number of the spins $N=50$, number of the randomly
distributed bath modes $k_{max}=200$. The time scale of the problem is the
picosecond order. Figure 3:  Time dependence of the derivative of trace
distance $R$ without magnetic field $B=0$. Negative parts $R<0$ correspond to
the flow of information from the system to the bath, and the positive parts
$R>0$ correspond to flow of information in the reverse direction. The values
of exchange constants $J=-1$, $J_{2}=1$. The number of the spins $N=50$,
number of the randomly distributed bath modes $k_{max}=200$. The time scale of
the problem is the picosecond order.

The negative sign in Eq.(12) means that process is Markovian and the positive
sign signals the non-Markovianity of the process. We note that the problem of
three spins admits the exact analytical solution presented in the appendix.
For the extended system, we carried out numeric calculations that involved 50
spins and 200 modes of the bath. The results of numeric calculations are shown
in Fig.(2) and Fig.(3). As we see from Fig.(2), applied external magnetic
field leads to the fast transitions between positive and negative periods of
evolution. The applied external electric field, DM interaction and
magnetoelectric coupling leads to the reshuffling of the positive and negative
periods - see Fig.(3). In particular $D=1$ and $D=0$ cases are shifted by
$\pi$ in time oscillations. The temporary backflow of information may occur in
many physical systems, and it does not always imply strong arguments about
thermalization. However, negative sign of the derivative of trace distance
decreases the distance between states and reduces their distinguishability,
while the positive sign increases this feature. Through the external electric
field one can control the distinguishability of the states. The effect we
noticed is related to the deterministic parameter, such as a magnetic field.
The effect we found characterizes the low-temperature properties of helical
systems with broken inversion symmetry.

##  4 Conclusions

In the present project, we study dissipative dynamics of the helical quantum
ring, coupled to the non-Markovian magnonic reservoir. The chiral spin order
in the system is formed due to the magnetoelectric coupling and can be
controlled through the applied electric field. We explored the sign of the
derivative of trace distance and found that the system’s time evolution
consists of the alternating positive-negative periods. The negative (positive)
sign is associated with the flow (backflow) of information from the system
(bath) to the bath (system). The negative sign of the derivative of trace
distance leads to the reduction of distinguishability between the states,
while the positive sign enhances it. Our main findings are as follows: applied
external magnetic field leads to the swift positive vs. negative transitions,
see Fig. (2), while the DM interaction and external electric field lead to the
reshuffling of the positive and negative periods, see Fig.(3).

##  5 Acknowledgment

This work is supported by the Grant No. FR-19-4049 from Shota Rustaveli
National Science Foundation of Georgia. S.S. acknowledges support from the
Norwegian Financial Mechanism under the Polish-Norwegian Research Project NCN
GRIEG “2Dtronics,” Project No. 2019/34/H/ST3/00515

##  Appendix A Lindblad equation

From Eq.(2) we deduce the explicit form of the normalization condition:

$\displaystyle\frac{d}{dt}\left\\{\sum\limits_{n}|c_{n}(t)|^{2}+\sum\limits_{k}|f_{k}(t)|^{2}\right\\}=0,$
$\displaystyle\sum\limits_{n}|c_{n}(t)|^{2}+\sum\limits_{k}|f_{k}(t)|^{2}=1-c_{0}^{2}.$
(14)

Integral in Eq.(2) in the explicit form reads:

$\displaystyle
f\left(t-t_{1}\right)=\frac{\gamma_{0}\lambda^{2}}{2\pi}\exp\left[-i\left(\omega_{m}+\omega_{c}\right)\left(t-t_{1}\right)\right]\int\limits_{-\infty}^{\infty}\frac{\exp\left[-i(t-t_{1})\Omega\right]}{\Omega^{2}+\lambda^{2}}d\Omega,$
$\displaystyle\Omega=\omega-\omega_{c}.$ (15)

After contour integration, from Eq.(A) we deduce the following expression:

$\displaystyle f(t-t_{1})$ $\displaystyle=$
$\displaystyle-\frac{\gamma_{0}\lambda}{2}\left\\{\exp\left[-(i(\omega_{m}+\omega_{c})+\lambda)(t-t_{1})\right]\theta(t-t_{1})+\right.$
(16)
$\displaystyle\left.+\exp\left[-(i(\omega_{m}+\omega_{c})+\lambda)(t_{1}-t)\right]\theta(t_{1}-t)\right\\},$

where $\theta(t-t_{1})$ is the Heaviside step function. For solving Eq.(2) we
exploit the Laplace transform. The Laplace transform of the function (Eq.(16))
reads:

$\displaystyle
f_{m}(p)=-\frac{\gamma_{0}\lambda}{2}\frac{1}{p+\lambda+i(\omega_{m}+\omega_{c})}.$
(17)

We present analytic solution Eq.(2) in the case of three spin ring. Solutions
for the larger systems are cumbersome:

$\displaystyle
c_{1}(p)=\displaystyle{\frac{\begin{bmatrix}c_{1}(0)&-f_{2}(p)&-f_{3}(p)\\\
c_{2}(0)&p&-f_{3}(p)\\\ c_{3}(0)&-f_{2}(p)&p\\\
\end{bmatrix}}{D(p)}},~{}~{}c_{2}(p)=\dfrac{\begin{bmatrix}p&c_{1}(0)&-f_{3}(p)\\\
-f_{1}(p)&c_{2}(0)&-f_{3}(p)\\\ -f_{1}(p)&c_{3}(0)&p\\\
\end{bmatrix}}{D(p)},~{}~{}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}c_{3}(p)=\dfrac{\begin{bmatrix}p&-f_{2}(p)&c_{1}(0)\\\
-f_{1}(p)&p&c_{2}(0)\\\ -f_{1}(p)&-f_{2}(p)&c_{3}(0)\\\ \end{bmatrix}}{D(p)}$
(18)

where:

$D(p)=\begin{bmatrix}p&-f_{2}(p)&-f_{3}(p)\\\ -f_{1}(p)&p&-f_{3}(p)\\\
-f_{1}(p)&-f_{2}(p)&p\\\ \end{bmatrix}.$

For inverse Laplace transform we need to find the roots of the equation

$\displaystyle
D(p)=p^{3}-p(f_{1}(p)f_{2}(p)+f_{1}(p)f_{3}(p)+f_{2}(p)f_{3}(p))-$
$\displaystyle 2f_{1}(p)f_{2}(p)f_{3}(p)=0.$ (19)

For the sake of simplicity we assume that $|c_{1}(0)|^{2}=1-|c_{0}|^{2}$ and
$c_{2}(0)=c_{3}(0)=0$. The expression for $D(p)$ can be rewritten in the other
form:

$\displaystyle
D(p)=\frac{ABCp^{3}-\frac{\gamma_{0}^{2}\lambda^{2}}{4}(A+B+C)p+\frac{\gamma_{0}^{3}\lambda^{3}}{4}}{ABC},$
(20)

where for brevity we introduce the following notations:

$\displaystyle A=p+\lambda+i(\omega_{1}+\omega_{c}),$ $\displaystyle
B=p+\lambda+i(\omega_{2}+\omega_{c}),$ $\displaystyle C\equiv
p+\lambda+i(\omega_{3}+\omega_{c}).$ (21)

One can notice that the numerator in the expression above is the sextic
polynomial. Depending on the values of parameters, all roots can be different
or multiple. At first, we assume that roots are different and write down the
following formulae for the coefficients $c_{i},i=1,2,3$:

$\displaystyle
c_{1}(p)=\frac{\alpha(p,\omega_{1},\omega_{2},\omega_{3})}{(p-p_{1})\times(p-p_{2})\times...\times(p-p_{6})},$
(22) $\displaystyle
c_{2}(p)=\frac{\beta(p,\omega_{1},\omega_{2},\omega_{3})}{(p-p_{1})\times(p-p_{2})\times...\times(p-p_{6})},$
$\displaystyle
c_{3}(p)=\frac{\delta(p,\omega_{1},\omega_{2},\omega_{3})}{(p-p_{1})\times(p-p_{2})\times...\times(p-p_{6})},$

where the functions $\alpha,\beta,\delta$ introduced for convenience, take the
form:

$\displaystyle\alpha(p,\omega_{1},\omega_{2},\omega_{3})=c_{1}(0)\left(ABCp^{2}-\left(\frac{\gamma_{0}\lambda}{2}\right)^{2}-A\right),$
$\displaystyle\beta(p,\omega_{1},\omega_{2},\omega_{3})=c_{1}(0)\left(\frac{\gamma_{0}\lambda}{2}\left(B\frac{\gamma_{0}\lambda}{2}-BCp\right)\right),~{}~{}~{}~{}$
$\displaystyle\delta(p,\omega_{1},\omega_{2},\omega_{3})=c_{1}(0)\left(\frac{\gamma_{0}\lambda}{2}\left(C\frac{\gamma_{0}\lambda}{2}-BCp\right)\right).~{}~{}~{}~{}$
(23)

After calculating the corresponding residua at the simple poles of the
integrand (see Eq. 10) and some algebra, we deduce the explicit fromulae for
the coefficients $c_{i},i=1,2,3$:

$\displaystyle c_{1}(t)$ $\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{6}\prod\limits_{m=1,m\neq
n}^{5}\left(p_{n}-p_{m}\right)^{-1}\times\sqrt{1-c_{0}^{2}}\left((p_{n}+\lambda+i(\omega_{1}+\omega_{c}))\right.$
$\displaystyle\left.(p_{n}+\lambda+i(\omega_{2}+\omega_{c}))(p_{n}+\lambda+i(\omega_{3}+\omega_{c}))p_{n}^{2}-\frac{\gamma_{0}^{2}\lambda^{2}}{4}\right)\exp\left(p_{n}t\right),$
$\displaystyle c_{2}(t)$ $\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{6}\prod\limits_{m=1,m\neq
n}^{5}\left(p_{n}-p_{m}\right)^{-1}\times\sqrt{1-c_{0}^{2}}\left(\frac{\gamma_{0}^{2}\lambda^{2}}{4}(p_{n}+\lambda+i(\omega_{2}+\omega_{c}))\right.$
$\displaystyle\left.-\frac{\gamma_{0}\lambda}{2}p_{n}(p_{n}+\lambda+i(\omega_{2}+\omega_{c}))(p_{n}+\lambda+i(\omega_{3}+\omega_{c}))\right)\exp\left(p_{n}t\right),$
$\displaystyle c_{3}(t)$ $\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{6}\prod\limits_{m=1,m\neq
n}^{5}\left(p_{n}-p_{m}\right)^{-1}\times\sqrt{1-c_{0}^{2}}\left(\frac{\gamma_{0}^{2}\lambda^{2}}{4}(p_{n}+\lambda+i(\omega_{3}+\omega_{c}))\right.$
$\displaystyle\left.-\frac{\gamma_{0}\lambda}{2}p_{n}(p_{n}+\lambda+i(\omega_{2}+\omega_{c}))(p_{n}+\lambda+i(\omega_{3}+\omega_{c}))\right)\exp\left(p_{n}t\right).$

Here $p_{n}$ are the roots of the expression:

$\displaystyle
D(p)=\frac{ABCp^{3}-\frac{\gamma_{0}^{2}\lambda^{2}}{4}(A+B+C)p+\frac{\gamma_{0}^{3}\lambda^{3}}{4}}{ABC}.$
(24)

Now we consider the case of two roots with the multiplicity two and four. Then
doing as previously, calculating the residua at the poles of the order of four
and two, we get:

$\displaystyle c_{1}(t)=4\sqrt{1-c_{0}^{2}}\left(\frac{ABC\cdot
p_{1}^{2}-\left(\frac{\gamma_{0}\lambda}{2}\right)^{2}A}{(p_{2}-p_{1})^{5}}e^{p_{1}t}+\frac{ABC\cdot
p_{2}^{2}-\left(\frac{\gamma_{0}\lambda}{2}\right)^{2}A}{(p_{1}-p_{2})^{5}}e^{p_{2}t}\right),$
$\displaystyle
c_{2}(t)=4\sqrt{1-c_{0}^{2}}\left(\frac{\frac{\gamma_{0}\lambda}{2}\left(B\frac{\gamma_{0}\lambda}{2}-BCp_{1}\right)}{(p_{2}-p_{1})^{5}}e^{p_{1}t}+\frac{\frac{\gamma_{0}\lambda}{2}\left(B\frac{\gamma_{0}\lambda}{2}-BCp_{2}\right)}{(p_{1}-p_{2})^{5}}e^{p_{2}t}\right),$
$\displaystyle
c_{3}(t)=4\sqrt{1-c_{0}^{2}}\left(\frac{\frac{\gamma_{0}\lambda}{2}\left(C\frac{\gamma_{0}\lambda}{2}-BCp_{1}\right)}{(p_{2}-p_{1})^{5}}e^{p_{1}t}+\frac{\frac{\gamma_{0}\lambda}{2}\left(C\frac{\gamma_{0}\lambda}{2}-BCp_{2}\right)}{(p_{1}-p_{2})^{5}}e^{p_{2}t}\right),$

where $A,B,C$ defined as previous and $p_{1},p_{2}$ are the corresponding
polynomial roots. Finally consider the third possibility when the polynomial
has three roots, all of the degree 2. Then doing in the same way, one can get
the following expression for the coefficient $c_{1}$:

$\displaystyle
c_{1}(t)=\frac{2\alpha(p_{1},\omega_{1},\omega_{2},\omega_{3})e^{p_{1}t}}{(p_{1}-p_{2})^{3}(p_{1}-p_{3})^{2}+(p_{1}-p_{2})^{2}(p_{1}-p_{3})^{3}}+$
$\displaystyle\frac{2\alpha(p_{2},\omega_{1},\omega_{2},\omega_{3})e^{p_{2}t}}{(p_{2}-p_{1})^{3}(p_{2}-p_{3})^{2}+(p_{2}-p_{1})^{2}(p_{1}-p_{3})^{3}}+~{}~{}$
$\displaystyle\frac{2\alpha(p_{3},\omega_{1},\omega_{2},\omega_{3})e^{p_{3}t}}{(p_{3}-p_{1})^{3}(p_{3}-p_{2})^{2}+(p_{3}-p_{1})^{2}(p_{3}-p_{2})^{3}}.~{}~{}~{}~{}~{}~{}$
(25)

The other two coefficients, $c_{2}(t),c_{3}(t)$ can be obtained from the
previous formula substituting $\alpha(p_{i},\omega_{1},\omega_{2},\omega_{3})$
by corresponding values of the functions $\beta(...)$ and $\delta(...)$.

Taking into account normalization condition:

$\displaystyle|c_{0}(0)|^{2}+\sum\limits_{n}|c_{n}(t)|^{2}+\sum\limits_{k}|f_{k}(t)|^{2}=1,$
(26)

we obtain the explicit expression of the reduced density matrix of the system:

$\displaystyle\hat{\varrho}_{R}^{s}=\left(1-\sum\limits_{n}^{N}|c_{n}(t)|^{2}\right)|g\rangle\langle
g|+\sum_{n=1}^{N}c_{0}c_{n}^{*}|g\rangle\langle n|+$
$\displaystyle\sum\limits_{n,m}^{N}c_{n}(t)c_{m}^{*}(t)|n\rangle\langle
m|+\sum_{n=1}^{N}c_{n}c_{0}^{*}|n\rangle\langle g|.~{}~{}~{}~{}$ (27)

In the matrix form $N=3$:

$\displaystyle\hat{\varrho}_{R}^{s}=\begin{bmatrix}(1-\sum_{n=1}^{N}|c_{n}(t)|^{2})&c_{0}(t)c_{1}^{*}(t)&c_{0}(t)c_{2}^{*}(t)&c_{0}(t)c_{3}^{*}(t)\\\
c_{1}(t)c_{0}^{*}(t)&|c_{11}(t)|^{2}&c_{1}(t)c_{2}^{*}(t)&c_{1}(t)c_{3}^{*}(t)\\\
c_{2}(t)c_{0}^{*}(t)&c_{2}(t)c_{1}^{*}(t)&|c_{22}(t)|^{2}&c_{2}(t)c_{3}^{*}(t)\\\
c_{3}(t)c_{0}^{*}(t)&c_{3}(t)c_{1}^{*}(t)&c_{3}(t)c_{2}^{*}(t)&|c_{33}(t)|^{2}\\\
\end{bmatrix},$ (28)

where time dependent coefficients are defined in Eq.(20).

## References

  * Diehl et al. [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. Büchler, P. Zoller,  Quantum states and phases in driven open quantum systems with cold atoms,  Nature Physics 4 (2008) 878–883.
  * Campisi et al. [2009] M. Campisi, P. Talkner, P. Hänggi,  Fluctuation theorem for arbitrary open quantum systems,  Physical review letters 102 (2009) 210401.
  * Zhang et al. [2012] W.-M. Zhang, P.-Y. Lo, H.-N. Xiong, M. W.-Y. Tu, F. Nori,  General non-markovian dynamics of open quantum systems,  Physical review letters 109 (2012) 170402.
  * Deffner and Lutz [2011] S. Deffner, E. Lutz,  Nonequilibrium entropy production for open quantum systems,  Physical review letters 107 (2011) 140404.
  * Sieberer et al. [2016] L. M. Sieberer, M. Buchhold, S. Diehl,  Keldysh field theory for driven open quantum systems,  Reports on Progress in Physics 79 (2016) 096001.
  * Laine et al. [2012] E.-M. Laine, H.-P. Breuer, J. Piilo, C.-F. Li, G.-C. Guo,  Nonlocal memory effects in the dynamics of open quantum systems,  Physical review letters 108 (2012) 210402.
  * Rivas [2020] Á. Rivas,  Strong coupling thermodynamics of open quantum systems,  Physical Review Letters 124 (2020) 160601.
  * Song et al. [2019] F. Song, S. Yao, Z. Wang,  Non-hermitian skin effect and chiral damping in open quantum systems,  Physical review letters 123 (2019) 170401.
  * Xu et al. [2019] X. Xu, J. Thingna, C. Guo, D. Poletti,  Many-body open quantum systems beyond lindblad master equations,  Physical Review A 99 (2019) 012106.
  * Tamascelli et al. [2018] D. Tamascelli, A. Smirne, S. F. Huelga, M. B. Plenio,  Nonperturbative treatment of non-markovian dynamics of open quantum systems,  Physical review letters 120 (2018) 030402.
  * Carrega et al. [2016] M. Carrega, P. Solinas, M. Sassetti, U. Weiss,  Energy exchange in driven open quantum systems at strong coupling,  Physical review letters 116 (2016) 240403.
  * Singh et al. [2020] A. Singh, L. Chotorlishvili, S. Srivastava, I. Tralle, Z. Toklikishvili, J. Berakdar, S. Mishra,  Generation of coherence in an exactly solvable nonlinear nanomechanical system,  Physical Review B 101 (2020) 104311.
  * Rabl et al. [2009] P. Rabl, P. Cappellaro, M. G. Dutt, L. Jiang, J. Maze, M. D. Lukin,  Strong magnetic coupling between an electronic spin qubit and a mechanical resonator,  Physical Review B 79 (2009) 041302.
  * Chotorlishvili et al. [2013] L. Chotorlishvili, D. Sander, A. Sukhov, V. Dugaev, V. Vieira, A. Komnik, J. Berakdar,  Entanglement between nitrogen vacancy spins in diamond controlled by a nanomechanical resonator,  Physical Review B 88 (2013) 085201.
  * Mishra et al. [2014] S. K. Mishra, L. Chotorlishvili, A. Rau, J. Berakdar,  Three-level spin system under decoherence-minimizing driving fields: Application to nitrogen-vacancy spin dynamics,  Physical Review A 90 (2014) 033817.
  * Breuer et al. [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, B. Vacchini,  Colloquium: Non-markovian dynamics in open quantum systems,  Reviews of Modern Physics 88 (2016) 021002.
  * De Vega and Alonso [2017] I. De Vega, D. Alonso,  Dynamics of non-markovian open quantum systems,  Reviews of Modern Physics 89 (2017) 015001.
  * Piilo et al. [2008] J. Piilo, S. Maniscalco, K. Härkönen, K.-A. Suominen,  Non-markovian quantum jumps,  Phys. Rev. Lett. 100 (2008) 180402. URL: https://link.aps.org/doi/10.1103/PhysRevLett.100.180402. doi:10.1103/PhysRevLett.100.180402.
  * Breuer et al. [2009] H.-P. Breuer, E.-M. Laine, J. Piilo,  Measure for the degree of non-markovian behavior of quantum processes in open systems,  Phys. Rev. Lett. 103 (2009) 210401. URL: https://link.aps.org/doi/10.1103/PhysRevLett.103.210401. doi:10.1103/PhysRevLett.103.210401.
  * Li et al. [2019] C.-F. Li, G.-C. Guo, J. Piilo,  Non-markovian quantum dynamics: What does it mean?,  EPL (Europhysics Letters) 127 (2019) 50001.
  * Li et al. [2020] C.-F. Li, G.-C. Guo, J. Piilo,  Non-markovian quantum dynamics: What is it good for?,  EPL (Europhysics Letters) 128 (2020) 30001.
  * Azimi et al. [2014] M. Azimi, L. Chotorlishvili, S. Mishra, S. Greschner, T. Vekua, J. Berakdar,  Helical multiferroics for electric field controlled quantum information processing,  Physical Review B 89 (2014) 024424.
  * Cheong and Mostovoy [2007] S.-W. Cheong, M. Mostovoy,  Multiferroics: a magnetic twist for ferroelectricity,  Nature materials 6 (2007) 13–20.
  * Seki et al. [2008] S. Seki, Y. Yamasaki, M. Soda, M. Matsuura, K. Hirota, Y. Tokura,  Correlation between spin helicity and an electric polarization vector in quantum-spin chain magnet licu 2 o 2,  Physical review letters 100 (2008) 127201.
  * Katsura et al. [2005] H. Katsura, N. Nagaosa, A. V. Balatsky,  Spin current and magnetoelectric effect in noncollinear magnets,  Physical review letters 95 (2005) 057205.
  * Menzel et al. [2012] M. Menzel, Y. Mokrousov, R. Wieser, J. E. Bickel, E. Vedmedenko, S. Blügel, S. Heinze, K. von Bergmann, A. Kubetzka, R. Wiesendanger,  Information transfer by vector spin chirality in finite magnetic chains,  Physical review letters 108 (2012) 197204.
  * Azimi et al. [2016] M. Azimi, M. Sekania, S. Mishra, L. Chotorlishvili, Z. Toklikishvili, J. Berakdar,  Pulse and quench induced dynamical phase transition in a chiral multiferroic spin chain,  Physical Review B 94 (2016) 064423.
  * Wang et al. [2018] X.-g. Wang, L. Chotorlishvili, G.-h. Guo, J. Berakdar,  Electric field controlled spin waveguide phase shifter in yig,  Journal of Applied Physics 124 (2018) 073903.
  * Khomeriki et al. [2015] R. Khomeriki, L. Chotorlishvili, B. Malomed, J. Berakdar,  Creation and amplification of electromagnon solitons by electric field in nanostructured multiferroics,  Physical Review B 91 (2015) 041408.
  * Wang et al. [2019] X.-g. Wang, L. Chotorlishvili, G.-h. Guo, C.-L. Jia, J. Berakdar,  Thermally assisted skyrmion drag in a nonuniform electric field,  Physical Review B 99 (2019) 064426.
  * Lorenzo et al. [2013] S. Lorenzo, F. Plastina, M. Paternostro,  Tuning non-markovianity by spin-dynamics control,  Physical Review A 87 (2013) 022317.
  * Li et al. [2010] J.-G. Li, J. Zou, B. Shao,  Non-markovianity of the damped jaynes-cummings model with detuning,  Physical Review A 81 (2010) 062124.
  * Man et al. [2014] Z.-X. Man, N. B. An, Y.-J. Xia,  Non-markovianity of a two-level system transversally coupled to multiple bosonic reservoirs,  Physical Review A 90 (2014) 062104.
  * Breuer et al. [2002] H.-P. Breuer, F. Petruccione, et al., The theory of open quantum systems, Oxford University Press, 2002.