# Self-bound crystals of antiparallel dipolar mixtures

Maria Arazo  Departament de Física Quàntica i Astrofísica, Universitat de
Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain  Institut de Ciències
del Cosmos, Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona,
Spain  Albert Gallemí  Institut für Theoretische Physik, Leibniz Universität
Hannover, 30167 Hannover, Germany  Montserrat Guilleumas  Departament de
Física Quàntica i Astrofísica, Universitat de Barcelona, Martí i Franquès 1,
08028 Barcelona, Spain  Institut de Ciències del Cosmos, Universitat de
Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain  Ricardo Mayol
Departament de Física Quàntica i Astrofísica, Universitat de Barcelona, Martí
i Franquès 1, 08028 Barcelona, Spain  Institut de Ciències del Cosmos,
Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain  Luis
Santos  Institut für Theoretische Physik, Leibniz Universität Hannover, 30167
Hannover, Germany

###### Abstract

Recent experiments have created supersolids of dipolar quantum droplets. The
resulting crystals lack, however, a genuine cohesive energy and are maintained
by the presence of an external confinement, bearing a resemblance to the case
of ion Coulomb crystals. We show that a mixture of two antiparallel dipolar
condensates allows for the creation of potentially large, self-bound crystals
which, resembling ionic crystals in solid-state physics, are maintained by the
mutual dipolar attraction between the components, with no need of transversal
confinement. This opens intriguing novel possibilities, including three-
dimensionally self-bound droplet-ring structures, stripe/labyrinthic patterns,
and self-bound crystals of droplets surrounded by an interstitial superfluid,
resembling the case of superfluid Helium in porous media.

††preprint: APS/123-QED

Solid-state crystals are held together by the interplay between different
forms of attractive and repulsive interactions between their constituents
Kittel (2004). This interplay results in a finite cohesive or binding energy,
defined as the energy that must be added to the crystal to separate its
components infinitely apart. In the presence of an external confinement,
crystals may form even in the absence of genuine cohesion. A prominent example
is provided by trapped ions, which form crystals due to the combination of
repulsive Coulomb interactions and external confinement Thompson (2015). There
is, however, no cohesive energy, and ion Coulomb crystals unravel in the
absence of the trap.

This feature is shared by recently created crystals of quantum droplets in
dipolar Bose-Einstein condensates Böttcher  _et al._ (2020); Chomaz  _et al._
(2022). Self-bound droplets, elongated along the dipole direction, result from
the quasi-cancellation of contact and dipolar interactions, and the
stabilizing effect of quantum fluctuations Petrov (2015); Kadau  _et al._
(2016); Chomaz  _et al._ (2016); Schmitt  _et al._ (2016). In the presence of
confinement along the dipole direction, energy is minimized by the creation of
multiple droplets, which, in the presence of an external confinement
perpendicular to the dipole orientation (transversal trap), arrange forming a
crystal Wenzel  _et al._ (2017) that may present supersolid properties Tanzi
_et al._ (2019a); Böttcher  _et al._ (2019); Chomaz  _et al._ (2019); Natale
_et al._ (2019); Tanzi  _et al._ (2019b); Guo  _et al._ (2019); Norcia  _et
al._ (2021); Bland  _et al._ (2022a). Similar to the case of ions in Coulomb
crystals, droplets repel each other. There is hence no genuine cohesive energy
of the droplet crystal (or of any other possible density pattern Zhang  _et
al._ (2019, 2021); Hertkorn  _et al._ (2021); Poli  _et al._ (2021)). The
transversal trap is crucial to keep it bound.

Recent experiments have created a mixture of two dipolar components Trautmann
_et al._ (2018); Durastante  _et al._ (2020); Politi  _et al._ (2022). These
mixtures are expected to present rich physics due to the competition between
intra- and inter-component contact and dipolar interactions, including
immiscible droplets Smith  _et al._ (2021); Bisset  _et al._ (2021), doping-
induced droplet nucleation Politi  _et al._ (2022); Scheiermann  _et al._
(2023), two-fluid supersolidity Scheiermann  _et al._ (2023), and the
formation of alternating-domain supersolids Li  _et al._ (2022); Bland  _et
al._ (2022b); Kirkby  _et al._ (2023). Interestingly, the dipoles of the two
components may be antiparallel, and hence the inter- and intra-component
interactions may have opposite sign Bland  _et al._ (2022b) (Fig. 1(a)).

Figure 1: (a) In an ADM, intra-component interactions are attractive (yellow
arrows) when the particles are head-to-tail, and repulsive when they are side-
by-side, whereas the opposite is true for the inter-component ones. (b)
Dipolar interaction $V^{1\to 2}_{\rm dd}({\bf
r})=\frac{2\mu_{0}\mu_{1}^{2}}{3\pi l_{z}^{2}}U^{1\to 2}_{\rm dd}$ that
component 1 exerts on component 2, as a function of $z$ and
$\rho=\sqrt{x^{2}+y^{2}}$. For simplicity, we have assumed a Gaussian droplet
$e^{-z^{2}/l_{z}^{2}}e^{-\rho^{2}/2l^{2}}$. The dashed line indicates the
half-width-at-half-maximum of the droplet. The inter-component dipolar
interaction results in an energy minimum on the $xy$ plane at a given radius
$\rho_{0}$ well outside the droplet.

In this letter, we investigate crystal formation in an antiparallel dipolar
mixture (ADM). As for a parallel one Smith  _et al._ (2021); Bisset  _et al._
(2021), in the absence of any confinement, an ADM may form an immiscible
three-dimensionally self-bound mixture, although with a markedly different
topology in which one of the components may eventually form a ring around a
droplet of the other. The presence of confinement along the dipole direction
results in crystal formation. In stark contrast to both single-component
dipolar condensates and parallel binary mixtures, in an ADM the crystal has a
genuine cohesive energy, remaining self-bound in the absence of a transversal
trap due to the mutual attraction between the components. This resembles the
case of ionic crystals in solid-state physics, where ions of opposite charge
arrange in an intertwined crystalline structure bound by their mutual
electrostatic interaction Kittel (2004). However, the resulting self-bound ADM
is not given by two intertwined droplet arrays. Symmetric ADMs with similar
intra-component interaction strengths form self-bound stripe/labyrinthic
density patterns. In contrast, in sufficiently asymmetric ADMs, one of the
components forms an incoherent droplet crystal with an approximate triangular
structure, whereas the second one remains superfluid and fills the lattice
interstitials, resembling to some extent superfluid Helium in porous media
Reppy (1992).

#### Model.–

We consider a bosonic ADM, with dipoles oriented, respectively, along and
antiparallel to the $z$ axis. The components may belong to the same species or
to two different ones. In order to illustrate the possible physics, we
consider a dysprosium mixture, with magnetic dipoles $\mu_{1}=10\,\mu_{B}$ and
$\mu_{2}=-10\,\mu_{B}$, with $\mu_{B}$ the Bohr magneton. Short-range
interactions are characterized by the intra- and inter-component scattering
lengths: $a_{11}$, $a_{22}$, and $a_{12}$. The physics of the mixture is well
described by the extended Gross-Pitaevskii equation Smith  _et al._ (2021);
Bisset  _et al._ (2021):

$\displaystyle i\hbar\dot{\Psi}_{\sigma}({\bf r},t)$ $\displaystyle=$
$\displaystyle\left[\frac{-\hbar^{2}\nabla^{2}}{2m}+V_{\rm trap}({\bf
r})+\sum_{\sigma^{\prime}}g_{\sigma\sigma^{\prime}}|\Psi_{\sigma^{\prime}}({\bf
r},t)|^{2}\right\delimiter 00$ (1) $\displaystyle+$
$\displaystyle\sum_{\sigma^{\prime}}\int
d^{3}r^{\prime}V^{\sigma\sigma^{\prime}}_{\rm dd}({\bf r}-{\bf
r}^{\prime})|\Psi_{\sigma^{\prime}}({\bf r}^{\prime},t)|^{2}$ $\displaystyle+$
$\displaystyle\mu_{{\rm LHY},\sigma}[n_{1,2}({\bf
r},t)]\Big{]}\Psi_{\sigma}({\bf r},t)\,,$

where $\Psi_{\sigma}({\bf r},t)$ is the condensate wavefunction of component
$\sigma=1,2$, $n_{\sigma}=|\Psi_{\sigma}|^{2}$, and
$g_{\sigma\sigma^{\prime}}=4\pi\hbar^{2}a_{\sigma\sigma^{\prime}}/m$, with $m$
the mass of the bosons. The atoms are confined, if at all, only along the $z$
axis by a potential $V_{\rm trap}({\bf r})=\frac{1}{2}m\omega_{z}^{2}z^{2}$.
The dipole-dipole interaction is given by the potential
$V^{\sigma\sigma^{\prime}}_{\rm dd}({\bf
r})=\frac{\mu_{0}\mu_{\sigma}\mu_{\sigma^{\prime}}}{4\pi
r^{3}}\left(1-3\cos^{2}\theta\right)$, with $\theta$ the angle sustained by
the $z$ axis and ${\bf r}$. The effect of quantum fluctuations is provided by
the Lee-Huang-Yang (LHY) term $\mu_{{\rm LHY},\sigma}[n_{1,2}({\bf
r},t)]=\delta E_{\rm LHY}/\delta n_{\sigma}$, where

$E_{\rm
LHY}=\frac{8}{15\sqrt{2\pi}}\Big{(}\frac{m}{4\pi\hbar^{2}}\Big{)}^{3/2}\int
d\theta_{k}\sin\theta_{k}\sum_{\lambda=\pm}V_{\lambda}(\theta_{k})^{5/2}\,$
(2)

is the LHY energy correction, with

$V_{\pm}(\theta_{k})=\sum_{\sigma=1,2}\eta_{\sigma\sigma}n_{\sigma}\pm\sqrt{(\eta_{11}n_{1}-\eta_{22}n_{2})^{2}+4\eta_{12}^{2}n_{1}n_{2}}$
(3)

and
$\eta_{\sigma\sigma^{\prime}}=g_{\sigma\sigma^{\prime}}+g^{d}_{\sigma\sigma^{\prime}}(3\cos^{2}\theta_{k}-1)$,
being
$g^{d}_{\sigma\sigma^{\prime}}=\mu_{0}\mu_{\sigma}\mu_{\sigma^{\prime}}/3$ and
$\theta_{k}$ the angle sustained by ${\bf k}$ with the $z$ axis. Since the
dipole moments of the components are antiparallel, the inter-component dipolar
potential is repulsive (attractive) when the components are placed head-with-
tail (side-by-side) (see Fig. 1(a)). As a result, the dipolar interaction
strongly favors immiscibility, and a very large and negative $a_{12}$ is
needed to drive the system miscible. In the following, we consider
$a_{12}=150\,a_{0}$, but the actual value is irrelevant as long as the inter-
component overlapping remains negligible.

#### Three-dimensionally self-bound ADM.–

We first consider the case of fully unconfined mixtures ($\omega_{z}=0$). As
for parallel dipolar mixtures Smith  _et al._ (2021); Bisset  _et al._ (2021),
an immiscible ADM may present a three-dimensionally self-bound solution, but
of a markedly different nature. This is best understood in the impurity limit
($N_{1}\gg N_{2}$). Let us assume that component 1 forms a self-bound droplet
with density $n_{1}({\bf r})$. The droplet exerts a potential $V^{1\to 2}_{\rm
dd}({\bf r})=\int d^{3}r^{\prime}V^{12}_{\rm dd}({\bf r}-{\bf
r}^{\prime})n_{1}({\bf r}^{\prime})$ on component 2, which, as seen in Fig.
1(b), is characterized by a marked minimum at a given radius $\rho_{0}$, well
outside the droplet. Particles in component 2 are trapped in this mexican-hat
potential.

Figure 2: Three-dimensionally self-bound ADMs. Ground-state configuration as a
function of the total atom number $N$ and of $a_{22}$, for $a_{11}=50\,a_{0}$
and $N_{1,2}=N/2$. Whereas component $1$ always forms a single elongated
droplet, component 2 may acquire different topologies, which we characterize
using the separation $\Delta r_{\mathrm{CM}}$ between the center of masses of
the two components (color code). The different topologies are illustrated in
the insets, where we depict the column density (integrated over $z$) of the
components, with red (blue) indicating component 1 (2).

In a more balanced mixture, the argument remains valid, but component $2$ also
induces a similar potential $V^{2\to 1}_{\rm dd}({\bf r})$ on component $1$.
Hence the two components confine each other mutually on the $xy$ plane,
resulting in self-bound ADMs, as illustrated in Fig. 2 for $N_{1,2}=N/2$,
$a_{11}=50\,a_{0}$, and different values of $a_{22}$ and $N$. For asymmetric
intra-component interactions $a_{11}<a_{22}$, component $1$ remains a compact
droplet, whereas the second component accommodates on the ring potential
around the droplet. For low enough $a_{22}$, the energy is minimized by the
formation of a single droplet in component 2, which for growing $N$ and
$a_{22}$ spreads around the mexican-hat minimum until eventually forming a
ring-like configuration. For intermediate $a_{22}$ values, there is a second
possible topology with two droplets of component 2 placed at opposite sides of
the annular potential.

#### Self-bound droplet crystals.–

When $\omega_{z}=0$, increasing the particle number $N$ results in more
elongated solutions along the $z$ direction. As for single-component (scalar)
dipolar condensates Böttcher  _et al._ (2020), this elongation is frustrated
in the presence of a trap along $z$ ($\omega_{z}>0$). In scalar condensates,
this frustration results in the formation of multiple droplets. Although the
droplets repel each other, the presence of a transversal trap on the $xy$
plane allows for the creation of 2D droplet crystals Norcia  _et al._ (2021);
Bland  _et al._ (2022a). These crystals have however no intrinsic cohesion,
and hence unravel in the absence of the $xy$ confinement.

Figure 3: Self-bound droplet crystals. (a) Phase diagram as a function of the
atom number $N$ and the trap frequency $f_{z}=\omega_{z}/2\pi$, for
$a_{11}=50\,a_{0}$ and $a_{22}=70\,a_{0}$. Colors correspond to configurations
with a different number of droplets $N_{D}$ in component 1. Figures (b–g) show
the column magnetization (integrated along $z$) of the lowest-energy solution
for selected cases, indicated with the corresponding symbol in Fig. (a). Red
(blue) regions are populated by component 1 (2). Figure 4: (a,b) Single-shot
realizations of the column magnetization for $a_{11}=50\,a_{0}$,
$a_{22}=70\,a_{0}$, $N_{1,2}=5\times 10^{4}$ and $\omega_{z}/2\pi=1300$ Hz.
Red (blue) regions are populated by component 1 (2). (c,d) Corresponding
momentum distribution for the second component, $\tilde{n}_{2}(k_{x},k_{y})$
in the $k_{z}=0$ plane, in the cases of Figs. (a) and (b), respectively. (e)
Momentum distribution $\tilde{n}_{2}(k_{x},k_{y})$ averaged over $10$
different realizations.

Remarkably, this is not the case in an ADM, as illustrated in Fig. 3 for a
balanced mixture $N_{1}=N_{2}$ and asymmetric intra-component interactions,
$a_{11}=50\,a_{0}$ and $a_{22}=70\,a_{0}$. For a low-enough $\omega_{z}$, the
three-dimensional solution (with a single droplet in component 1) remains
valid (Fig. 3(b)). For an $N$-dependent critical $\omega_{z}$ the droplet
splits into two. Each one of them exerts a mexican-hat potential on the second
component, which gets trapped in the combined energy minimum. At the same
time, crucially, the second component glues the two droplets together, forming
a self-bound ADM (Fig. 3(c)). As shown in Fig. 3(a), and illustrated for
particular cases in Figs. 3(d–g), further increasing $\omega_{z}$ results in a
growing number of droplets of component $1$ surrounded by a bath of component
2. In a scalar condensate, each droplet requires a minimal atom number to
remain self-bound (otherwise kinetic energy unbinds it), drastically limiting
the total number of droplets. In contrast, in an ADM, droplets remain confined
by the inter-component interaction, allowing for droplets with a much smaller
number of atoms Li  _et al._ (2022); Bland  _et al._ (2022b). As a result,
increasing $\omega_{z}$ results in 2D crystals with much more droplets
compared to scalar condensates with the same total number of atoms.

We should emphasize that our results, based on imaginary-time evolution of Eq.
(1) with random initial conditions, reveal many possible solutions with very
similar energy, which differ in the exact number and arrangement of the
droplets (see SM ). We hence expect a significant experimental shot-to-shot
variability, similar to that recently observed in experiments on 2D
supersolids Norcia  _et al._ (2021).

#### Interstitial superfluid.–

Due to the lack of overlapping, the droplets are mutually incoherent. In
contrast, the component filling the crystal interstitials forms a superfluid
that resembles, to some extent, the case of Helium in a porous medium
(although, in contrast to that scenario, droplets of component 1 do not form a
rigid structure). The approximately triangular crystalline structure of the
droplets is inherited as well by the interstitial component 2, which builds
hence a peculiar form of supersolid. The coherence and spatial density
modulation of component 2 may be revealed in time-of-flight measurements.
Figure 4 shows the momentum distribution $\tilde{n}_{2}(k_{x},k_{y})$ in the
$k_{z}=0$ plane. The approximate triangular structure (Fig. 4(a)) results in
an hexagonal pattern in the $\tilde{n}_{2}$ distribution (Fig. 4(c)), although
the above-mentioned variability of the exact droplet arrangement may result in
a significant shot-dependent distortion (see Figs. 4(b,d)). Note as well that,
due to the lack of any confinement on the $xy$ plane, the patterns
spontaneously break the polar symmetry and hence experience a random rotation
from shot to shot. In any case, as expected from the theory of roton
immiscibility Wilson  _et al._ (2012); Bland  _et al._ (2022b), the inter-
droplet distance $R$ is fixed by the oscillator length
$a_{z}=\sqrt{\hbar/m\omega_{z}}$. For the case of Fig. 3, $R\simeq 3\,a_{z}$
for all values of $N$ and $\omega_{z}$. This periodicity becomes evident from
the average of the momentum distribution over many realizations, which shows a
marked ring at $1/R$ (see Fig. 4(e)).

Figure 5: (a) Phase diagram for $N_{1,2}=5\times 10^{4}$ and
$\omega_{z}/2\pi=1200$ Hz. Two different self-bound solutions are found: a
droplet crystal (D-CRYSTAL), illustrated in Fig. (b) for $a_{11}=50\,a_{0}$
and $a_{22}=70\,a_{0}$, and a stripe/labyrinthic (STR-LAB) phase, illustrated
in Fig. (c) for $a_{11}=55\,a_{0}$ and $a_{22}=60\,a_{0}$. The case of a well-
defined stripe phase is illustrated in Fig. (d), which has been evaluated for
$\omega_{z}/2\pi=180$ Hz, $a_{11}=80\,a_{0}$, and $a_{22}=80\,a_{0}$. In Figs.
(b–d) we depict the column magnetization. Red (blue) regions are populated by
component 1 (2).

#### Crystal sublimation.–

For a fixed total number of particles, the cohesive energy decreases when the
droplet number grows, since lowering the density reduces the inter-component
dipolar attraction. Eventually, at a critical frequency
$\omega_{z}^{\mathrm{cr}}$, the crystal unbinds, and both components
evaporate. The critical frequency ($\omega_{z}^{\mathrm{cr}}/2\pi\simeq 1400$
Hz for the case on Fig. 3) is approximately determined as that for which the
energy per particle reaches $\hbar\omega_{z}/2$, corresponding to an
infinitely spread solution on the $xy$ plane. Interestingly, when $\omega_{z}$
approaches $\omega_{z}^{\mathrm{cr}}$, mutual attraction may still be enough
to maintain a stable crystal, but insufficient to bind the whole interstitial
superfluid, which hence partially evaporates (see SM  for a more detailed
discussion).

#### Self-bound stripe/labyrinthic patterns.–

Up to this point, we have considered a mixture with markedly asymmetric intra-
component interactions. Interestingly, when $a_{11}\simeq a_{22}$, the mixture
arranges in a different form of self-bound pattern (note that $a_{11}=a_{22}$
if we consider a mixture of two maximally stretched magnetic states of the
same atomic species). This is illustrated by the phase diagram of Fig. 5 (a),
obtained for $\omega_{z}/2\pi=1200$ Hz and $N_{1,2}=5\times 10^{4}$. For
sufficiently large $|a_{11}-a_{22}|$, we obtain the above-mentioned droplet
crystal (Fig. 5(b)), which, as mentioned above, presents partial evaporation
of the interstitial component in the vicinity of the unbinding threshold. In
contrast, when $a_{11}\simeq a_{22}$ the mixture arranges in a labyrinthic
phase, with a large shot-to-shot variability, formed by stripes with different
orientations (Fig. 5(c)). For lower trap frequencies, the ground-state
configuration is given by a well-defined stripe crystal (Fig. 5(d)). Note that
in the labyrinthic/stripe phase both components form mutually incoherent
domains.

#### Summary and Outlook.–

Antiparallel dipolar mixtures allow for the formation of crystals with a
genuine cohesive energy that remain self-bound in the absence of a transversal
trap. The mutual confinement stems from the attractive inter-component
interactions, and results in incoherent stripe/labyrinthic crystals in
mixtures with symmetric intra-component interactions, and self-bound droplet
crystals in asymmetric mixtures. The latter are particularly interesting,
since while one component forms an approximately triangular array of
incoherent droplets, the other component builds a superfluid in the
interstitials, forming a peculiar form of supersolid that may be readily
probed using time-of-flight measurements. Although we have considered the
particular example of a dysprosium mixture, our results generally apply to
other antiparallel magnetic or electric dipolar mixtures, including those of
polar molecules.

The possibility of creating self-bound dipolar crystals opens intriguing
perspectives for future studies, including the character of lattice
excitations, which may remain self-bound or result in phonon evaporation
(resembling droplet evaporation in non-dipolar mixtures Petrov (2015)), the
probing (e.g. by vortex formation) of the superfluidity of the interstitial
component, as well as in general the exploration of the dynamics of self-bound
crystals.

###### Acknowledgements.

This work has been supported by Grant No. PID2020-114626GB-I00 (Ministerio de
Ciencia e Innovación), by the European Union Regional Development Fund within
the ERDF Operational Program of Catalunya (project QUASICAT/QuantumCat), by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy–EXC-2123 QuantumFrontiers–390837967, and FOR
2247. M. A. is supported by FPI Grant PRE2018-084091.

## References

  * Kittel (2004) C. Kittel,  _Introduction to Solid State Physics_, 8th ed. (Wiley, 2004).
  * Thompson (2015) R. C. Thompson, Contemporary Physics 56, 63 (2015), https://doi.org/10.1080/00107514.2014.989715 .
  * Böttcher  _et al._ (2020) F. Böttcher, J.-N. Schmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, Reports on Progress in Physics 84, 012403 (2020).
  * Chomaz  _et al._ (2022) L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Reports on Progress in Physics 86, 026401 (2022).
  * Petrov (2015) D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * Kadau  _et al._ (2016) H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier, I. Ferrier-Barbut, and T. Pfau, Nature 530, 194 (2016).
  * Chomaz  _et al._ (2016) L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Phys. Rev. X 6, 041039 (2016).
  * Schmitt  _et al._ (2016) M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, Nature 539, 259 (2016).
  * Wenzel  _et al._ (2017) M. Wenzel, F. Böttcher, T. Langen, I. Ferrier-Barbut, and T. Pfau, Phys. Rev. A 96, 053630 (2017).
  * Tanzi  _et al._ (2019a) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019a).
  * Böttcher  _et al._ (2019) F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051 (2019).
  * Chomaz  _et al._ (2019) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Phys. Rev. X 9, 021012 (2019).
  * Natale  _et al._ (2019) G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 123, 050402 (2019).
  * Tanzi  _et al._ (2019b) L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Nature 574, 382 (2019b).
  * Guo  _et al._ (2019) M. Guo, F. Böttcher, J. Hertkorn, J.-N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, Nature 574, 386 (2019).
  * Norcia  _et al._ (2021) M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. N. Bisset, L. Santos, and F. Ferlaino, Nature 596, 357 (2021).
  * Bland  _et al._ (2022a) T. Bland, E. Poli, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, L. Santos, and R. N. Bisset, Phys. Rev. Lett. 128, 195302 (2022a).
  * Zhang  _et al._ (2019) Y.-C. Zhang, F. Maucher, and T. Pohl, Phys. Rev. Lett. 123, 015301 (2019).
  * Zhang  _et al._ (2021) Y.-C. Zhang, T. Pohl, and F. Maucher, Phys. Rev. A 104, 013310 (2021).
  * Hertkorn  _et al._ (2021) J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. Res. 3, 033125 (2021).
  * Poli  _et al._ (2021) E. Poli, T. Bland, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, R. N. Bisset, and L. Santos, Phys. Rev. A 104, 063307 (2021).
  * Trautmann  _et al._ (2018) A. Trautmann, P. Ilzhöfer, G. Durastante, C. Politi, M. Sohmen, M. J. Mark, and F. Ferlaino, Phys. Rev. Lett. 121, 213601 (2018).
  * Durastante  _et al._ (2020) G. Durastante, C. Politi, M. Sohmen, P. Ilzhöfer, M. J. Mark, M. A. Norcia, and F. Ferlaino, Phys. Rev. A 102, 033330 (2020).
  * Politi  _et al._ (2022) C. Politi, A. Trautmann, P. Ilzhöfer, G. Durastante, M. J. Mark, M. Modugno, and F. Ferlaino, Phys. Rev. A 105, 023304 (2022).
  * Smith  _et al._ (2021) J. C. Smith, D. Baillie, and P. B. Blakie, Phys. Rev. Lett. 126, 025302 (2021).
  * Bisset  _et al._ (2021) R. N. Bisset, L. A. P. n. Ardila, and L. Santos, Phys. Rev. Lett. 126, 025301 (2021).
  * Scheiermann  _et al._ (2023) D. Scheiermann, L. A. P. n. Ardila, T. Bland, R. N. Bisset, and L. Santos, Phys. Rev. A 107, L021302 (2023).
  * Li  _et al._ (2022) S. Li, U. N. Le, and H. Saito, Phys. Rev. A 105, L061302 (2022).
  * Bland  _et al._ (2022b) T. Bland, E. Poli, L. A. P. n. Ardila, L. Santos, F. Ferlaino, and R. N. Bisset, Phys. Rev. A 106, 053322 (2022b).
  * Kirkby  _et al._ (2023) W. Kirkby, T. Bland, F. Ferlaino, and R. N. Bisset,  (2023), 10.48550/ARXIV.2301.08007.
  * Reppy (1992) J. D. Reppy, Journal of Low Temperature Physics 87, 205 (1992).
  * (32) See Supplemental Material at [URL will be inserted by publisher] for additional details on XXX.
  * Wilson  _et al._ (2012) R. M. Wilson, C. Ticknor, J. L. Bohn, and E. Timmermans, Phys. Rev. A 86, 033606 (2012).

##  Appendix A Supplementary Information

The supplementary information contains additional details on the variability
of the density patterns and the partial evaporation of the mixture mentioned
in the main text.

###  A.1 Shot-to-shot variability

There is a large shot-to-shot variability of the exact number of droplets and
their arrangement in the droplet crystal. We illustrate this point with Fig.
S1, where we show different configurations for the same parameters
$a_{11}=50\,a_{0}$, $a_{22}=70\,a_{0}$, $\omega_{z}/2\pi=1200$ Hz, and
$N_{1,2}=5\times 10^{4}$. The configurations, which have an energy per
particle $E/N\simeq 0.22\,\hbar\omega_{z}$, differ in energy by less than
$1\%$, and have a number of droplets ranging from $23$ to $31$.

Figure S1: Shot-to-shot variability. Different droplet lattice configurations
obtained for the same parameters $a_{11}=50\,a_{0}$, $a_{22}=70\,a_{0}$,
$\omega_{z}/2\pi=1200$ Hz, $N_{1,2}=5\times 10^{4}$. The plots show the column
magnetization, with red (blue) indicating component 1 (2).

Note that the droplet arrangement is only approximately triangular, and the
actual distribution, which is generally non-uniform, may significantly depart
from a triangular lattice. This is reflected in the structure of the
interstitial component, which results in the distorted momentum distributions
depicted in Fig. 4 of the main text. A similar shot-to-shot variability is
observed in the stripe/labyrinthic phase.

###  A.2 Evaporation of the droplet crystal

When $\omega_{z}$ increases, the number of droplets grows and the density
decreases. As a result, the cohesive energy is reduced, and the self-bound
solution eventually unbinds. Figure S2 shows the energy per particle as a
function of the trap frequency for $a_{11}=50\,a_{0}$, $a_{22}=70\,a_{0}$ and
$N_{1,2}=5\times 10^{4}$. The unbinding of the droplet crystal occurs
approximately when the energy per particle $E/N$ reaches $\hbar\omega_{z}/2$,
which, for the case of Fig. S2, occurs at $\omega_{z}^{\mathrm{cr}}/2\pi\simeq
1400$ Hz. Indeed, beyond that value we do not find well defined self-bound
solutions in our simulations.

When $\omega_{z}$ approaches $\omega_{z}^{\mathrm{cr}}$, the crystal remains
bound, but the interstitial component may present partial evaporation. In
order to take this into account we considered absorbing boundary conditions in
our imaginary-time simulations. We fixed a given radius $\rho_{c}$ on the $xy$
plane, such that the crystal is well contained in a circle of radius
$\rho<\rho_{c}$. Particles that reach $\rho>\rho_{c}$ during imaginary-time
evolution are considered as evaporated. We indicate in Fig. S2 the proportion
of the interstitial component that is evaporated, which as expected grows when
approaching $\omega_{z}^{\mathrm{cr}}$. In the unbound regime, the whole
mixture eventually leaves (in imaginary time) the region $\rho<\rho_{c}$.

Figure S2: Energy per particle (blue circles) as a function of the trap
frequency for $a_{11}=50\,a_{0}$, $a_{22}=70\,a_{0}$ and $N_{1,2}=5\times
10^{4}$. The dashed line depicts the energy per particle ($\hbar\omega_{z}/2$)
corresponding to an infinitely spread mixture on the $xy$ plane. The red
squares indicate the proportion of atoms in component 2 which are evaporated
(see text).


# Heat-conserving three-temperature model for ultrafast demagnetization of
$3d$ ferromagnets

M. Pankratova<EMAIL_ADDRESS>Department of Physics and
Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden  I. P.
Miranda  Department of Physics and Astronomy, Uppsala University, Box 516,
SE-75120 Uppsala, Sweden  D. Thonig  School of Science and Technology, Örebro
University, SE-701 82, Örebro, Sweden  Department of Physics and Astronomy,
Uppsala University, Box 516, SE-75120 Uppsala, Sweden  M. Pereiro  Department
of Physics and Astronomy, Uppsala University, Box 516, SE-75120 Uppsala,
Sweden  E. Sjöqvist  Department of Physics and Astronomy, Uppsala University,
Box 516, SE-75120 Uppsala, Sweden  A. Delin  Department of Applied Physics,
School of Engineering Sciences, KTH Royal Institute of Technology, AlbaNova
University Center, SE-10691 Stockholm, Sweden  SeRC (Swedish e-Science
Research Center), KTH Royal Institute of Technology, SE-10044 Stockholm,
Sweden  P. Scheid  Université de Lorraine, LPCT, CNRS, UMR 7019, BP 70239,
54506 Vandoeuvre-lés-Nancy Cedex, France  Université de Lorraine, IJL, CNRS,
UMR 7198, BP 70239, 54000 Nancy Cedex, France  Department of Physics and
Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden  O. Eriksson
Department of Physics and Astronomy, Uppsala University, Box 516, SE-75120
Uppsala, Sweden  A. Bergman  Department of Physics and Astronomy, Uppsala
University, Box 516, SE-75120 Uppsala, Sweden

###### Abstract

We study the ultrafast magnetization dynamics of bcc Fe and fcc Co using the
recently suggested heat-conserving three-temperature model (HC3TM), together
with atomistic spin- and lattice dynamics simulations. It is shown that this
type of Langevin-based simulation is able to reproduce observed trends of the
ultrafast magnetization dynamics of fcc Co and bcc Fe, in agreement with
previous findings for fcc Ni. The simulations are performed by using
parameters that to as large extent as possible are obtained from electronic
structure theory. The one parameter that was not calculated in this way, was
the damping term used for the lattice dynamics simulations, and here a range
of parameters were investigated. It is found that this term has a large
influence on the details of the magnetization dynamics. The dynamics of iron
and cobalt is compared with previous results for nickel and similarities and
differences in the materials’ behavior are analysed following the absorption
of a femtosecond laser pulse. Importantly, for all elements investigated so
far with this model, we obtain a linear relationship between the value of the
maximally demagnetized state and the fluence of the laser pulse, which is in
agreement with experiments.

##  I INTRODUCTION

Ultrafast demagnetization was discovered by Beaurepaire and coauthors in 1996
[1]. They observed demagnetization in a nickel film on picosecond timescales
following the absorption of a femtosecond laser pulse. From the point of view
of applications, ultrafast demagnetization is an important process in all-
optical magnetization switching as well as for new applications in magnetic
data storage and spintronics [2]. In the same pioneering work of Beuarepaire
[1], these experimental observations were interpreted using a three-
temperature model (3TM) [1], which assumes three thermalised reservoirs, in
particular, spin-, lattice, and electron reservoirs, that can exchange energy
through coupling parameters (via electron-phonon, electron-spin, and spin-
lattice coupling). The 3TM is often used to interpret ultrafast magnetization
dynamics processes [1, 3, 4]. Recently, many other models have been proposed
[2, 3, 5, 6] to describe possible mechanisms of ultrafast demagnetization,
such as the importance of spin-dependent transport of laser-excited electrons
[7], the optical inter-site spin transfer effect (OISTR) that was considered
in Ref. [8], and the Elliott–Yafet electron-phonon spin-flip scattering,
studied in Ref. [9]. In [10] a dynamic spin-lattice-electron model was
proposed. Using this model, the authors calculated laser-induced
demagnetization of iron thin film and obtained very accurate agreement with
the experimental observations. One of the important outcomes of the [10] is
also establishing the relation between the dissipative parameters entering the
Langevin equations for the lattice and spin degrees of freedom and the heat
transfer coefficients of 3TM.

In addition, Zahn et al. [11, 12] proposed an energy-conserving model based on
atomistic spin dynamics simulations. Using this model they studied ultrafast
demagnetization of nickel [11], iron and cobalt [12]. However, the model of
Zahn et al. relies on the microscopic electron-lattice heat transfer
coefficient which is hard to estimate. Literature values for the electron-
phonon coupling parameter Gep can vary for nickel up to one order of magnitude
[11, 4]. For iron, reported values of the electron-phonon coupling parameter
differs from $7\cdot 10^{17}$ to $5.48\cdot 10^{18}\,\mathrm{W/m^{3}K}$ [12],
and for cobalt from $6\cdot 10^{17}$ to $4.05\cdot 10^{18}\,\mathrm{W/m^{3}K}$
(please see [12] and references therein).

One of the reasons for this uncertainty of the electron-lattice heat transfer
coefficient is that the electron-phonon coupling is often deduced from
observables, without taking into account the energy cost of demagnetization.
Also, these parameters are often estimated indirectly using, for example, a
two-temperature model [13]. It was demonstrated that $G_{ep}$ extracted from
the data employing the two-temperature model is twice smaller than that
calculated from ab initio theory [13]. These differences in reported $G_{ep}$
values makes the interpretation of the experimental observations very
challenging.

To address this, a previous work proposed a heat-conserving three-temperature
model (HC3TM) [4] for calculation on spin, lattice and electron temperatures
during atomistic spin-dynamics simulations. The benefit of this model is that
it reduces the dependence on heat-transfer parameters that are hard to
estimate, such as electron-phonon, $G_{ep}$, electron-spin, $G_{es}$, and the
spin-lattice coupling, $G_{sl}$. In addition, this model allows for a
description of the nickel demagnetization rate on sub-picosecond timescales.

In the present work, we apply the HC3TM model to study ultrafast
demagnetization of bcc Fe and fcc Co and compare with previously published
data on ultrafast dynamics of fcc Ni [4]. The ultrafast demagnetization of
these three ferromagnets was experimentally studied in Ref. [14]. It was
observed that the amplitude of the demagnetization value increases linearly
with the fluence of the light for Fe, Ni, and Co. A model to explain the
experimental observations was proposed based on the assumption that the linear
dependence of demagnetization on fluence is driven by the increase in
temperature, the electron-phonon coupling, the electron-magnon scattering, and
a reduction of the inter-atomic exchange. In the present work, we aim to study
the impact of fluence on the demagnetization of 3d ferromagnets using the
HC3TM. Moreover, we investigate in more detail the impact of various factors
on ultrafast demagnetization, such as Gilbert damping, lattice damping,
microscopic spin-lattice coupling, and laser fluence.

##  II Atomistic spin-lattice dynamics simulations

Spin and lattice dynamics is obtained from coupled atomistic spin- and lattice
dynamics (SLD) simulations, using Langevin dynamics simulations, as proposed
in Ref. [15]:

$\begin{split}\frac{d\bm{m}_{i}}{dt}=-\frac{\gamma}{(1+\alpha^{2})}\bm{m}_{i}\times(\bm{B}_{i}+\bm{B}_{i}^{fl})\\\
-\frac{\gamma}{(1+\alpha^{2})}\frac{\alpha}{m_{i}}\bm{m}_{i}\times(\bm{m}_{i}\times[\bm{B}_{i}+\bm{B}_{i}^{fl}])\end{split}$
(1) $\frac{d\bm{u}_{k}}{dt}=\bm{v}_{k}$ (2)
$\frac{d\bm{v}_{k}}{dt}=\frac{\bm{F}_{k}}{M_{k}}+\frac{\bm{F}_{k}^{fl}}{M_{k}}-\nu\bm{v}_{k},$
(3)

where $\bm{m_{i}}$ represents an atomic magnetic moment, $m_{i}$ and $\gamma$
are the saturation magnetization and the gyromagnetic ratio correspondingly.
Atomic displacements are denoted by $\bm{u_{k}}$, and velocities are
$\bm{v_{k}}$. We obtain an effective exchange field $\bm{B}_{i}=-\partial
H_{\mathrm{SLD}}/\partial\bm{m}_{i}$ from the spin-lattice Hamiltonian,
$H_{SLD}$. The Hamiltonian used in this work includes magnetic, lattice and
spin-lattice coupling parts following Ref. [15]:

$H_{\mathrm{SLD}}=H_{\mathrm{LL}}+H_{\mathrm{SS}}+H_{\mathrm{SSL}}$ (4)

Concretely:

$H_{\mathrm{LL}}=\frac{1}{2}\sum_{kl}\Phi_{kl}^{\mu\nu}u_{k}^{\mu}u_{l}^{\nu}+\frac{1}{2}\sum_{k}M_{k}\nu_{k}^{\mu}\nu_{k}^{\mu}$
(5)

where the force constant tensor $\Phi_{kl}^{\mu\nu}$ is a rank 2 tensor in
real space, and $M_{k}$ is the mass of atom $k$. The magnetic Hamiltonian is
described by

$H_{\mathrm{SS}}=-\frac{1}{2}\sum_{ij}J_{ij}^{\alpha\beta}(0)m_{i}^{\alpha}m_{j}^{\beta}$
(6)

where $J_{ij}^{\alpha\beta}(0)$ is the exchange tensor at the equilibrium
lattice positions. $\alpha,\beta\in\\{x,y,z\\}$ denote Cartesian components in
spin space, while $\mu,\nu\in\\{x,y,z\\}$ corresponds to Cartesian components
in real space.

It was recently suggested that spin-lattice coupling can be represented in
coupled spin-lattice simulations, by a dependence of the exchange interaction
on atomic displacements $\bm{u_{k}}$, resulting in functions of the type
$J_{ij}^{\alpha\beta}(\bm{u_{k}})$ (see Ref. [15]). As the values of
$\bm{u_{k}}$ are usually small, the spin-lattice term in the Hamiltonian can
be obtained by Taylor expanding the magnetic bi-linear Hamiltonian with
respect to the lattice displacements. This results in a spin-lattice coupling
term that is bi-linear in spin and linear in displacements [15]. Concretely,

$H_{\mathrm{SSL}}=-\frac{1}{2}\sum_{ijk}\Gamma_{ijk}^{\alpha\beta\mu}m_{i}^{\alpha}m_{j}^{\beta}u_{k}^{\mu},$
(7)

where we introduce the coupling constant
$\Gamma_{ijk}^{\alpha\beta\mu}\equiv\Gamma_{ijk}^{\alpha\beta\mu}(u_{k}^{\mu})=\partial
J_{ij}^{\alpha\beta}(u_{k}^{\mu})/\partial u_{k}^{\mu}$.

All parameters for the spin-lattice dynamics simulations were obtained from
ab-initio calculations (see Appendix B for details). This includes exchange
interactions, magnetic moments, inter atomic forces, and spin-lattice
couplings. The force at site $k$ is defined by $\bm{F}_{k}=-\partial
H_{\mathrm{SLD}}/\partial\bm{u}_{k}$. Gilbert and lattice damping constants
are denoted $\alpha$ and $\nu$, respectively. In these types of Langevin
simulations one employs stochastic fields, $\bm{B}_{i}^{fl}$ and
$\bm{F}_{k}^{fl}$, as white noise with properties $\langle
B_{i,\mu}^{fl}(t)B_{j,\nu}^{fl}(t^{\prime})\rangle=2D_{M}\delta_{ij}\delta_{\mu\nu}\delta(t-t^{\prime})$
and $\langle
F_{i,\mu}^{fl}(t)B_{j,\nu}^{fl}(t^{\prime})\rangle=2D_{L}\delta_{kl}\delta_{\mu\nu}\delta(t-t^{\prime})$.
In our simulations, we use $D_{M}=\alpha k_{B}T/\gamma m$, $D_{L}=\nu Mk_{B}T$
, where $T$ and $k_{B}$ are temperature and Boltzmann constant respectively
(for details see e.g. Ref. [16]).

##  III Heat-conserving three-temperature model

The heat-conserving three-temperature model was proposed in Ref. [4] for the
study of ultrafast demagnetisation of fcc Ni. Similar to Beurepaire’s three-
temperature model, it is assumed that lattice, and spin systems are connected
to a finite electronic heat bath. Both spin and lattice subsystems have their
own damping; Gilbert or lattice damping, correspondingly. The model does not
explicitly rely on the heat transfer parameters of the conventional 3TM [1].
The difference between HC3TM and the conventional 3TM is easy to see from the
schematic picture in Fig. 1. Overall, instead of microscopic heat transfer
coefficients, as used in 3TM, the HC3TM is formulated in terms of atomistic
damping coefficients. Thus, the HC3TM calculates the heat transfer between the
three reservoirs (spin, lattice and electron)[1] directly from the
simulations. In fact, all parameters of this theory can be calculated using,
for example, density functional theory, although we note that in the present
work we varied the lattice damping parameter in a range of realistic values,
as opposed to performing an explicit calculation of it. The relation between
the lattice damping and $G_{ep}$ was established in [10] using spin-lattice-
electron three-temperature model. It was shown that the lattice damping is
directly proportional to $G_{ep}$, which was previously also reported by Duffy
et el. in [17, 18]. Another distinction of HC3TM is that (unlike 3TM) the
temperature of the different subsystems is calculated at every time step of
the spin-lattice simulation using the expression:

$\Delta
T_{e}(t)=-\frac{C_{l}(T_{l})}{C_{e}(T_{e})}T_{l}(t)-\frac{C_{s}(T_{s})}{C_{e}(T_{e})}T_{s}(t)+\frac{W(t)}{C_{e}(T_{e})},$
(8)

where $T_{l}$ is calculated from the average kinetic energy of the lattice
vibrations; $\langle E^{kin}_{l}\rangle/k_{B}$. The spin temperature, $T_{s}$
is calculated according to [19], i.e. using
$T_{s}=\frac{\langle\sum_{i}\left|\hat{\bm{m}}_{i}\times\bm{B}_{i}\right|^{2}\rangle}{2k_{B}\langle\sum_{i}\hat{\bm{m}}_{i}\cdot\bm{B}_{i}\rangle}$,
where $\hat{\bm{m}_{i}}$ is the normalized local spin moment. The definition
of an instantaneous measure of the temperatures in an out-of-equilibrium
system is not obvious [20], but is still typically assumed to hold in these
kinds of models. Introducing a time average of the temperatures, for narrow
time windows, does however not change the results significantly. In addition,
Eq. (8) contains the temperature dependent specific heats of the electron-,
lattice- and spin systems, $C_{e}(T)$, $C_{l}(T)$ and $C_{s}(T)$,
respectively. The calculations details of the specific heat capacities are
given in Appendix A. The possibility for an external stimulus, e.g. a laser,
to momentarily increase the temperature of the electronic subsystem is
captured by the source term $W(t)$ of Eq. (8), which is modelled as a Gaussian
function: $W(t)=W_{0}\exp(-(t-t_{0})^{2}/2\sigma^{2})$, where $t_{0}$ is the
center, $\sigma=0.02$ ps is the width of the pulse, and $W_{0}$ is the
magnitude (in Joules), correspondingly. More details concerning the HC3TM and
its comparison with 3TM can be found in Ref. [4].

(a) |
---|---
(b) |
Figure 1:  The main idea behind the HC3TM (a) and 3TM (b) models. The text
along the arrows indicates the parameters responsible for heat transfer in
each model. For the parameters of the 3TM, see Ref. [1].

##  IV Results

###  IV.1 HC3TM for cobalt and iron

In this section we present results of spin-lattice dynamics simulations of fcc
Co and bcc Fe. In both cases, we run simulations for cells with a $60\times
60\times 60$ repetition of the fcc (bcc) unit cell, using periodic boundary
conditions. In addition, we used $N_{t}=1\times 10^{6}$ time steps of
$dt=10^{-16}\,\mathrm{s}$, in combination with other parameters presented in
Table 1.

We start by describing results from ultrafast demagnetization of fcc Co.
First, we present temperature profiles (Fig. 2) for spin, lattice, and
electron subsystems, for a laser fluence of $45\,\mathrm{J/m^{2}}$, together
with the corresponding magnetization dynamics. Other parameters relevant for
these simulations are listed in Table 1. It can be seen from the figure that
similar to experimental observations [21] and theoretical studies [12]
demagnetization happens on subpicosecond timescales, followed by a
remagnetization that has a faster recovery after $\sim$ 1 ps and a slower
remagnetization after that. We obtain values of the position of the
magnetization minima that are close to experimental values, and we also obtain
demagnetization/remagnetization timescales that are consistent with
experimental observations. It is important to note that in Ref. [12], using
the parametrized three-temperature model, features of the magnetization
dynamics similar to the data of Fig. 2 were also observed, however for an
adjusted Gilbert damping value of $0.1$, which is unrealistically large. In
contrast, in the calculations presented here we use a value of damping from
ab-inito calculations (that is close to experiment), that as seen in Table 1
is much smaller than the value used in Ref. [12].

Table 1: List of the parameters used in the simulations.  Parameter | Value
---|---
Iron |
Gilbert damping $\alpha$ | 0.008 [10]
Lattice damping $\nu$ |  $1\times 10^{-13}$ kg/s
Cobalt |
Gilbert damping $\alpha$ | 0.0024 ([22]) or 0.0014 [23]
Lattice damping $\nu$ |  $1\times 10^{-13}$ kg/s
Nickel |
Gilbert damping $\alpha$ | 0.024 ([4])
Lattice damping $\nu$ |  $1\times 10^{-13}$ kg/s
Figure 2:  Spin, lattice and electron temperatures (a) and magnetization
dynamics (b) of fcc Co, obtained with the HC3TM model. The pulse fluence is 45
J/m2.

Similar spin-lattice simulations using the HC3TM, were performed for bcc Fe.
Just like for fcc Co we first present magnetization dynamics after the
absorption of a laser pulse and the corresponding spin, lattice, and
electronic temperatures, see Fig. 3. Using the HC3TM model we obtain realistic
demagnetization/remagnetization times in comparison with experimental studies
and in this case even the maximally reduced magnetization is quite close to
the experimentally observed value [9], for the same pulse fluence ($M/M_{0}$
is equal around $0.9$ in our calculations).

Figure 3:  Spin, lattice and electron temperatures (a) and magnetization
dynamics (b) of bcc Fe, obtained with the HC3TM model, pulse fluence is 15
J/m2.

###  IV.2 Results from coupled spin-lattice simulations

To study the impact of direct spin-lattice coupling on the magnetization
dynamics, given by Eq. (7), we consider explicitly the spin-lattice coupling,
where we obtain the elements of the tensor, $\Gamma_{ijk}$ in the spin-lattice
coupling, from calculations based on density functional theory (e.g. as
described in Ref. [15]). Calculated values of these parameters were used to
study the impact of direct spin-lattice coupling for bcc Fe, and fcc Co. This
explicit, calculated spin-lattice coupling allows for a direct exchange of
heat between the spin and lattice subsystems, outside of the channel provided
by the HC3TM. Note that the HC3TM can be applied whether or not the direct
spin-lattice term in Eq. (7) is included or neglected, and the results that
follow were obtained from the HC3TM with Eq. (7) included explicitly in the
simulations.

Figure 4:  Impact of spin-lattice coupling on magnetization dynamics of bcc
Fe. The pulse fluence is 20 J/m2.

In Fig. 4 we show for bcc Fe the demagnetization profile from coupled spin-
lattice dynamics simulations in combination with the HC3TM. The simulations
are, as noted, done with and without the coupling term in Eq. (7). It can be
seen from Fig. 4 that spin-lattice coupling impacts the magnetization dynamics
very weakly for bcc Fe. If any effect can be observed, it would be a slightly
slower remagnetization rate when the spin-lattice coupling is included. The
difference compared to the scenario without the spin-lattice coupling is
however almost of the same order of magnitude as the numerical noise of the
simulations. We note that a simulation which completely ignores this term
still has ability to transfer heat from the spin system to the lattice, via
Eq. (8). A possible explanation for the low influence of the spin-lattice
couplings can be drawn as follows: even at the maximum spin temperature in the
simulations ($T_{s}\lesssim 600$ K for Fe), the spin system is far away from
$T_{c}$, i.e. with a general strong coupling between spin moments. Thus, the
small average displacements (in the order of $10^{-2}a$) and their
distribution throughout the sample result in the summed effect of Eq. (7)
being only a minor perturbation on the total exchange field.

Also, for bcc Fe the resulting profile is in rather good agreement with
experimental observations. For materials with larger spin-lattice coupling
[24] it can be expected that the importance of Eq. (7) increases, and our
simulations show that for bcc Fe the spin-lattice coupling becomes more
important for higher pulse fluences (data not shown).

While the impact of spin-lattice coupling, via Eq. (7), on the magnetization
dynamics is only marginal for the elemental ferromagnts, the impact of the
dynamical properties of the lattice itself, and in particular, the lattice
damping, is significant. We illustrate this for fcc Co, in Fig. 5. Here it is
seen that a decrease of the lattice damping ($\nu$ in Eq. (3) leads to an
increase of the demagnetization; both the magnetization minimum as well as the
reduced magnetization ($M\over M_{0}$) are influenced by the lattice damping
parameter. This can best be understood from the fact that the coupled spin-
and lattice system allows for heat to be dissipated from both reservoirs, when
the lattice (and spin) damping parameter is larger. This causes the
temperature profile to reach lower values and to equilibrate quicker both for
the spin- and lattice system, as shown in Fig. 5, with a resulting stronger
impact on the magnetization profile for lower values of the lattice damping.
This is detailed clearly when comparing the spin temperatures (Fig. 5) for
various values of lattice damping. Lower lattice damping values lead to higher
spin temperature and therefore, to larger magnetization drop and a longer
remagnetization time. One of the important conclusions from this observation
is that it is essential to take the lattice dynamics properly into account
while studying ultrafast demagnetization. If the lattice is not considered at
all, then one needs to note, while comparing simulations with experimental
data, that this will lead to an overestimation of the magnetization drop and
to longer remagnetization times, for a given laser pulse.

Figure 5:  Spin temperature (a) and magnetization dynamics (b) of fcc Co for
various lattice damping values used in HC3TM ($\nu$ in Eq. (3)), with
temperature-dependent heat capacities (see Ref. [4] for details). The pulse
fluence is 45 J/m2.

###  IV.3 Comparison of magnetization dynamics: bcc Fe, fcc Co, fcc Ni

In this section we discuss magnetization dynamics of bcc Fe and fcc Co. We
compare the results from the two elements to previous results reported for fcc
Ni, which also has been studied using the HC3TM [4]. The results are
summarized in Figs. 6 and 7. We begin with comparing ultrafast magnetization
dynamics of fcc Co, fcc Ni, and bcc Fe, for the same laser pulse fluence,
presented in Fig. 6. It can be seen that in this case, the demagnetization for
fcc Ni is the most prominent among the three systems, while the smallest
effect is found for fcc Co. This difference in demagnetization properties is
consistent with the difference in observed values of $T_{c}$, which is highest
for fcc Co and smallest for fcc Ni.

Figure 6: Ultrafast demagnetization in fcc Ni, fcc Co and bcc Fe for the same
value of the laser pulse fluence $8\,\mathrm{J/m^{2}}$.

The next step in our comparison is to choose values of the pulse fluence for
all three materials resulting in the same (or very similar) magnetization
dynamics. As one can see in Fig. 7, due to very different $T_{c}$ values, we
need a pulse fluence as high as $45\,\mathrm{J/m^{2}}$ for cobalt to exhibit a
demagnetization that is similar to nickel for a fluence of
$4.2\,\mathrm{J/m^{2}}$. It can also be seen from Fig. 7 that in this case the
magnetization dynamics of all three systems are very similar. The remaining
differences are minute, for example, the minima of cobalt’s magnetization
curve appears somewhat later than the one of iron and nickel.

Figure 7:  Ultrafast demagnetization in fcc Ni, fcc Co and bcc Fe for a
different values of laser pulse fluence. The fluence is chosen to obtain the
same demagnetization amplitude.

It is important to take into account that this very similar dynamics is
observed for specific parameter values, such as the Gilbert and lattice
damping. For example, to obtain the curve for cobalt presented in Fig. 7, we
used the damping value from Ref. [22], see Table 1. However, with the damping
value obtained from ab initio calculations, the magnetization dynamics will be
somewhat different. In particular, we then observe much slower remagnetization
of cobalt (as, for example, shown in Fig. 5 obtained with the Gilbert damping
set to 0.0014).

Figure 8:  The dependence of the maximally demagnetized value versus laser
fluence for bcc Fe, fcc Co and fcc Ni. The shaded area demonstrate a range of
reported [11, 4, 12] lattice damping parameters obtained using the relation
between lattice damping and $G_{ep}$ from [10].

In Fig. 8 we plot the maximally demagnetized state (minimum value of the
magnetization curve, $M_{min}/M_{0}$) as function of absorbed laser pulse
fluence for bcc Fe, fcc Co, and fcc Ni for the same value of the lattice
damping (solid lines). The first observation is that one needs the lowest
laser fluence to demagnetize fcc Ni and the strongest for fcc Co. This can, as
mentioned earlier in this paper, be connected to the corresponding $T_{c}$.
Secondly, $M_{min}/M_{0}$ in the simulations depends linearly on laser
fluence, something which was observed experimentally in Ref. [14] and also for
other materials in Refs. [25, 26].

In our work we use Gilbert damping that is the highest for nickel and the
lowest for cobalt, which might add to the observed trend in Fig. 8. In
particular, a higher Gilbert damping, which in the HC3TM corresponds to the
electron-magnon coupling, leads to a stronger heat transfer from the
electronic to the spin subsystem, and therefore faster demagnetization, as in
the case of nickel.

In the literature, the reported Gilbert damping parameters of iron, nickel,
and cobalt deviate between different experiments. For example, reported
$\alpha$ values for iron are in the range of $0.0019-0.0072$ [23]. Nickel does
have a consistently larger Gilbert damping than the other two elements,
stemming from the $1/m_{Ni}$ scaling and the peak in the minority spin channel
at the Fermi level [23]. Varying the Gilbert damping in our HC3TM model
results in curves that deviate slightly from those presented in Fig. 8 (data
not shown) but the overall behaviour still remains the same.

To demonstrate the impact of lattice damping on the trend of the
demagnetization more in detail, we have performed spin-lattice simulations for
all three elements using several different values of the lattice damping. In
addition to the simulations for the same lattice damping value ($\nu=1\times
10^{-13}$ kg/s - solid lines in Fig. 8), we also used lattice damping values
corresponding to the upper and lower limits of the electron-phonon coupling
$G_{ep}$ found in the literature [11, 4, 12]. These intervals for the reported
electron-phonon couplings are shown as the shaded areas in Fig. 8 where it can
be noticed that the trend of largest demagnetization for nickel and smallest
demagnetization for cobalt holds over almost the full interval of the shaded
areas. For all three elements, the lower bounds of the shaded areas in Fig. 8
correspond to the lowest lattice damping while the upper bounds indicate the
result of the largest lattice damping of the considered values. This further
exemplifies the general picture that a stronger electron-lattice coupling
transfers more heat from the electron to the lattice sub-system thus
effectively reducing the amount of heat that gets pumped into the electron
system.

In Ref. [14] Scheid et al. experimentally showed that the maximum
demagnetization amplitude in Fe and Co are surprisingly similar. Various
reasons behind the distribution of demagnetization rates in iron, cobalt,
nickel were discussed, including the impact of electron-phonon coupling, the
electron-magnon scattering rate, and the ultrafast light-induced quenching of
the inter-atomic exchange. In our model, as seen Fig. 8. this similarity could
be explained by a very large lattice damping in Fe relatively to Co. Indeed,
in HC3TM, and by using first principle values of Gilbert, to obtain similar
slopes for iron and cobalt, we need to use a lattice damping for iron 10 times
larger than for cobalt. This is illustrated in Fig. 8 where the dashed curves
indicate a lattice damping of $\nu=5\times 10^{-13}$ kg/s for iron and
$\nu=5\times 10^{-14}$ kg/s for cobalt. However recent experimental [12] and
first principles [27] results suggest that the lattice damping in cobalt is in
fact larger than in iron. This would amplify the demagnetization in iron,
while minimizing the one of cobalt. That conclusion of Scheid et al. [14] is
therefore in line with our findings, because even for the same lattice damping
for all systems (which in our model corresponds to electron-phonon coupling,
see Fig. 1) we clearly obtain the highest demagnetization rate for nickel
followed by iron and then cobalt (see solid lines in Fig. 8). One may further
improve the agreement with the experimental observations by considering
additionally ultrafast light-induced quenching of the interatomic exchange as
suggested in Ref. [14].

##  V Conclusions

Using a heat-conserving three-temperature model we have calculated spin,
lattice, and electron temperatures in simulations of the ultrafast
magnetization dynamics of cobalt and iron. We have studied in detail the
impact of Gilbert and lattice damping and compared results of the
magnetization dynamics for bcc Fe, fcc Co and fcc Ni. It was found that in
studies of ultrafast magnetic phenomena, the lattice dynamics plays a
surprisingly important role, even in cases when the direct spin-lattice
coupling is not significant (i.e., when Eq. 7 is neglected). Simulations for
various laser fluences were performed to study and compare ultrafast
demagnetization of iron, cobalt, and nickel in various regimes. We have
demonstrated that the simulations are consistent with recent observations,
which show a linear trend between the maximally demagnetized state and laser
fluence. We furthermore show that the experimentally found trend that for the
same fluence, fcc Ni demagnetizes the most and fcc Co the least, holds for a
wide range of realistic choices of the electron-phonon coupling strength.

##  VI ACKNOWLEDGEMENT

This work was financially supported by the Knut and Alice Wallenberg
Foundation (grant numbers 2018.0060, 2021.0246, and 2022.0108),
Vetenskapsrådet (grant numbers 2019-03666, 2016-05980, and 2019-05304), the
European Research Council (grant number 854843-FASTCORR), the foundation for
Strategic Research SSF, and Olle Engkvist foundation. Support from STandUP and
eSSENCE is also acknowledged. Computations were enabled by resources provided
by the Swedish National Infrastructure for Computing (SNIC) at NSC, partially
funded by the Swedish Research Council through grant agreement no. 2018-05973.
P.S. acknowledges support from the ANR-20-CE09-0013 UFO.

##  Appendix A Heat capacities of the spin, lattice, and electron subsystems

In our calculations we use temperature-dependent heat capacities for the
lattice, electron, and spin subsystems. The corresponding curves are presented
in Fig. 9.

Figure 9:  Heat capacities for the electronic, spin, and lattice subsystems
used in the calculations for bcc Fe.

The lattice heat capacity is obtained from the Debye model. The spin heat
capacity is calculated using a recently proposed approach from quantum
statistics [28]. The method more accurately describes heat capacities at low
temperatures than what is obtained from Boltzmann statistics. However, its
disadvantage is an overestimation of the heat capacity around $T_{c}$, as can
be seen from Fig. 9. Since in our calculations the spin temperature does not
reach $T_{c}$, this shortcoming of the method does not impact our results.
Moreover, our calculations show that using this approach does not change any
of our main conclusions or results. We can confirm that all results and
conclusions such as linear dependence of demagnetization amplitude on fluence,
etc remain the same regardless of heat capacities used in calculations (please
see Appendix A in [4] for details of heat capacities calculations), including
constant (i.e., temperature independent) heat capacities. The only thing that
is affected by the change of capacities calculation is demagnetization
amplitude. For the electronic subsystem, we assume that the electronic heat
capacity is proportional to the electronic temperature
$C_{e}=\gamma_{e}T_{e}$, and for bcc Fe we use $\gamma_{e}=4.9\cdot
10^{-3}\,\mathrm{J\>mol^{-1}\>K^{-2}}$.

##  Appendix B Details of the ab-initio calculations

Similarly to Ref. [4], the calculation of the interactions to parameterize
Eqs. 1 and 3 were performed in the framework of Density Functional Theory
(DFT), using two different methods: (i) the real-space linear muffin-tin
orbital method in the atomic sphere approximation (RS-LMTO-ASA) [29, 30] and
(ii) the plane-wave pseudopotential-based Quantum ESPRESSO (QE) package [31]
in combination with the PHONOPY [32] software. On one hand, this choice is
motivated, essentially, on the well-known shortcoming of the LMTO-ASA method
on obtaining accurate Hellmann-Feynman forces [33] – and, by consequence,
phonon frequencies. On the other hand, to calculate the spin-lattice couplings
(SLC), a supercell approach is necessary because of the breaking of inversion
symmetry and long-range interactions [4]. In this sense, RS-LMTO-ASA is more
suitable for this calculation because it deals directly with real-space.

For the computation of the force constants using a combination of QE and
PHONOPY, the fcc Ni case used the same definitions as described in Ref. [4].
Analogously, in the case of bcc Fe and fcc Co, in QE the scalar relativistic
schemes based on the projector augmented-wave method (PAW) [34] and the
ultrasoft pseudopotentials (USPP) [35], respectively, are employed. As the
exchange-correlation (XC) functional, we use the generalized gradient
approximation (GGA) with the Perdew-Burke-Ernzerhof parametrization [36].
Here, the choice of a spin-polarized GGA XC term is motivated by the fact that
the calculated lattice constants as well as the phonon frequencies better
reproduce the experimental data in comparison with the local spin density
approximation [37]. The iron and cobalt atoms were described by standard
pseudopotentials from the QE library, with $3s$, $3p$, $4s$, $4p$ and $3d$
valence electrons. A cutoff of 100 (1080) and 90 (500) Ry for the kinetic
energy (charge density) was considered for fcc Co and bcc Fe, in this order.
In turn, a Monkhorst-Pack (MP) [38] grid of $24\times 24\times 24$
$\vec{k}$-points was set for the first Brillouin zone integration. The self-
consistent calculations with a $10^{-10}$ Ry convergence threshold were
carried out using the Marzari-Vanderbilt [39] cold smearing with a spreading
of 0.01 Ry. The force constants and phonon frequencies, computed by PHONOPY,
were based on $6\times 6\times 6$ supercells (216 atoms), where we considered
a $4\times 4\times 4$ $\vec{k}$-points mesh in the self-consistent cycle.

The exchange parameters, Gilbert damping ($\alpha$), density of states and
spin-lattice couplings (SLC) were obtained using the RS-LMTO-ASA method. In
this real-space formalism, the eigenvalue problem is solved with the help of
the recursion method [40], where the recursion chain is ended after $LL$ steps
by making use of the Beer-Pettifor terminator [41]. For bcc Fe we considered
$LL=31$, while for fcc Co a much higher $LL$ value ($LL=51$) is needed to
better describe the density of states and Green’s functions at the Fermi
level. In line with recent experience on finding a model to calculate reliable
$\alpha$’s from first principles [23], the bulk systems consist of a big
cluster in real-space containing $\sim 55000$ (bcc) and $\sim 696000$ (fcc)
atoms, located in the perfect crystal positions. Differently from the QE
calculations (in which the structures were optimized), the lattice parameters
were fixed to the experimental values of $a=2.87\,$Å (bcc Fe) and $a=3.54\,$Å
(fcc Co) [42]. In this setup, the local density approximation (LSDA), with the
von Barth and Hedin parametrization [43], was used. In particular, this choice
of XC was based on the fact that some of quantities that parameterize Eqs. 1-7
were already investigated using LSDA [23, 44, 45, 46]. It is a known trait,
however, that GGA produces quantities with similar quality for the $3d$
magnetic elements [46, 47, 23]. Finally, the spin-orbit coupling (SOC) is
included as a $l\cdot s$ term, computed in each variational step [48, 49].

## References

  * Beaurepaire  _et al._ [1996] E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Phys. Rev. Lett. 76, 4250 (1996).
  * Kirilyuk  _et al._ [2010] A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * Scheid  _et al._ [2022] P. Scheid, Q. Remy, S. Lebègue, G. Malinowski, and S. Mangin, Journal of Magnetism and Magnetic Materials 560, 169596 (2022).
  * Pankratova  _et al._ [2022] M. Pankratova, I. P. Miranda, D. Thonig, M. Pereiro, E. Sjöqvist, A. Delin, O. Eriksson, and A. Bergman, Phys. Rev. B 106, 174407 (2022).
  * Locht  _et al._ [2015] I. L. M. Locht, I. Di Marco, S. Garnerone, A. Delin, and M. Battiato, Phys. Rev. B 92, 064403 (2015).
  * Shim  _et al._ [2020] J.-H. Shim, A. A. Syed, J.-I. Kim, H.-G. Piao, S.-H. Lee, S.-Y. Park, Y. S. Choi, K. M. Lee, H.-J. Kim, J.-R. Jeong, _et al._ , Scientific Reports 10, 6355 (2020).
  * Battiato  _et al._ [2010] M. Battiato, K. Carva, and P. M. Oppeneer, Phys. Rev. Lett. 105, 027203 (2010).
  * Willems  _et al._ [2020] F. Willems, C. von Korff Schmising, C. Strüber, D. Schick, D. W. Engel, J. Dewhurst, P. Elliott, S. Sharma, and S. Eisebitt, Nature communications 11, 1 (2020).
  * Carpene  _et al._ [2008] E. Carpene, E. Mancini, C. Dallera, M. Brenna, E. Puppin, and S. De Silvestri, Phys. Rev. B 78, 174422 (2008).
  * Ma  _et al._ [2012] P.-W. Ma, S. L. Dudarev, and C. H. Woo, Phys. Rev. B 85, 184301 (2012).
  * Zahn  _et al._ [2021] D. Zahn, F. Jakobs, Y. W. Windsor, H. Seiler, T. Vasileiadis, T. A. Butcher, Y. Qi, D. Engel, U. Atxitia, J. Vorberger, and R. Ernstorfer, Phys. Rev. Research 3, 023032 (2021).
  * Zahn  _et al._ [2022] D. Zahn, F. Jakobs, H. Seiler, T. A. Butcher, D. Engel, J. Vorberger, U. Atxitia, Y. W. Windsor, and R. Ernstorfer, Phys. Rev. Research 4, 013104 (2022).
  * Waldecker  _et al._ [2016] L. Waldecker, R. Bertoni, R. Ernstorfer, and J. Vorberger, Phys. Rev. X 6, 021003 (2016).
  * Scheid  _et al._ [2023] P. Scheid, J. Hohlfeld, G. Malinowski, A. Bergman, O. Eriksson, S. Lebègue, and S. Mangin, arXiv preprint arXiv:2301.01055 10.48550/arXiv.2301.01055 (2023).
  * Hellsvik  _et al._ [2019] J. Hellsvik, D. Thonig, K. Modin, D. Iuşan, A. Bergman, O. Eriksson, L. Bergqvist, and A. Delin, Phys. Rev. B 99, 104302 (2019).
  * Eriksson  _et al._ [2017] O. Eriksson, A. Bergman, L. Bergqvist, and J. Hellsvik,  _Atomistic spin dynamics: foundations and applications_ (Oxford university press, 2017).
  * Duffy and Rutherford [2006] D. Duffy and A. Rutherford, Journal of Physics: Condensed Matter 19, 016207 (2006).
  * Duffy and Rutherford [2009] D. Duffy and A. Rutherford, Journal of nuclear materials 386, 19 (2009).
  * Ma  _et al._ [2010] P.-W. Ma, S. L. Dudarev, A. A. Semenov, and C. H. Woo, Phys. Rev. E 82, 031111 (2010).
  * Casas-Vázquez and Jou [1994] J. Casas-Vázquez and D. Jou, Phys. Rev. E 49, 1040 (1994).
  * Unikandanunni  _et al._ [2021] V. Unikandanunni, R. Medapalli, E. E. Fullerton, K. Carva, P. M. Oppeneer, and S. Bonetti, Applied Physics Letters 118, 232404 (2021).
  * Schoen  _et al._ [2017] M. A. W. Schoen, J. Lucassen, H. T. Nembach, B. Koopmans, T. J. Silva, C. H. Back, and J. M. Shaw, Phys. Rev. B 95, 134411 (2017).
  * Lu  _et al._ [2022] Z. Lu, I. Miranda, S. Streib, M. Pereiro, E. Sjöqvist, O. Eriksson, A. Bergman, D. Thonig, and A. Delin, arXiv https://doi.org/10.48550/arXiv.2211.13486 (2022).
  * Sadhukhan  _et al._ [2022] B. Sadhukhan, A. Bergman, Y. O. Kvashnin, J. Hellsvik, and A. Delin, Phys. Rev. B 105, 104418 (2022).
  * Windsor  _et al._ [2022] Y. W. Windsor, S.-E. Lee, D. Zahn, V. Borisov, D. Thonig, K. Kliemt, A. Ernst, C. Schüßler-Langeheine, N. Pontius, U. Staub, _et al._ , Nature materials 21, 514 (2022).
  * Mishra  _et al._ [2021] K. Mishra, A. Ciuciulkaite, M. Zapata-Herrera, P. Vavassori, V. Kapaklis, T. Rasing, A. Dmitriev, A. Kimel, and A. Kirilyuk, Nanoscale 13, 19367 (2021).
  * Ritzmann  _et al._ [2020] U. Ritzmann, P. M. Oppeneer, and P. Maldonado, Phys. Rev. B 102, 214305 (2020).
  * Bergqvist and Bergman [2018] L. Bergqvist and A. Bergman, Phys. Rev. Materials 2, 013802 (2018).
  * Peduto  _et al._ [1991] P. R. Peduto, S. Frota-Pessa, and M. S. Methfessel, Phys. Rev. B 44, 13283 (1991).
  * Frota-Pessoa [1992] S. Frota-Pessoa, Phys. Rev. B 46, 14570 (1992).
  * 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, _et al._ , Journal of physics: Condensed matter 29, 465901 (2017).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Wills  _et al._ [2010] J. M. Wills, M. Alouani, P. Andersson, A. Delin, O. Eriksson, and O. Grechnyev, Linear muffin-tin orbital method in the atomic sphere approximation, in  _Full-Potential Electronic Structure Method: Energy and Force Calculations with Density Functional and Dynamical Mean Field Theory_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010) pp. 35–46.
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Vanderbilt [1990] D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Dal Corso and de Gironcoli [2000] A. Dal Corso and S. de Gironcoli, Phys. Rev. B 62, 273 (2000).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Marzari  _et al._ [1999] N. Marzari, D. Vanderbilt, A. De Vita, and M. C. Payne, Phys. Rev. Lett. 82, 3296 (1999).
  * Haydock [1980] R. Haydock, in  _Solid state physics_, Vol. 35 (Elsevier, 1980) pp. 215–294.
  * Beer and Pettifor [1984] N. Beer and D. Pettifor, in  _The Electronic Structure of Complex Systems_ (Springer, 1984) pp. 769–777.
  * Billins and Gray [1972] B. H. Billins and D. E. Gray, eds.,  _American Institute of Physics Handbook (3rd Ed)_ (McGraw-Hill, 1972).
  * Von Barth and Hedin [1972] U. Von Barth and L. Hedin, Journal of Physics C: Solid State Physics 5, 1629 (1972).
  * Pajda  _et al._ [2001] M. Pajda, J. Kudrnovský, I. Turek, V. Drchal, and P. Bruno, Phys. Rev. B 64, 174402 (2001).
  * Kvashnin  _et al._ [2015] Y. O. Kvashnin, O. Grånäs, I. Di Marco, M. I. Katsnelson, A. I. Lichtenstein, and O. Eriksson, Phys. Rev. B 91, 125133 (2015).
  * Szilva  _et al._ [2022] A. Szilva, Y. Kvashnin, E. A. Stepanov, L. Nordström, O. Eriksson, A. I. Lichtenstein, and M. I. Katsnelson, arXiv preprint arXiv:2206.02415 (2022).
  * Grechnev  _et al._ [2007] A. Grechnev, I. Di Marco, M. I. Katsnelson, A. I. Lichtenstein, J. Wills, and O. Eriksson, Phys. Rev. B 76, 035107 (2007).
  * Andersen [1975] O. K. Andersen, Phys. Rev. B 12, 3060 (1975).
  * Frota-Pessôa [2004] S. Frota-Pessôa, Phys. Rev. B 69, 104401 (2004).