# Pressure-dependent electronic superlattice in the Kagome-superconductor
CsV3Sb5

F. Stier  These authors contributed equally to this work  Institut für
Festkörper- und Materialphysik, Technische Universität Dresden, 01062 Dresden,
Germany  A.-A. Haghighirad  These authors contributed equally to this work
Institute for Quantum Materials and Technologies, Karlsruhe Institute of
Technology, Kaiserstr. 12, D-76131 Karlsruhe, Germany  G. Garbarino  ESRF, The
European Synchrotron, 71, avenue des Martyrs, CS 40220 F-38043 Grenoble Cedex
9  S. Mishra  Institut für Festkörper- und Materialphysik, Technische
Universität Dresden, 01062 Dresden, Germany  Würzburg-Dresden Cluster of
Excellence ct.qmat, Technische Universität Dresden, 01062 Dresden, Germany  N.
Stilkerich  Institut für Festkörper- und Materialphysik, Technische
Universität Dresden, 01062 Dresden, Germany  Max Planck Institute for Chemical
Physics of Solids, Dresden, Germany  D. Chen  Max Planck Institute for
Chemical Physics of Solids, Dresden, Germany  College of Physics, Qingdao
University, Qingdao 266071, China  C. Shekhar  Max Planck Institute for
Chemical Physics of Solids, Dresden, Germany  T. Lacmann  Institute for
Quantum Materials and Technologies, Karlsruhe Institute of Technology,
Kaiserstr. 12, D-76131 Karlsruhe, Germany  C. Felser  Max Planck Institute for
Chemical Physics of Solids, Dresden, Germany  Würzburg-Dresden Cluster of
Excellence ct.qmat, Technische Universität Dresden, 01062 Dresden, Germany  T.
Ritschel  Institut für Festkörper- und Materialphysik, Technische Universität
Dresden, 01062 Dresden, Germany  J. Geck<EMAIL_ADDRESS>Institut
für Festkörper- und Materialphysik, Technische Universität Dresden, 01062
Dresden, Germany  Würzburg-Dresden Cluster of Excellence ct.qmat, Technische
Universität Dresden, 01062 Dresden, Germany  M. Le Tacon
<EMAIL_ADDRESS>Institute for Quantum Materials and Technologies,
Karlsruhe Institute of Technology, Kaiserstr. 12, D-76131 Karlsruhe, Germany

###### Abstract

We present a high-resolution single crystal x-ray diffraction study of Kagome
superconductor CsV3Sb5, exploring its response to variations in pressure and
temperature. We discover that at low temperatures, the structural modulations
of the electronic superlattice, commonly associated with charge density wave
order, undergo a transformation around $p\sim$ 0.7 GPa from the familiar
$2\times 2$ pattern to a long-range ordered modulation at wavevector
$q=(0,3/8,1/2)$. Our observations align with inferred changes in the CDW
pattern from prior transport and nuclear magnetic resonance studies, providing
new insights into these transitions. Interestingly, the pressure-induced
variations in the electronic superlattice correlate with two peaks in the
superconducting transition temperature as pressure changes, hinting that
fluctuations within the electronic superlattice could be key to stabilizing
superconductivity. However, our findings contrast with the minimal pressure
dependency anticipated by ab initio calculations of the electronic structure.
They also challenge prevailing scenarios based on a Peierls-like nesting
mechanism involving van Hove singularities.

Understanding the correlation between unconventional superconductivity (SC)
and spatial electronic modulations, encompassing magnetic, nematic, or charge
density wave (CDW) orders, represents one of the key challenges in
contemporary condensed matter research [1, 2, 3]. Consequently, extensive
endeavors are underway to unravel the underlying physics in quantum materials
like cuprates [4], nickelates [5, 6], iron or nickel pnictides [7, 8, 9],
transition metal dichalcogenides [10] or heavy fermion compounds [11], where
the relation between SC and CDW phenomena can be examined in different
settings.

In this framework, a new and extremely intriguing class of materials has
recently been discovered: the layered Kagome-compounds AV3Sb5 (A=K, Cs, Rb)
[12]. Characterized by an electronic structure featuring flat bands, van Hove
singularities (vHs), Dirac cones, and non-trivial band topology, these
materials offer a unique platform for exploring novel electronic states of
matter with intertwined orders [3]. Particularly noteworthy is the detection
of an electronic superlattice below 94 K, [12] which is commonly associated
with a CDW, and that is succeeded by superconductivity below $T_{c}\simeq 2.7$
K in CsV3Sb5– a revelation that has sparked considerable excitement. This
excitement has been further substantiated by recent accounts of the absence of
Kohn-anomalies in the phonon dispersion associated with the formation of the
CDW [13, 14], chirality [15], time-reversal symmetry breaking or electronic
nematicity [16] of the electronic superlattice. Consequently, while we refer
to the electronic modulation in CsV3Sb5 as a charge density wave (CDW) in the
following, we recognize that it is often rightly considered unconventional and
that the underlying electronic order may be substantially more complex. As a
matter of fact, numerous important aspects of this CDW and of its relationship
with SC in CsV3Sb5 remain elusive and subject to intense debate [16, 17, 18,
19, 20]. The origin of the CDW and the nature of the interplay between
electronic and lattice degrees of freedom in these materials therefore
requires further elucidation [13, 21, 14, 22].

Figure 1: Temperature and pressure dependence of the superlattice modulation
observed for sample $K$ (top) and $D$ (bottom). $hl$\- and the $hk$-cuts
through XRD intensity distributions reciprocal space are shown. At 40 K and
$p=0.23$ GPa, the two samples both display a $2\times 2$ modulation along $ab$
but a different superstructure along $c$. Increasing the pressure at 40 K
drives both samples into the same 3/8-phase.

To address such issues, pressure tuning has proven to be a valuable tool in
quantum materials research [23, 24, 25], providing a reversible method to
adjust the balance between competing energy scales and enabling the
investigation of complex phase diagrams in these materials without introducing
chemical disorder. This approach is without doubt most relevant in the case of
kagome superconductors, since both SC and CDW transition temperatures have
been reported to strongly depend on hydrostatic pressure, displaying
intriguing behavior already for relatively modest pressures [26, 27, 17]. This
is demonstrated across a diverse array of measurements encompassing
electrical- and magnetotransport, magnetic susceptibility, nuclear magnetic
and quandrupole resonances (NMR and NQR), x-ray diffraction (XRD), and muon
spin resonance [26, 27, 17, 28, 29] which have revealed that both the
superconducting transition temperature ($T_{c}$) and the upper critical field
$H_{c2}$ as a function of pressure $p$ exhibit a distinctive double-dome
feature below $p<2.0$ GPa , characterized by a first maximum of $\sim$ 7 K at
$p_{1}\simeq 0.7$ GPa, succeeded by a second maximum at $p_{2}\simeq 2$ GPa
with $T_{c}\simeq$ 8 K [17, 27, 30, 31]. The CDW transition temperature
appears monotonically suppressed with pressure, but from a marked reduction of
the magnetoresistance [26] and the appearance of new NMR lines above $p_{1}$,
changes in the electronic modulation pattern have been indirectly inferred.

Nonetheless, a thorough examination of the pressure dependence of CDW
superstructures across various temperatures –essential for clarifying the
microscopic physics at play– has yet to be undertaken. Here we close this gap,
by reporting a detailed high resolution single crystal x-ray diffraction (XRD)
study of CsV3Sb5 as function of pressure and temperature. We report the
gradual suppression of the familiar 2x2 CDW modulation with pressure,
alongside the appearance of a new type of modulation with a completely
different pattern and periodicity above $p_{1}$. The two orders coexist for a
while before the disappearance of the former. At 25 K, no CDW survives above
$\sim$ 1.7 GPa. The behavior of these CDWs is compatible with their ultimate
vanishing at a quantum critical point located near $p_{2}$ [27].

However, the striking transformation of the CDW pattern within this pressure
range contrasts with the modest impact on the electronic structure found in
density functional theory (DFT). In particular, the van Hove singularities
below the Fermi energy show minimal change at these pressures. This
discrepancy challenges weak coupling scenarios that attribute the formation of
electronic superlattices in Kagome metals AV3Sb5 (A=K, Cs, Rb) to a Peierls-
like nesting of van Hove singularities [32, 33].

Figure 2: Pressure dependence of the CDW order in CsV3Sb5. The contour plot
shows the XRD-intensity along the $A-L$ reciprocal direction (obtained from
the averaging of the intensity along the six equivalent reciprocal directions
$A-L$ from the (3,0,1/2) A point) measured at 25 K. The intensity at the L
point below 0.7 GPa signals the $2\times 2\times 2$ order, which transforms
into the 3/8 modulation with increasing pressure revealed by the intensity at
$q=3/8$ and $q=5/8$ above 0.7 GPa.

We investigated single crystals from two distinct batches, referred to as
batch $D$ and batch $K$ in the following, synthesized and characterized, as
described in the supplementary information [34]. High-resolution XRD
measurements as a function of $T$ and $p$ were performed at the ID15B beamline
of the ESRF in Grenoble, France, using a monochromatic x-ray beam with an
energy of 30.17 keV, focused down to 2x4 $\mu$m2 at the sample position. The
diffraction patterns were measured with an EIGER2 X 9M CdTe flat panel
detector. Diamond anvil cells (DACs) with helium as a pressure transmitting
medium were used to apply pressures of up to 27 GPa under close to hydrostatic
conditions. In each DAC, one sample of batch $D$ and one of batch $K$ were
loaded into the sample space, so that both samples were in exactly the same
experimental condition during our pressure and temperature dependent XRD
measurements. Further details are given in the supplementary information [34].

At room temperature and ambient pressure, both samples exhibit a hexagonal
$P6/mmm$ structure, in full agreement with previously published results [35].
In the following, reflections will be always indexed referring to this
structure, even though there are structural transitions as a function of $p$
and $T$ [34]. In agreement with previous reports [12], both samples develop a
$2\times 2$ orthorhombic superstructure within the $ab$-planes when cooled
below 94 K at ambient or small pressure. This CDW transition is signalled by
additional satellite reflections, as demonstrated in panels (a) and (c) of
Fig. 1. Interestingly, as can also be observed in these panels, the structural
modulation of the two set of samples differs along the $c$-direction. At 40 K
and $p=0.23$ GPa, despite special care taken to ensure they are precisely in
the same experimental condition and sharing the same $pT$-history, sample $K$
exhibits a $2\times 2\times 4$ modulation, whereas sample $D$ displays a
$2\times 2\times 2$ order. This observation supports earlier conclusions that
the orders along $c$ are metastable [36], while the $2\times 2$ modulation
within the $ab$-layers is a defining feature of the electronic ordering
instability of CsV3Sb5 at ambient pressure [37].

The major experimental finding of our study is represented by the data in the
right panels (b) and (d) of Fig. 1: Irrespective of the periodicity of the
modulation along the $c$-axis, an increase of the pressure $p$ to 0.7 GPa at
constant T=40 K induces in both samples a transition into a new phase
characterized by the same commensurate modulation vector
$q=(0,\frac{3}{8},\frac{1}{2})$. We will refer to this as the ’3/8’ phase in
the following. This observation aligns very well with recent findings from
${}^{51}V$ NMR [17] and ${}^{121/123}Sb$ NQR [27], which reported in this very
pressure range the appearance of new spectral features that are distinct from
those of the $2\times 2$ modulation and have been interpreted as evidence of a
new type of CDW order. Although these local probes cannot directly determine
the ordering wave vector, proposals for an incommensurate structure [27] or a
unidirectional order with a period of $4\times a$ [17] have been suggested.
While our XRD data confirms a pressure-driven change in the electronic order,
we unambiguously reveal here that the $p$-induced ordering is commensurate and
different from the suggested $4\times a$.

With a further increase in $p$ at constant $T$, the intensities of the
superlattice peaks of the 3/8-phase gradually diminish and vanish above 1.7
GPa, again in good agreement with the nuclear resonance data [17, 27].
According to our structural refinement, the 3/8 phase is accompanied with a
small monoclinic distortion [34].

To further elucidate the evolution of the CDW order with respect $p$, we show
in Fig. 2 the pressure dependence of the XRD-intensity along the high symmetry
$A-L$ line (starting from the reciprocal A point at (3,0,0.5)) as a function
of $p$ at constant $T=25$ K. With $p$ increasing towards 0.7 GPa, the
satellite reflections of the $2\times 2$ in-plane order are gradually
suppressed. As can be observed in Fig. 2, the intensity of this peak weakens
significantly around 0.7 GPa, but remains detectable up to 0.9 GPa. At the
same time, the new 3/8 modulation emerges around 0.7 GPa, i.e., the $2\times
2$ and the 3/8 phases coexist between 0.7 and 0.9 GPa. Beyond 0.9 GPa, the
$2\times 2$ disappears completely, and by further increasing $p$ also the
scattering due to the 3/8 modulation is progressively reduced until it falls
below our detection limit at 1.7 GPa. The data shown in Fig. 2 (a) further
reveals that, within the error of our experiment the positions of the
satellite peaks corresponding to the $2\times 2$ and the 3/8 order do not
change with varying $p$ and that, despite the strong suppression of the
intensity, the satellite reflections always remain as sharp as the Bragg
reflections and are essentially resolution limited [34].This is true both for
in- and out-of-plane directions.

Figure 3: Electronic phase diagram of CsV3Sb5 in the low-pressure and low-
temperature region. The orthorhombic $2\times 2\times 2$ and $2\times 2\times
4$ phases, summarized as $2\times 2$, exist in the green region below $\approx
0.7$ GPa. The new monoclinic phase with the 3/8-order is indicated in blue.
Both phases coexist in a narrow pressure range as indicated by the striped
area. The red dashed line sketches the evolution of the superconducting
$T_{c}$ as a function of pressure (after Ref. 17).

The results of our XRD studies below 100 K and up to 2.5 GPa are summarized in
the $pT$-phase diagram presented in Fig. 3, where the stability regions of the
$2\times 2$ and the 3/8 orders are indicated by the green and blue areas,
respectively. Comparing the pressure dependence of the electronic order to the
pressure dependence of the superconducting $T_{c}$ [38, 39] reveals an
interesting correlation: The pressure $p_{1}\simeq 0.7$ GPa, where $T_{c}(p)$
shows its first maximum, aligns perfectly to the critical pressure where the
$2\times 2$ order is suppressed and the 3/8 order appears. Furthermore, the
pressure $p_{2}\simeq 2$ GPa of the second maximum of $T_{c}(p)$ coincides
well with that at which the 3/8 phase disappears. The suppression of the 3/8
phase also further supports the possible existence of quantum critical point
below this dome.

The observation that every time an electronic order is suppressed with
increasing $p$, the superconducting transition temperature $T_{c}(p)$ exhibits
a maximum is remarkable. One obvious explanation to be considered is, of
course, the competition between CDW and SC. Direct experimental evidence for
this competition, for example in the form of a suppression of the CDW in the
superconducting state is still lacking. Nevertheless, the observation of a
peak in $T_{c}(p)$ near the critical point of a distinct electronically
ordered phase (here the 2x2 or the 3/8 phase) suggests that the corresponding
fluctuations may contribute to supporting superconductivity. This conjecture
certainly warrants careful examination in future investigations.

It is important to note that the consistently resolution-limited satellite
reflections (cf. Fig. 2) do not necessarily imply the absence of fluctuations.
The resolution-limited width of these reflections simply indicates that the
coherent contribution to a given satellite reflection originates from a volume
of about $10^{4}$ unit cells . However, it remains indeterminate whether this
volume constitutes a single contiguous entity or a series of smaller, distinct
volumes maintaining a fixed phase relation (see Ref. [40]). In essence,
despite the sharpness of the satellite peaks, regions experiencing
fluctuations may still exist. These fluctuations would manifest as an
exceedingly weak diffuse scattering signal, beyond the detection capabilities
of the current experimental setup.

The comparison of the present $p$-dependent XRD results to those of a recent
doping dependent study [41] is also very interesting: In the doping dependent
study, it was observed that replacing Sb with Sn destabilizes the $2\times 2$
superlattice as well, leading to the emergence of a new superlattice at low
temperatures with a modulation wave vector similar to that of the 3/8 phase
observed here. However, there are notable distinctions between the doping and
$p$-induced orders: The doping-induced phase exhibits short and highly
anisotropic correlation lengths within the $ab$-layers and signs of electronic
nematicity. In contrast, the superlattice modulations observed in pristine
CsV3Sb5, including the $p$-induced 3/8-phase, always display full long-range
order in all three dimensions with an isotropic in-plane correlation length.
The latter is evident from the high-resolution data presented in Fig. 1, and
confirmed by further analysis which revealed that all satellites reflections
are as narrow (resolution limited) as the neighboring Bragg reflections with
correlation lengths amounting to 100 Å both in and out of the $hk$ plane. The
2x2 and the pressure-induced 3/8-order therefore both lack the characteristics
discussed earlier in terms of electronic nematicity.

In principle, Sn-substitution and $p$ do not represent equivalent
modifications of CsV3Sb5. On a qualitative level, however, they both affect
the $c/a$ ratio, which has been found to play a critical role for controlling
the electronic properties [19, 22]. Due to the anisotropic layered structure
of CsV3Sb5, hydrostatic pressure has indeed a pronounced effect on the c/a
ratio, which can be determined precisely by our measurements. We find that a
pressure as modest as 1 GPa is sufficient to reduce c/a by about 3% [34].

In good agreement with DFT calculations [22], our structural refinement
reveals that the distance between apical Sb and the kagome layers is barely
affected by $p=1$ GPa, in stark contrast to the inter unit-cell distance
between apical Sb-sites, which decreases by about 5%. These structural changes
primarily affect the Sb 5$p_{z}$-states, as clearly demonstrated in the DFT
band structures shown in supplementary Fig. 7, which are based on our
structural data. As can be observed in this figure, the bottom of the
$\Gamma$-centered electron pocket and the top of the hole-like band at the
A-point move significantly closer to the Fermi level with increasing $p$.
However, the vHs at the M-point, which are below the Fermi-level, remain
essentially unaffected. This behaviour is in very good agreement with previous
theoretical studies [42, 22, 43], where similar trends have been reported.

The drastic pressure dependence of the electronic superlattice reported here
therefore implies that the nesting of vHs cannot be the primary stabilization
mechanism in CsV3Sb5. If this was the case, the $2\times 2$ order would be
expected to remain entirely stable below 1 GPa, at odds with our experimental
observation. This conclusion is further supported by a detailed analysis of
the Lindhard function from DFT [44], which shows that there is no unique
nesting vector and, in particular, that the vHs do not cause well-defined
peaks in the Lindhard function. Within DFT, the electronic system alone does
therefore not favor a specific ordering wave vector [44]. The emergence of the
3/8 order, which has not been predicted by any first principles calculations
to our knowledge, highlights that the mechanism underlying charge order
instabilities in these kagome materials remains to be understood.

A strongly momentum-dependent electron-phonon coupling has for instance been
considered as a driving force behind the stabilization of a specific long-
range order [45, 44, 46]. The $p$-driven transition from the $2\times 2$ to
the 3/8 phase would then imply a strongly $p$-dependent electron-phonon
coupling at pressures below 1 GPa, which is hardly consistent with weak
coupling models. This hints at strong coupling effects beyond DFT.

In conclusion, the present study yields two major findings namely the
correlation between the critical pressures of the electronic modulations and
superconductivity in CsV3Sb5 and a transition from the $2\times 2$ to a long-
range charge ordered 3/8 phase at pressures below 1 GPa. Furthermore, our
findings reveal pressure-induced structural modulations that appear to
coincide with the trend in the superconductivity manifested in other
experiments. These observations strongly substantiate the notion that critical
fluctuations and complex electron-lattice interactions, extending well beyond
conventional scenarios and models, play a prominent role in shaping the
complex phase diagram of AV3Sb5.

## Acknowledgements

We acknowledge fruitful discussions with R. Fernandes, R. Heid, M. Merz and S.
M. Souliou. This research received support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation), under projects
422213477 (TRR 288 project B03) and 247310070 (CRC 1143 project C06 and B05).
T.R., F.S., C.F. and J.G. express gratitude to the Würzburg-Dresden Cluster of
Excellence on Complexity and Topology in Quantum Matter (ct.qmat EXC-2147,
Project No. 390858490) for financial backing. Additionally, we extend our
appreciation to the European Synchrotron Radiation Facility (ESRF) for
providing access to their synchrotron radiation facilities.

## References

  * Fernandes  _et al._ [2019] R. M. Fernandes, P. P. Orth, and J. Schmalian, Intertwined vestigial order in quantum materials: Nematicity and beyond, Annual Review of Condensed Matter Physics 10, 133 (2019).
  * Agterberg  _et al._ [2020] D. F. Agterberg, J. C. S. Davis, S. D. Edkins, E. Fradkin, D. J. Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, The physics of pair-density waves: Cuprate superconductors and beyond, Annual Review of Condensed Matter Physics 11, 231 (2020).
  * Neupert  _et al._ [2022] T. Neupert, M. M. Denner, J.-X. Yin, R. Thomale, and M. Z. Hasan, Charge order and superconductivity in kagome materials, Nature Physics 18, 137 (2022).
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, From quantum matter to high-temperature superconductivity in copper oxides, Nature 518, 179 (2015).
  * Tam  _et al._ [2022] C. C. Tam, J. Choi, X. Ding, S. Agrestini, A. Nag, M. Wu, B. Huang, H. Luo, P. Gao, M. García-Fernández, L. Qiao, and K.-J. Zhou, Charge density waves in infinite-layer NdNiO2 nickelates, Nature Materials 21, 1116 (2022).
  * Krieger  _et al._ [2022] G. Krieger, L. Martinelli, S. Zeng, L. E. Chow, K. Kummer, R. Arpaia, M. Moretti Sala, N. B. Brookes, A. Ariando, N. Viart, M. Salluzzo, G. Ghiringhelli, and D. Preziosi, Charge and spin order dichotomy in NdNiO2 driven by the capping layer, Physical Review Letters 129, 027002 (2022).
  * Yi  _et al._ [2018] M. Yi, A. Frano, D. H. Lu, Y. He, M. Wang, B. A. Frandsen, A. F. Kemper, R. Yu, Q. Si, L. Wang, M. He, F. Hardy, P. Schweiss, P. Adelmann, T. Wolf, M. Hashimoto, S. K. Mo, Z. Hussain, M. Le Tacon, A. E. Bohmer, D. H. Lee, Z. X. Shen, C. Meingast, and R. J. Birgeneau, Spectral evidence for emergent order in Ba1-xNaxFe2As2, Physical Review Letters 121, 10.1103/PhysRevLett.121.127001 (2018).
  * Souliou  _et al._ [2022] S. M. Souliou, T. Lacmann, R. Heid, C. Meingast, M. Frachet, L. Paolasini, A. A. Haghighirad, M. Merz, A. Bosak, and M. Le Tacon, Soft-phonon and charge-density-wave formation in nematic BaNi2As2, Physical Review Letters 129, 247602 (2022).
  * Lacmann  _et al._ [2023] T. Lacmann, A.-A. Haghighirad, S.-M. Souliou, M. Merz, G. Garbarino, K. Glazyrin, R. Heid, and M. Le Tacon, High-pressure phase diagram of ${\mathrm{bani}}_{2}{\mathrm{as}}_{2}$: Unconventional charge density waves and structural phase transitions, Physical Review B 108, 224115 (2023).
  * Leroux  _et al._ [2015] M. Leroux, I. Errea, M. Le Tacon, S. M. Souliou, G. Garbarino, L. Cario, A. Bosak, F. Mauri, M. Calandra, and P. Rodiere, Strong anharmonicity induces quantum melting of charge density wave in 2h-nbse2 under pressure, Physical Review B 92, 140303 (2015).
  * Aishwarya  _et al._ [2023] A. Aishwarya, J. May-Mann, A. Raghavan, L. Nie, M. Romanelli, S. Ran, S. R. Saha, J. Paglione, N. P. Butch, E. Fradkin, and V. Madhavan, Magnetic-field-sensitive charge density waves in the superconductor ute2, Nature 618, 928 (2023).
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, CsV3Sb5: A ${\mathbb{z}}_{2}$ topological kagome metal with a superconducting ground state, Physical Review Letters 125, 247002 (2020).
  * Li  _et al._ [2021] H. Li, T. T. Zhang, T. Yilmaz, Y. Y. Pai, C. E. Marvinney, A. Said, Q. W. Yin, C. S. Gong, Z. J. Tu, E. Vescovo, C. S. Nelson, R. G. Moore, S. Murakami, H. C. Lei, H. N. Lee, B. J. Lawrie, and H. Miao, Observation of unconventional charge density wave without acoustic phonon anomaly in kagome superconductors AV3Sb5 (A=Rb, Cs), Physical Review X 11, 031050 (2021).
  * Subires  _et al._ [2023] D. Subires, A. Korshunov, A. H. Said, L. Sánchez, B. R. Ortiz, S. D. Wilson, A. Bosak, and S. Blanco-Canosa, Order-disorder charge density wave instability in the kagome metal (Cs,Rb)V3Sb5, Nature Communications 14, 1015 (2023).
  * Guo  _et al._ [2022] C. Guo, C. Putzke, S. Konyzheva, X. Huang, M. Gutierrez-Amigo, I. Errea, D. Chen, M. G. Vergniory, C. Felser, M. H. Fischer, T. Neupert, and P. J. W. Moll, Switchable chiral transport in charge-ordered kagome metal csv3sb5, Nature 611, 461 (2022).
  * Nie  _et al._ [2022] L. Nie, K. Sun, W. Ma, D. Song, L. Zheng, Z. Liang, P. Wu, F. Yu, J. Li, M. Shan, D. Zhao, S. Li, B. Kang, Z. Wu, Y. Zhou, K. Liu, Z. Xiang, J. Ying, Z. Wang, T. Wu, and X. Chen, Charge-density-wave-driven electronic nematicity in a kagome superconductor, Nature 604, 59 (2022).
  * Zheng  _et al._ [2022] L. Zheng, Z. Wu, Y. Yang, L. Nie, M. Shan, K. Sun, D. Song, F. Yu, J. Li, D. Zhao, S. Li, B. Kang, Y. Zhou, K. Liu, Z. Xiang, J. Ying, Z. Wang, T. Wu, and X. Chen, Emergent charge order in pressurized kagome superconductor CsV3Sb5, Nature 611, 682 (2022).
  * Asaba  _et al._ [2023] T. Asaba, A. Onishi, Y. Kageyama, T. Kiyosue, K. Ohtsuka, S. Suetsugu, Y. Kohsaka, T. Gaggl, Y. Kasahara, H. Murayama, K. Hashimoto, R. Tazai, H. Kontani, B. R. Ortiz, S. D. Wilson, Q. Li, H. H. Wen, T. Shibauchi, and Y. Matsuda, Evidence for an odd-parity nematic phase above the charge density wave transition in kagome metal CsV3Sb5 (2023), arXiv:2309.16985 [cond-mat.str-el] .
  * Frachet  _et al._ [2023] M. Frachet, L. Wang, W. Xia, Y. Guo, M. He, N. Maraytta, R. Heid, A.-A. Haghighirad, M. Merz, C. Meingast, and F. Hardy, Colossal c-axis response and lack of rotational symmetry breaking within the kagome plane of the CsV3Sb5 superconductor (2023), arXiv:2310.06102 [cond-mat.supr-con] .
  * Liu  _et al._ [2023] Z. Liu, Y. Shi, Q. Jiang, E. W. Rosenberg, J. M. DeStefano, J. Liu, C. Hu, Y. Zhao, Z. Wang, Y. Yao, D. Graf, P. Dai, J. Yang, X. Xu, and J.-H. Chu, Absence of nematic instability in the kagome metal CsV3Sb5 (2023), arXiv:2309.14574 [cond-mat.supr-con] .
  * Zhong  _et al._ [2023] Y. Zhong, S. Li, H. Liu, Y. Dong, K. Aido, Y. Arai, H. Li, W. Zhang, Y. Shi, Z. Wang, S. Shin, H. N. Lee, H. Miao, T. Kondo, and K. Okazaki, Testing electron–phonon coupling for the superconductivity in kagome metal CsV3Sb5, Nature Communications 14, 1945 (2023).
  * Ritz  _et al._ [2023] E. T. Ritz, H. S. Røising, M. H. Christensen, T. Birol, B. M. Andersen, and R. M. Fernandes, Superconductivity from orbital-selective electron-phonon coupling in AV3Sb5, Phys. Rev. B 108, L100510 (2023).
  * Souliou  _et al._ [2020] S. M. Souliou, A. Bosak, G. Garbarino, and M. L. Tacon, Inelastic x-ray scattering studies of phonon dispersions in superconductors at high pressures, Superconductor Science and Technology 33, 124004 (2020).
  * Pimenta Martins  _et al._ [2023] L. G. Pimenta Martins, R. Comin, M. J. S. Matos, M. S. C. Mazzoni, B. R. A. Neves, and M. Yankowitz, High-pressure studies of atomically thin van der waals materials, Applied Physics Reviews 10, 011313 (2023).
  * Mao  _et al._ [2018] H.-K. Mao, X.-J. Chen, Y. Ding, B. Li, and L. Wang, Solids, liquids, and gases under high pressure, Rev. Mod. Phys. 90, 015007 (2018).
  * Yu  _et al._ [2021] F. H. Yu, D. H. Ma, W. Z. Zhuo, S. Q. Liu, X. K. Wen, B. Lei, J. J. Ying, and X. H. Chen, Unusual competition of superconductivity and charge-density-wave state in a compressed topological kagome metal, Nature Communications 12, 3645 (2021).
  * Feng  _et al._ [2023] X. Y. Feng, Z. Zhao, J. Luo, J. Yang, A. F. Fang, H. T. Yang, H. J. Gao, R. Zhou, and G.-q. Zheng, Commensurate-to-incommensurate transition of charge-density-wave order and a possible quantum critical point in pressurized kagome metal CsV3Sb5, npj Quantum Materials 8, 23 (2023).
  * Mielke  _et al._ [2022] C. Mielke, D. Das, J.-X. Yin, H. Liu, R. Gupta, Y.-X. Jiang, M. Medarde, X. Wu, H. C. Lei, J. Chang, P. Dai, Q. Si, H. Miao, R. Thomale, T. Neupert, Y. Shi, R. Khasanov, M. Z. Hasan, H. Luetkens, and Z. Guguchia, Time-reversal symmetry-breaking charge order in a kagome superconductor, Nature 602, 245 (2022).
  * Guguchia  _et al._ [2023] Z. Guguchia, R. Khasanov, and H. Luetkens, Unconventional charge order and superconductivity in kagome-lattice systems as seen by muon-spin rotation, npj Quantum Materials 8, 41 (2023).
  * Chen  _et al._ [2021a] K. Y. Chen, N. N. Wang, Q. W. Yin, Y. H. Gu, K. Jiang, Z. J. Tu, C. S. Gong, Y. Uwatoko, J. P. Sun, H. C. Lei, J. P. Hu, and J. G. Cheng, Double superconducting dome and triple enhancement of ${T}_{c}$ in the kagome superconductor CsV3Sb5 under high pressure, Physical Review Letters 126, 247001 (2021a).
  * Du  _et al._ [2021a] F. Du, S. Luo, B. R. Ortiz, Y. Chen, W. Duan, D. Zhang, X. Lu, S. D. Wilson, Y. Song, and H. Yuan, Pressure-induced double superconducting domes and charge instability in the kagome metal KV3Sb5, Physical Review B 103, L220504 (2021a).
  * Tan  _et al._ [2021] H. Tan, Y. Liu, Z. Wang, and B. Yan, Charge density waves and electronic properties of superconducting kagome metals, Physical Review Letters 127, 046401 (2021).
  * Kang  _et al._ [2022] M. Kang, S. Fang, J.-K. Kim, B. R. Ortiz, S. H. Ryu, J. Kim, J. Yoo, G. Sangiovanni, D. Di Sante, B.-G. Park, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, S. D. Wilson, J.-H. Park, and R. Comin, Twofold van hove singularity and origin of charge order in topological kagome superconductor CsV3Sb5, Nature Physics 18, 301 (2022).
  * [34] See Supplemental Material at XXX for information on the samples and experimental methods as well as on the detailed structural analysis, which includes also references [47, 48, 49, 50, 51, 51, 52, 53, 54].
  * Ortiz  _et al._ [2019] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. H. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, E. Ertekin, T. M. McQueen, and E. S. Toberer, New kagome prototype materials: discovery of KV3Sb5, RbV3Sb5, and CsV3Sb5, Phys. Rev. Mater. 3, 094407 (2019).
  * Xiao  _et al._ [2023] Q. Xiao, Y. Lin, Q. Li, X. Zheng, S. Francoual, C. Plueckthun, W. Xia, Q. Qiu, S. Zhang, Y. Guo, J. Feng, and Y. Peng, Coexistence of multiple stacking charge density waves in kagome superconductor CsV3Sb5, Physical Review Research 5, L012032 (2023).
  * Ortiz  _et al._ [2021] B. R. Ortiz, S. M. L. Teicher, L. Kautzsch, P. M. Sarte, N. Ratcliff, J. Harter, J. P. C. Ruff, R. Seshadri, and S. D. Wilson, Fermi surface mapping and the nature of charge-density-wave order in the kagome superconductor CsV3Sb5, Phys. Rev. X 11, 041030 (2021).
  * Chen  _et al._ [2021b] K. Y. Chen, N. N. Wang, Q. W. Yin, Y. H. Gu, K. Jiang, Z. J. Tu, C. S. Gong, Y. Uwatoko, J. P. Sun, H. C. Lei, J. P. Hu, and J.-G. Cheng, Double superconducting dome and triple enhancement of ${T}_{c}$ in the kagome superconductor CsV3Sb5 under high pressure, Phys. Rev. Lett. 126, 247001 (2021b).
  * Du  _et al._ [2021b] F. Du, S. Luo, B. R. Ortiz, Y. Chen, W. Duan, D. Zhang, X. Lu, S. D. Wilson, Y. Song, and H. Yuan, Pressure-induced double superconducting domes and charge instability in the kagome metal KV3Sb5, Phys. Rev. B 103, L220504 (2021b).
  * Hossain  _et al._ [2013] M. A. Hossain, I. Zegkinoglou, Y. D. Chuang, J. Geck, B. Bohnenbuck, A. G. C. Gonzalez, H. H. Wu, C. Schüßler-Langeheine, D. G. Hawthorn, J. D. Denlinger, R. Mathieu, Y. Tokura, S. Satow, H. Takagi, Y. Yoshida, Z. Hussain, B. Keimer, G. A. Sawatzky, and A. Damascelli, Electronic superlattice revealed by resonant scattering from random impurities in Sr3Ru2O7, Scientific Reports 3, 2299 (2013).
  * Kautzsch  _et al._ [2023] L. Kautzsch, Y. M. Oey, H. Li, Z. Ren, B. R. Ortiz, G. Pokharel, R. Seshadri, J. Ruff, T. Kongruengkit, J. W. Harter, Z. Wang, I. Zeljkovic, and S. D. Wilson, Incommensurate charge-stripe correlations in the kagome superconductor CsV3Sb5-xSnx, npj Quantum Materials 8, 37 (2023).
  * Wenzel  _et al._ [2023] M. Wenzel, A. A. Tsirlin, F. Capitani, Y. T. Chan, B. R. Ortiz, S. D. Wilson, M. Dressel, and E. Uykur, Pressure evolution of electron dynamics in the superconducting kagome metal CsV3Sb5, npj Quantum Materials 8, 45 (2023).
  * Tsirlin  _et al._ [2022] A. A. Tsirlin, P. Fertey, B. R. Ortiz, B. Klis, V. Merkl, M. Dressel, S. D. Wilson, and E. Uykur, Role of sb in the superconducting kagome metal CsV3Sb5 revealed by its anisotropic compression, SciPost Phys. 12, 049 (2022).
  * Kaboudvand  _et al._ [2022] F. Kaboudvand, S. M. L. Teicher, S. D. Wilson, R. Seshadri, and M. D. Johannes, Fermi surface nesting and the lindhard response function in the kagome superconductor CsV3Sb5, Applied Physics Letters 120, 111901 (2022).
  * Xie  _et al._ [2022] Y. Xie, Y. Li, P. Bourges, A. Ivanov, Z. Ye, J.-X. Yin, M. Z. Hasan, A. Luo, Y. Yao, Z. Wang, G. Xu, and P. Dai, Electron-phonon coupling in the charge density wave state of CsV3Sb5, Phys. Rev. B 105, L140501 (2022).
  * Ye  _et al._ [2022] Z. Ye, A. Luo, J.-X. Yin, M. Z. Hasan, and G. Xu, Structural instability and charge modulations in the kagome superconductor $a{\mathrm{v}}_{3}{\mathrm{sb}}_{5}$, Phys. Rev. B 105, 245121 (2022).
  * [47] Crysalispro single crystal x-ray diffraction data collection and processing software.
  * Sheldrick [2015] G. M. Sheldrick, SHELXT – Integrated space-group and crystal-structure determination, Acta Crystallographica Section A 71, 3 (2015).
  * Petříček  _et al._ [2014] V. Petříček, M. Dušek, and L. Palatinus, Crystallographic Computing System JANA2006: General features, Zeitschrift für Kristallographie - Crystalline Materials 229, 345 (2014).
  * Stahl  _et al._ [2022] Q. Stahl, D. Chen, T. Ritschel, C. Shekhar, E. Sadrollahi, M. C. Rahn, O. Ivashko, M. v. Zimmermann, C. Felser, and J. Geck, Temperature-driven reorganization of electronic order in CsV3Sb5, Phys. Rev. B 105, 195136 (2022).
  * Zhang  _et al._ [2021] Z. Zhang, Z. Chen, Y. Zhou, Y. Yuan, S. Wang, J. Wang, H. Yang, C. An, L. Zhang, X. Zhu, Y. Zhou, X. Chen, J. Zhou, and Z. Yang, Pressure-induced reemergence of superconductivity in the topological kagome metal CsV3Sb5, Phys. Rev. B 103, 224513 (2021).
  * Yu  _et al._ [2022] F. Yu, X. Zhu, X. Wen, Z. Gui, Z. Li, Y. Han, T. Wu, Z. Wang, Z. Xiang, Z. Qiao, J. Ying, and X. Chen, Pressure-induced dimensional crossover in a kagome superconductor, Phys. Rev. Lett. 128, 077001 (2022).
  * Tsirlin  _et al._ [2023] A. A. Tsirlin, B. R. Ortiz, M. Dressel, S. D. Wilson, S. Winnerl, and E. Uykur, Effect of nonhydrostatic pressure on the superconducting kagome metal CsV3Sb5, Phys. Rev. B 107, 174107 (2023).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Phys. Rev. B 59, 1743 (1999).


# Multiparticle Collision Dynamics for Ferrofluids

Patrick Ilg  School of Mathematical, Physical, and Computational Sciences,
University of Reading, Reading, RG6 6AX, United Kingdom<EMAIL_ADDRESS>
###### Abstract

Detailed studies of the intriguing field-dependent dynamics and transport
properties of confined flowing ferrofluids require efficient mesoscopic
simulation methods that account for fluctuating ferrohydrodynamics. Here, we
propose such a new mesoscopic model for the dynamics and flow of ferrofluids,
where we couple the multi-particle collision dynamics method as a solver for
the fluctuating hydrodynamics equations to the stochastic magnetization
dynamics of suspended magnetic nanoparticles. This hybrid model is validated
by reproducing the magnetoviscous effect in Poiseuille flow, obtaining the
rotational viscosity in quantitative agreement with theoretical predictions.
We also illustrate the new method for the benchmark problem of flow around a
square cylinder. Interestingly, we observe that the length of the
recirculation region is increased whereas the drag coefficient is decreased in
ferrofluids when an external magnetic field is applied, compared with the
field-free case at the same effective Reynolds number. The presence of thermal
fluctuations and the flexibility of this particle-based mesoscopic method
provides a promising tool to investigate a broad range of flow phenomena of
magnetic fluids and could also serve as an efficient way to simulate solvent
effects when colloidal particles are immersed in ferrofluids.

##  I Introduction

Colloidal suspensions of magnetic nanoparticles, also known as ferrofluids,
are fascinating model systems that combine superparamagnetic and
magnetoviscous effects Rosensweig (1985); Odenbach (2009). The classical
experiments of McTague McTague (1969) demonstrated that pipe flow of
ferrofluids can be manipulated by external magnetic fields, showing
anisotropic response depending on the relative orientation of the magnetic
field with respect to the flow direction. Shliomis and co-workers Martsenyuk,
Raikher, and Shliomis (1974) successfully explained these experiments by a
kinetic theory of hindered particle rotations and the resulting changes in
rotational viscosity. These ground-breaking works sparked numerous subsequent
studies on the magnetoviscous effect, leading to more refined experiments and
a better theoretical understanding, as well as several novel applications of
ferrofluids that rely on their field-dependent flow behaviour (see e.g.
Odenbach (2002); Kröger, Ilg, and Hess (2003); Colombo  _et al._ (2012); Ilg
and Odenbach (2008); Felicia, Vinod, and Philip (2016) and references
therein). Also various computer simulation methods for ferrofluid flow have
been developed in recent years. The simulations range from simple channel or
pipe flows Rinaldi and Zahn (2002); Schumacher  _et al._ (2003); Hirabayashi,
Chen, and Ohashi (2001) to more complicated geometries Sheikholeslami  _et
al._ (2018) including free-surface flows Huang, Hädrich, and Michels (2019).
Different numerical methods were used in these studies, ranging from
perturbative solutions Rinaldi and Zahn (2002) to adaption of finite-element
and finite-volume computational fluid dynamics codes Schumacher  _et al._
(2003); Sheikholeslami  _et al._ (2018), the Lattice Boltzmann method
Hirabayashi, Chen, and Ohashi (2001), and smooth particle hydrodynamics
simulations Huang, Hädrich, and Michels (2019). Such flow simulations are
helpful for various applications, e.g. planning magnetic drug targeting
treatments Kayal  _et al._ (2011). Also micro magnetofluidics shows promising
potential for contactless mixing, separation, and trapping of particles and
polymers on small scales, where flow simulations help to improve their
effectiveness Munaz, Shiddiky, and Nguyen (2018).

Thermal fluctuations are important in small-scale flows which, however, are
neither captured by traditional computational fluid dynamics codes, nor in the
Lattice Boltzmann method. Instead, multi-particle collision (MPC) dynamics
provides a very versatile and flexible method for simulating fluid flow
including thermal fluctuations Malevanets and Kapral (1999); Ihle and Kroll
(2001); Gompper  _et al._ (2009). In this off-lattice method, each particle
represents a small fluid element and evolves through a sequence of streaming
and simplified collision steps. By locally conserving mass, momentum and
energy, the correct hydrodynamic behavior is generated. In addition, thermal
fluctuations are naturally present in this particle-based scheme. Recently,
the MPC method has been extended to model anisotropic fluids like nematic
liquid crystals Shendruk and Yeomans (2015); Lee and Mazza (2015); Mandal and
Mazza (2019). Here, we use similar reasoning to extend the MPC method to model
ferrofluid flow including stochastic magnetization dynamics. We note that none
of the previous approaches to ferrofluid flow simulations Rinaldi and Zahn
(2002); Schumacher  _et al._ (2003); Sheikholeslami  _et al._ (2018);
Hirabayashi, Chen, and Ohashi (2001); Huang, Hädrich, and Michels (2019)
includes thermal fluctuations in the magnetization dynamics. For simplicity,
we restrict ourselves to a quasi-two-dimensional system, i.e. two-dimensional
translational motion with three-dimensional magnetization dynamics. We
demonstrate that this model is able to reproduce fluctuating
ferrohydrodynamics. In particular, we verify that the field-dependent
effective viscosity resulting from this hybrid method quantitatively agrees
with theoretical predictions. Furthermore, we show the flexibility of this
approach by simulating the flow around a square cylinder. Analyzing the flow
field, we determine the length of the recirculation region and find it
decreasing with increasing field strength. Moreover, we also calculate the
drag coefficient and study its dependence on the magnetic field.

The new hybrid method presented here is a useful tool that also allows to
study the effect of fluctuations on ferrofluid flow on a mesoscopic scale.
Different models of magnetization dynamics can be incorporated in a
straightforward manner. Thereby, the flexibility of the MPC method also
facilitates the study of more complicated geometries. This paper is organized
as follows. In Sect. II, we briefly review the ferrohydrodynamic theory. The
hybrid MPC model coupled to stochastic magnetization dynamics is described in
Sect. III. Sect. IV describes the simulation setup, before results are
presented and discussed for channel flows in Sect. V and for flow around
square cylinder in Sect. VI. Finally some conclusions are offered in Sect.
VII.

##  II Ferrohydrodynamics

We here give a brief overview of the basic equations of ferrohydrodynamics to
make the paper self-contained. Further details can be found e.g. in Refs.
Rosensweig (1985, 2002); Shliomis (2002). The fluid momentum balance equation
reads

$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\bm{\nabla}p/\rho+\nu\nabla^{2}\mathbf{v}+\mathbf{f}_{\mathbf{M}}$
(1)

with $\mathbf{v}$ the fluid velocity field, $\mathrm{d}\mathbf{v}/\mathrm{d}t$
denotes the material derivative, $\rho$ and $\nu$ are the fluid density and
kinematic viscosity, respectively, and $p$ the scalar pressure. With allowance
for internal rotations, the influence of magnetization effects on fluid flow
is described by the force density

$\rho\mathbf{f}_{\mathbf{M}}=(\mathbf{M}\cdot\bm{\nabla})\mathbf{H}+\frac{1}{2}\bm{\nabla}\times(\mathbf{M}\times\mathbf{H})$
(2)

where the magnetic field and the magnetization are denoted by $\mathbf{H}$ and
$\mathbf{M}$, respectively. We assume the fluid to be incompressible and non-
conducting, therefore

$\bm{\nabla}\cdot\mathbf{v}=0,\quad\bm{\nabla}\times\mathbf{H}={\bf
0},\quad\bm{\nabla}\cdot{\bf B}=0,$ (3)

where ${\bf B}=\mu_{0}(\mathbf{H}+\mathbf{M})$ denotes the magnetic induction
and $\mu_{0}$ the permeability of free space.

Equations (1) – (3) are not closed due to the appearance of the magnetization
$\mathbf{M}$. Some of the previous simulation approaches have assumed quasi-
equilibrium conditions and locally approximate
$\mathbf{M}\approx\mathbf{M}_{\rm eq}(\mathbf{H})$ Sheikholeslami  _et al._
(2018); Huang, Hädrich, and Michels (2019). Such approaches, however, neglect
relaxation phenomena and a more accurate treatment of the magnetization
dynamics is often desirable. Unfortunately, despite many efforts and long
debates in the literature about the appropriate form of the magnetization
equation for ferrofluids, no consensus has been reached yet (see e.g. Refs.
Shliomis (2002); Leschhorn and Lücke (2006); Liu and Stierstadt (2008);
Rosensweig (2002); Ilg and Odenbach (2008) and references therein).

One of the advantages of the present method is that it can handle different
magnetization equations rather straightforwardly. Since the scheme naturally
includes thermal fluctuations in the hydrodynamic variables, we here choose to
implement a mesoscopic model that also accounts for thermal fluctuations in
the magnetization. In particular, we employ the classical model of ferrofluid
dynamics of rigid dipoles proposed in Ref. Martsenyuk, Raikher, and Shliomis
(1974) and widely studied since Ilg and Odenbach (2008); Ilg, Kröger, and Hess
(2002); Soto-Aquino, Rosso, and Rinaldi (2011). In the free-draining limit and
neglecting inter-particle interactions, the torque balance of viscous,
magnetic, and Brownian torques reads

$-\zeta(\bm{\omega}_{i}-\bm{\Omega}(\mathbf{r}_{i}))+\mu\mathbf{u}_{i}\times\mathbf{H}(\mathbf{r}_{i})+{\bf
T}^{\rm B}_{i}={\bf 0},$ (4)

where $\mu\mathbf{u}_{i}$ and $\bm{\omega}_{i}$ are the magnetic moment and
angular velocity of magnetic nanoparticle $i$, respectively. The magnitude of
the magnetic moment is denoted by $\mu$. The rotational viscous torque is
proportional to the rotational friction coefficient $\zeta$ and arises when
the particle rotation does not match the local vorticity,
$\bm{\Omega}=(1/2)\bm{\nabla}\times\mathbf{v}$. The magnetic torque is
generated by the magnetic field $\mathbf{H}$ at the location of the
nanoparticle. Brownian fluctuations are modelled by random torques ${\bf
T}^{\rm B}_{i}$ with Gaussian white noise, $\langle{\bf T}^{\rm
B}_{i}\rangle=0$ and $\langle{\bf T}^{\rm B}_{i}(t){\bf T}^{\rm
B}_{j}(t^{\prime})\rangle=2k_{\rm B}T\zeta\delta_{ij}\delta(t-t^{\prime}){\bf
I}$, $k_{\rm B}$ and $T$ Boltzmann’s constant and temperature, respectively,
and ${\bf I}$ the three-dimensional unit matrix. The fluid magnetization is
calculated from local averages over the magnetic moments of the nanoparticles
as
$\mathbf{M}(\mathbf{r},t)=\langle\sum_{i}\mu\mathbf{u}_{i}(t)w(\mathbf{r}-\mathbf{r}_{i}(t))\rangle$,
where $w$ is a suitable weight function. The magnetization arising from an
applied field is responsible for the backflow effect via the magnetic force
density given in Eq. (2). In equilibrium, the model describes
superparamagnetic behavior, $\mathbf{M}_{\rm eq}\sim
L_{1}(H)\mathbf{H}/|\mathbf{H}|$, with $L_{1}$ the Langevin function
Rosensweig (1985).

##  III Model and Numerical implementation

MPC is a very versatile and flexible, particle-based method to simulate fluid
flow. The method has recently been extended to complex fluids Gompper  _et
al._ (2009) like polymer solutions Kowalik and Winkler  and in particular also
to anisotropic fluids like nematic liquid crystals Shendruk and Yeomans
(2015); Lee and Mazza (2015); Mandal and Mazza (2019). Inspired by these
recent works, we here present an application to another class of anisotropic
fluids, namely ferrofluids.

Within MPC, the fluid is represented by a system of $N$ identical particles,
each with mass $m_{i}=m$, where $\mathbf{r}_{i}$ and $\mathbf{v}_{i}$ denote
the position and velocity of particle $i$ and $i=1,\ldots,N$. Here, we follow
previous works Kowalik and Winkler ; Shendruk and Yeomans (2015); Lee and
Mazza (2015); Mandal and Mazza (2019) and let an MPC particle represent a
coarse-grained description of a fluid element, containing solvent as well as
solute particles. Thus, each MPC particle also carries a magnetic moment
$\mu\mathbf{u}_{i}$, where $\mu$ denotes the magnitude and the three-
dimensional unit vector $\mathbf{u}_{i}$ describes the orientation of the
magnetic moment of particle $i$. The basic idea underlying MPC is that
particle dynamics can be split into a streaming and a simplified collision
step, both can be computed very efficiently. For anisotropic fluids, the
translational motion needs to be coupled to rotational dynamics as detailed
below for the case of ferrofluids.

In the streaming step, the positions and velocities of the particles are
updated according to Jonathan K Whitmer (2010)

$\displaystyle\mathbf{r}_{i}(t+\Delta t)$
$\displaystyle=\mathbf{r}_{i}(t)+\Delta t\,\mathbf{v}_{i}(t)+\frac{\Delta
t^{2}}{2m_{i}}\,{\bf F}_{i}(t)$ (5) $\displaystyle\mathbf{v}^{\prime}_{i}(t)$
$\displaystyle=\mathbf{v}_{i}(t)+\frac{\Delta t}{m_{i}}\,{\bf F}_{i}(t)$ (6)

with time step $\Delta t$. The force ${\bf F}_{i}=\mathbf{f}_{\rm
ext}+\mathbf{f}_{\mathbf{M}}(\mathbf{r}_{i})$ acting on particle $i$ consists
of a uniform external forcing $\mathbf{f}_{\rm ext}$ (e.g. due to an applied
pressure gradient) and the magnetic force $\mathbf{f}_{\mathbf{M}}$ defined in
Eq. (2).

The collision step ensures momentum exchange between particles and can be
expressed as Lee and Mazza (2015)

$\mathbf{v}_{i}(t+\Delta t)=\mathbf{V}_{\\!C_{i}}(t)+\beta_{\rm
th}\mathbf{R}\cdot[\mathbf{v}^{\prime}_{i}(t)-\mathbf{V}_{\\!C_{i}}(t)]$ (7)

where $\mathbf{V}_{\\!C_{i}}$ is the center-of-mass velocity of the collision
cell to which particle $i$ belongs, $\mathbf{V}_{\\!C_{i}}=\sum_{j\in{\cal
C}_{i}}m_{j}\mathbf{v}^{\prime}_{j}/M_{C_{i}}$. The total mass in this
collision cell is given by $M_{C_{i}}=\sum_{j\in{\cal C}_{i}}m_{j}$, and
${\cal C}_{i}$ labels all particles currently located in collision cell
$C_{i}$. In MPC, the collision step (7) is performed simultaneously between
all particles $j$ currently residing in the same collision cell. These cells
are defined by a regular square grid of length $a$. Collisions lead to
rotation of the relative velocities, which is described by the rotation matrix
${\mathbf{R}}$ in Eq. (7). In two spatial dimensions considered here, ${\bf
R}$ is completely specified by the rotational angle $\alpha$. Thus, ${\bf R}$
rotates the relative velocities by an angle $\pm\alpha$ with equal
probabilities.

As has been emphasized before Gompper  _et al._ (2009), angular momentum
conservation is important for the rotational dynamics of anisotropic fluids,
but is in general violated by Eq. (7). Angular momentum conservation can be
restored for this scheme, however, if the rotation angle $\alpha$ is not
fixed, but chosen in collision cell $C_{i}$ according to

$\cos\alpha=\frac{A_{1}^{2}-A_{2}^{2}}{A_{1}^{2}+A_{2}^{2}},\quad\sin\alpha=-\frac{2A_{1}A_{2}}{A_{1}^{2}+A_{2}^{2}}$
(8)

where $A_{1}=\sum_{j\in{\cal
C}_{i}}[\mathbf{r}_{j}\times\tilde{\mathbf{v}}_{j}]_{z}$,
$A_{2}=\sum_{j\in{\cal C}_{i}}\mathbf{r}_{j}\cdot\tilde{\mathbf{v}}_{j}$, with
relative velocities
$\tilde{\mathbf{v}}_{i}=\mathbf{v}^{\prime}_{i}-\mathbf{V}_{\\!C_{i}}$ Gompper
_et al._ (2009).

In flow simulations, a thermostat is generally needed in order to remove the
energy input into the system. Here, we follow common practice (see e.g. Ref.
Lee and Mazza (2015)) and use a simple rescaling of the relative velocities to
the bath temperature $T$ by $\beta_{\rm th}=\sqrt{T/T_{C_{i}}}$. The
instantaneous kinetic temperature in collision cell $C_{i}$ is defined in two
spatial dimensions by $k_{\rm B}T_{C_{i}}=(1/2N_{C_{i}})\sum_{j\in{\cal
C}_{i}}m_{j}\tilde{\mathbf{v}}_{j}^{2}$, with $N_{C_{i}}=\sum_{j\in{\cal
C}_{i}}1$ the number of particles in the collision cell.

Using a fixed grid of collision cells violates Galilean invariance and might
lead to spurious correlations. Therefore, following Ref. Ihle and Kroll
(2001), in each time step, the grid is shifted by a random vector, where each
component is uniformly distributed in $[-a/2,a/2]$.

The simulations discussed below employ two types of boundary conditions. Bulk
simulations employ periodic boundary conditions in all spatial dimensions. For
flow in a channel geometry, we employ periodic boundary conditions only in the
flow direction, while no-slip boundary conditions on the confining walls are
realized by the well-known bounce back rule, where the velocity of the
particles are reversed as a result of collision with the wall Jonathan K
Whitmer (2010). Since random grid shifting can lead to underpopulated cells at
the walls that would undermine the no-slip condition, additional, so-called
ghost particles are added in the collision step Jonathan K Whitmer (2010).

So far, we have specified the translational motion only. To simulate
ferrofluid flow, we also need to specify the dynamics of the magnetic moments
of the particles. From the kinematic equation for rotations,
$\dot{\mathbf{u}}_{i}=\bm{\omega}_{i}\times\mathbf{u}_{i}$, and Eq. (4), a
weak first-order scheme can be used, $\mathbf{u}_{i}(t+\Delta t_{\rm
B})=\mathbf{u}_{i}^{\rm P}/|\mathbf{u}_{i}^{\rm P}|$, where
$\mathbf{u}_{i}^{\rm
P}=\mathbf{u}_{i}(t)+\Delta\bm{\omega}_{i}\times\mathbf{u}_{i}(t)$ with

$\Delta\bm{\omega}_{i}=[\tau_{\rm
B}\bm{\Omega}_{C_{i}}+\frac{1}{2}\mathbf{u}_{i}\times\mathbf{h}_{C_{i}}]\frac{\Delta
t_{\rm B}}{\tau_{\rm B}}+\frac{1}{\sqrt{\tau_{\rm B}}}\Delta\mathbf{W}_{i}$
(9)

In Eq. (9) we introduced the dimensionless magnetic field, which is defined as
the ratio of the local Zeeman energy over the thermal energy,
$\mathbf{h}_{C_{i}}=\mu\mathbf{H}(\mathbf{r}_{C_{i}})/k_{\rm B}T$.
Furthermore, $\Delta\mathbf{W}_{i}$ denotes the Wiener increment over a time
interval $\Delta t_{\rm B}$ and $\tau_{\rm B}=\zeta/(2k_{\rm B}T)$ the
Brownian rotational relaxation time of a single nanoparticle. Equation (9)
requires $\Delta t_{\rm B}/\tau_{\rm B}\ll 1$. For the simulations shown
below, we use a second-order stochastic Heun algorithm. More details on the
algorithm are given in appendix A.

In Eq. (9), we denote the local vorticity and local magnetic field with
$\bm{\Omega}_{C_{i}}$ and $\mathbf{h}_{C_{i}}$, respectively, to indicate that
these fields are not evaluated on the precise location of particle $i$, but at
the center of the collision cell to which particle $i$ belongs at this time
step. Similarly, we define the instantaneous magnetization of the collision
cell as $\mathbf{M}_{C_{i}}=(M_{\rm sat}/N_{C_{i}})\sum_{j\in{\cal
C}_{i}}\mathbf{u}_{j}$, with the saturation magnetization $M_{\rm sat}=n\mu$
and $n$ the number density of magnetic nanoparticles. For the mean
magnetization, we recover the usual relation
$\langle\mathbf{M}_{C_{i}}\rangle=M_{\rm sat}\langle\mathbf{u}\rangle$.

The magnetization dynamics (4) is coupled to the flow field via the vorticity,
i.e. gradients of the velocity field,
$\bm{\Omega}=(1/2)\bm{\nabla}\times\mathbf{v}$, whereas gradients in the
magnetization or magnetic field lead to backflow effects via the magnetic
force density in the momentum balance, Eq. (2). In order to compute these
gradients, we evaluate the velocity and magnetization field in each time step
via a kernel smoothing method and determine spatial gradients from a finite-
difference approximation. Further details on the numerical procedure are
provided in appendix B.

##  IV Simulation setup

Consider a simple channel geometry of width $L$ with infinite plates
perpendicular to the $y$-direction, as sketched in Fig. 1. The externally
applied field is denoted by $\mathbf{H}_{0}$. We assume the exterior to be
non-magnetic such that $\mathbf{M}_{0}={\bf 0}$ with the external magnetic
induction $\mathbf{B}_{0}=\mu_{0}\mathbf{H}_{0}$. The magnetic field inside
the channel is denoted by $\mathbf{H}$. The presence of ferrofluid gives rise
to a magnetization $\mathbf{M}$, such that the magnetic induction inside the
channel is given by $\mathbf{B}=\mu_{0}(\mathbf{H}+\mathbf{M})$.

If $\hat{\mathbf{k}}$ denotes the normal vector on the interface, the
continuity conditions of magnetostatics require the normal component of the
magnetic induction $\mathbf{B}$ and the tangential component of the magnetic
field $\mathbf{H}$ are continuous at the interface Rosensweig (1985)

$\hat{\mathbf{k}}\cdot(\mathbf{B}_{0}-\mathbf{B})=0,\quad\hat{\mathbf{k}}\times(\mathbf{H}_{0}-\mathbf{H})={\bf
0}.$ (10)

The internal field can be expressed as
$\mathbf{H}=\mathbf{H}_{0}-\bar{\mathbf{D}}\cdot\mathbf{M}$ where the
demagnetization tensor $\bar{\mathbf{D}}$ depends on the geometry. For the
simple channel geometry of Fig. 1 with infinite plates and a spatially
homogeneous external field $\mathbf{H}_{0}$, we obtain $H_{x}=H_{0,x}$ and
$H_{y}=H_{0,y}-M_{y}$. The same result is obtained for an infinite needle-
shaped ellipsoidal geometry Osborn (1945). It is important to note that
Maxwell’s equation together with the continuity conditions (10) require for
our case the magnetization $\mathbf{M}$ to be independent of $x$-position,
thus $\mathbf{M}=\mathbf{M}(y)$ can depend only on the vertical position
within the channel. Further details are provided in Appendix C. In terms of
the dimensionless field $\mathbf{h}$, the relation between external and
internal field can be written as $h_{y}=h_{0,y}-3\chi_{\rm L}\langle
u_{y}\rangle$ with the Langevin susceptibility $\chi_{\rm L}=n\mu^{2}/(3k_{\rm
B}T)$.

(a)$x$$y$$L$$0$$\mathbf{H},\mathbf{M}$$\mathbf{H}_{0}$$\hat{\bf
k}$$\mathbf{v}$

(b)$L_{x}$$L_{x}/4$$L$$D$

Figure 1: (a) Schematic of channel geometry and velocity profile. (b) Channel
including square cylinder in Sect. VI.

In order to study Poiseuille flow, we model an applied pressure gradient as a
homogeneous external forcing in $x$-direction, $\mathbf{f}_{\rm ext}=f_{\rm
ext}\hat{\bf i}$. This constant force is applied to all particles, see Eq.
(5). For the channel geometry of Fig. 1, the dimensionless magnetic force
density $\mathbf{f}_{\mathbf{M}}^{\ast}=(\Delta
t^{2}/a)\mathbf{f}_{\mathbf{M}}$ from Eq. (2) becomes

$\displaystyle\mathbf{f}_{\mathbf{M}}^{\ast}$
$\displaystyle=-n^{\ast}\left(\begin{array}[]{c}\frac{1}{2}[(3\chi_{\rm
L}\langle u_{y}\rangle-h_{0,y})\frac{\mathrm{d}\langle
u_{x}\rangle}{\mathrm{d}y^{\ast}}+(h_{0,x}+3\chi_{\rm L}\langle
u_{x}\rangle)\frac{\mathrm{d}\langle u_{y}\rangle}{\mathrm{d}y^{\ast}}]\\\
3\chi_{\rm L}\langle u_{y}\rangle\frac{\mathrm{d}\langle
u_{y}\rangle}{\mathrm{d}y^{\ast}}\end{array}\right)$ (13)

where $n^{\ast}=nk_{\rm B}T/p_{\rm ref}$ is a density ratio and $p_{\rm
ref}=\rho a^{2}/\Delta t^{2}$ a reference pressure. In Eq. (13), we used the
fact that the magnetostatic fields are independent of the $x$-coordinate.

We follow common practice in MPC simulations and choose particle masses as
$m_{i}=1$, the linear size of collision cell $a=1$, and time step $\Delta
t=1$. We note that we limit ourselves to small forcings, otherwise smaller
values of $\Delta t$ are needed. Having set the basic MPC units, table 1 gives
an overview of the remaining model parameters and the values or range of
values used in subsequent simulations. A number of comments are in order.
Temperature is measured in units of $T_{\rm ref}=mv_{\rm max}^{2}/k_{\rm B}$,
with $v_{\rm max}=a/\Delta t$ the maximum propagation speed. The mean free
path is defined as $\lambda=\Delta t\sqrt{k_{\rm B}T/m}$. We choose
$T^{\ast}=T/T_{\rm ref}<1$ to ensure the mean free path is smaller than the
grid size. The mean number of MPC particles per collision cell $Q$ is chosen
between $20$ and $100$. Furthermore, we choose $\tau_{\rm B}^{\ast}=100$ to
model slower relaxation of nanoparticles compared to the base fluid. This then
allows us to set $\Delta t_{\rm B}=\Delta t$, satisfying $\Delta t_{\rm
B}/\tau_{\rm B}=10^{-2}\ll 1$.

Since every MPC particle represents a small volume element of fluid, there is
no contradiction of choosing $n^{\ast}$ independent of $Q$. Although particles
in the original MPC scheme are point-like, introducing the rotational friction
coefficient $\zeta$ in Eq. (4) allows us to associate an effective
hydrodynamic diameter $\sigma$ to the particles via Stokes’ formula
$\zeta=\pi\eta_{\rm s}\sigma^{3}$, where $\eta_{\rm s}$ denotes the viscosity
of the solvent. Here, we refer to the pure MPC fluid as “solvent”. Therefore,
it is possible to express the density ratio $n^{\ast}$ in terms of the volume
fraction of magnetic particles $\phi$ as $n^{\ast}=3\phi\nu_{\rm
s}^{\ast}/\tau_{\rm B}^{\ast}$, where we used the definition of the reference
pressure $p_{\rm ref}$ and the dimensionless kinetic viscosity $\nu_{\rm
s}^{\ast}=\nu_{\rm s}\Delta t/a^{2}$ with $\nu_{\rm s}=\eta_{\rm s}/\rho$.
Therefore, even though every particle in the simulation carries a magnetic
moment $\mu$, the condition $n^{\ast}\tau_{\rm B}^{\ast}/\nu_{\rm s}^{\ast}\ll
1$ still corresponds to dilute suspensions of magnetic nanoparticles, which,
strictly speaking, is needed for justifying the magnetization dynamics
described by Eq. (4).

variable . | $Q$ | $T^{\ast}$ |  ..$\tau_{\rm B}^{\ast}$ .. |  .$\chi_{\rm L}$ . | $h$ | $n^{\ast}$ | $f_{\rm ext}^{\ast}$ | $v_{\rm max}$
---|---|---|---|---|---|---|---|---
values . | 20–100 | 0.05–0.4 | 100 | $0$ | 0–5 | 0–0.005 | 0–$5\\!\times\\!10^{-5}$ | 0–0.1
Table 1: Overview of model parameters and their numerical values used in these
simulations: mean number of MPC particles per collision cell $Q=\langle
N_{C_{i}}\rangle$, reduced temperature $T^{\ast}=T/T_{\rm ref}$, dimensionless
Brownian relaxation time $\tau_{\rm B}^{\ast}=\tau_{\rm B}/\Delta t$, the
Langevin susceptibility $\chi_{\rm L}$, the Langevin parameter $h=\mu H/k_{\rm
B}T$, the density ratio $n^{\ast}=nk_{\rm B}T/p_{\rm ref}$, the strength of
the external forcing $f_{\rm ext}^{\ast}$, and the maximum velocity $v_{\rm
max}$.

##  V Results for channel flow

All simulations are started from an equilibrium initial state, where random
velocities are assigned to the particles according to the Maxwell-Boltzmann
distribution at temperature $T$. Initially, particles are placed at random
positions to uniformly fill the simulation box and the orientations
$\mathbf{u}_{i}$ of the particles’ magnetic moments are chosen randomly from
the three-dimensional isotropic distribution. Typically, simulations are run
for $10^{5}$ integration steps and averages are collected only in the second
half of the run where steady-state has been reached.

For $n^{\ast}=0$, backflow effects are absent and the standard translational
dynamics of the MPC fluid is recovered. We tested our numerical implementation
for this case against literature data, dropping angular momentum conservation
and choosing a fixed rotation angle $\alpha$. For equilibrium bulk simulations
with periodic boundary conditions in both spatial dimensions, we recover the
theoretical result for the self-diffusion coefficient,
$D^{\ast}=T^{\ast}\Delta t(1/2+b/[1-b])$ with $b=1/Q+(1-1/Q)\cos\alpha$, where
$Q=\rho a^{2}/m$ is the average number of particles in a collision cell Ihle
and Kroll (2001). For this test, we chose $Q=20$ and varied $\alpha$ between
70 and 120 degrees.

As a second test, we reproduced the channel flow simulations presented in Ref.
Lamura  _et al._ (2001) for $Q=35,\alpha=90$ degrees and $T^{\ast}=0.4$. We
verified that the no-slip condition is satisfied. Moreover, from a fit of the
velocity profile to the theoretical result for Poiseuille flow,

$v_{x}(y)=\frac{f_{\rm ext}}{2\nu}y(L-y),$ (14)

we extract the fluid viscosity $\nu^{\ast}\approx 0.089\pm 0.001$, in
agreement with the findings reported in Ref. Lamura  _et al._ (2001).

Having tested our numerical implementation in the uncoupled case, we now
proceed to investigate the fully coupled system including angular momentum
conservation and backflow effects ($n^{\ast}>0$). For the channel geometry
shown in Fig. 1, we choose the channel width $L=32$ or $L=64$, length $50$,
and $Q=100$ particles per collision cell, i.e. in total $1.6$ or $3.2\times
10^{5}$ particles. For $T^{\ast}=0.1$ and no externally applied field,
$h_{0}=0$, we find a parabolic velocity profile, from which we extract the
kinematic viscosity of the base fluid from Eq. (14) as $\nu_{\rm
s}^{\ast}=0.114\pm 0.001$. We refer to this as the solvent viscosity since the
magnetic force density (2) vanishes in this case. For the same conditions but
in the quiescent state, we extract a self-diffusion coefficient $D\approx
0.05$, in agreement with the $Q\gg 1$ limit of the theoretical expression,
$D^{\ast}\to T^{\ast}\Delta t/2$. Therefore, the Schmidt number measuring the
ratio of viscous over molecular diffusion becomes ${\rm Sc}=\nu_{\rm
s}/D\approx 2.3$. The condition ${\rm Sc}>1$ indicates that collisional
transport dominates over kinetic transport, which is the relevant regime for
fluid simulations Ripoll  _et al._ (2005). Particle-based flow simulations in
general and MPC in particular do not strictly observe the incompressibility
condition $\bm{\nabla}\cdot\mathbf{v}=0$. For small Mach numbers, ${\rm
Ma}=v_{\rm max}/c_{\rm s}\ll 1$, incompressibility is approximately restored.
Here, $v_{\rm max}$ is the maximal fluid velocity and $c_{\rm s}=\sqrt{5k_{\rm
B}T/(3m)}$ the speed of sound. For the present conditions,
$T^{\ast}=0.1,Q=100$, and $f_{\rm ext}=5\times 10^{-5}$, we find $v_{\rm
max}\approx 0.056$, leading to relatively small Mach numbers, ${\rm Ma}\approx
0.14$. Therefore, we can consider the fluid to be approximately
incompressible. The Reynolds number is defined as ${\rm Re}=UL/\nu_{\rm s}$,
where $U$ denotes a characteristic flow velocity. For $U=v_{\rm max}$, we
obtain ${\rm Re}\approx 15.7$, well in the laminar regime for channel flow.
The Weissenberg number ${\rm Wi}=\tau_{\rm B}U/L$ quantifies the ratio of
elastic to viscous forces. Choosing again $v_{\rm max}$ as the characteristic
flow velocity $U$, we find ${\rm Wi}\approx 0.18$, which can be considered to
be still in the Newtonian limit Ilg, Kröger, and Hess (2005).

Having ensured that the simulations are performed in the proper parameter
regime for laminar, viscous, near-incompressible flow, we now investigate the
effect of an external magnetic field on the flow behavior. For simplicity, we
choose $\chi_{\rm L}=0$, in agreement with the model (4), so that
demagnetization effects are absent and the internal field ${\mathbf{h}}$
coincides with the external field ${\mathbf{h}}_{0}$. Other parameters are
chosen as $\tau_{\rm B}^{\ast}=100$ and $n^{\ast}=10^{-3}$.

Figure 2: Velocity profiles for $f^{\ast}_{\rm ext}=5\times
10^{-5},n^{\ast}=10^{-3}$, $\chi_{\rm L}=0$ and various strengths of the
magnetic field $h$. Solid lines show fits to the parabolic profile, Eq. (14).

Figure 2 shows the flow profile for identical external forcing $f^{\ast}_{\rm
ext}=5\times 10^{-5}$ but different strengths of the magnetic field $h$
applied in the gradient direction. We observe from Fig. 2 that the maximum
flow velocity decreases with increasing magnetic field strength. We find that
the Poiseuille profile (14) fits the velocity field very accurately. From
these fits, we extract a field-dependent effective kinematic viscosity
$\nu(h)$. For $h=0$, we obtain the solvent viscosity, $\nu(0)=\nu_{\rm s}$, as
discussed above. For $h>0$, we find that the effective viscosity increases
with increasing magnetic field strength. This phenomenon was pioneered by
McTague McTague (1969) and is known in the literature as “magnetoviscous
effect” Odenbach (2002); Ilg and Odenbach (2008).

Figure 3: Dependence of the relative viscosity $\Delta\nu/\nu(0)$ on the
dimensionless magnetic field strength $h$ for $n^{\ast}=10^{-3}$ and
$n^{\ast}=2\times 10^{-3}$. Other parameters are chosen as $f^{\ast}_{\rm
ext}=5\times 10^{-5}$, $\chi_{\rm L}=0$. Symbols show simulation results,
whereas solid lines correspond to the theoretical result, Eq. (15).

Figure 3 shows the relative viscosity increase $\Delta\nu(h)/\nu(0)$, where
$\Delta\nu(h)=\nu(h)-\nu(0)$ denotes the viscosity change and $\nu(h)$ is
obtained from fits to the velocity profile (14) as described above. We
verified that we obtain identical results within numerical accuracy for weaker
forcing, $f^{\ast}_{\rm ext}=2\times 10^{-5}$. For weak flows, ${\rm Wi}\ll
1$, the model (4) can be solved analytically to give

$\frac{\Delta\nu(h)}{\nu(0)}=\frac{3}{2}\phi\frac{hL_{1}^{2}(h)}{h-L_{1}(h)},$
(15)

where $L_{1}(h)=\coth(h)-1/h$ denotes the Langevin function Martsenyuk,
Raikher, and Shliomis (1974); Ilg, Kröger, and Hess (2002). As can be seen
from Fig 3, the numerical simulations are in very good agreement with this
theoretical result. Also the prefactor $(3/2)\phi$, giving the maximum of the
relative viscosity increase, determined from fits to the numerical data
($\phi\approx 0.27\pm 0.01$ for $n^{\ast}=10^{-3}$) are in good agreement with
the relation $\phi=n^{\ast}\tau_{\rm B}^{\ast}/(3\nu_{s}^{\ast})$ obtained
above (giving here $\phi\approx 0.29$). Before commenting below on the slight
deviation in the value of $\phi$, we want to mention that ferrofluids are
typically dilute, corresponding to smaller values of $\phi$. For the purpose
of demonstrating the method, however, the current choice of parameters should
be sufficient to validate the correct implementation of the hybrid model
combining MPC with the stochastic magnetization dynamics.

Upon closer inspection, we find from the velocity profiles that the no-slip
boundary condition and therefore the parabolic velocity profile is not
perfectly satisfied. This effect is known in the literature and methods have
been proposed to deal with this effect Jonathan K Whitmer (2010).

Figure 4: (a) The profile of the flow vorticity $\Omega_{z}/\Omega_{\rm max}$
is shown across the channel for the same conditions as in Fig. 2. The position
$y$ within the channel is scaled with the channel width $L$ and $\Omega_{\rm
max}$ the maximum vorticity for $h=0$. Solid symbols correspond to $L=32$,
open symbols to $L=64$. Solid lines show a linear fit. (b) The perpendicular
magnetization component $M_{\perp}$ normalized by the saturation magnetization
$M_{\rm sat}$ is shown across the channel with $y$ for the same conditions as
in Fig. 2.

To investigate these deviations in more detail, we show in Fig. 4(a) the
profile of the flow vorticity
$\Omega_{z}=(1/2)[\bm{\nabla}\times\mathbf{v}]_{z}=-(1/2)\mathrm{d}v_{x}(y)/\mathrm{d}y$.
The same conditions as used for Fig. 2 are chosen. We scale $\Omega_{z}$ with
the maximum vorticity expected for Poiseuille flow in the field-free case,
$\Omega_{\rm max}=f_{\rm ext}L/(4\nu_{{\rm s}})$. While the profiles are
nicely linear in the center of the channel as expected for Poiseuille flow,
deviations near the channel walls are apparent. As expected, these deviations
are more pronounced in the vorticity compared to the velocity profile due to
spatial gradients. We performed additional simulations for wider channels
($L=64$) and verified that the perturbation of the profile does not grow with
$L$ but remains a boundary effect, propagating only a distance around $2a$
into the channel.

Figure 4(b) shows the perpendicular magnetization profile $M_{\perp}/M_{\rm
sat}=\langle u_{y}\rangle$ across the channel. As expected, $M_{\perp}$
vanishes in the center of the channel since the vorticity is zero there.
Overall, we find a linear profile of $M_{\perp}$ that follows the linear
profile of the vorticity of the Poiseuille flow shown in Fig. 4(a). However,
deviations from the linear profile are apparent near the channel walls which
result from deviations in the vorticity near the walls observed in Fig. 4(a).
It is worth noting that deviations from the linear magnetization profile occur
over a slightly broader boundary layer compared to the deviations from the
Poiseuille velocity profile.

Discussing these boundary effects in more detail is beyond the scope of the
present work. Here, we just want to mention that the reduced vorticity near
the wall leads to a corresponding reduction of the local perpendicular
magnetization component. As a consequence, the overall viscosity change
$\Delta\nu$ is slightly reduced compared to the theoretical value, as seen in
the slightly reduced value of $\phi$ found above.

##  VI Results for flow around cylinder

As a further demonstration of the flexibility of the MPC method, we here
consider the two-dimensional flow of a ferrofluid around a square cylinder of
diameter $D$ inside a planar channel of width $L$ and length $L_{x}$. No-slip
boundary conditions are imposed on the channel and cylinder walls by the same
bounceback algorithm as used in Sect. V. Following Ref. Lamura  _et al._
(2001), we impose flow by prescribing the average velocity
$v_{x}(y)=\frac{4v_{\rm max}}{L^{2}}y(L-y)$ for particles within the inlet
region $0\leq x<10$. In addition, periodic boundary in the flow direction are
imposed. The horizontal position of the cylinder center is chosen as
$L_{x}/4$, which was found to sufficiently reduce the influence of inflow and
outflow boundary conditions Breuer  _et al._ (2000). Following previous
studies Lamura  _et al._ (2001); Breuer  _et al._ (2000), we choose the
blockage ratio $D/L=1/8$.

Determining the flow around a cylinder is a classical problem in fluid
dynamics where different flow regimes can be distinguished according to the
value of the Reynolds number ${\rm Re}=v_{\rm max}D/\nu$. The creeping flow
regime for ${\rm Re}<1$ is dominated by viscous forces and no separation is
observed. For larger Reynolds numbers, the flow field separates at the
downstream side of the cylinder, forming two steady, counter-rotating vortices
behind the cylinder. The size of this recirculation region $L_{\rm r}$
increases with ${\rm Re}$ until the onset of the van-Kármán vortex street at a
critical Reynolds number Breuer  _et al._ (2000).

We here focus on the regime $5\leq{\rm Re}\leq 30$ where we expect to see
steady recirculation regimes. In order to obtain reasonable spatial
resolution, we choose $L=200,L_{x}=600$. For the average number of particles
per collision cell, we set $Q=30$, which means the simulations contain around
$3.6\times 10^{6}$ particles. The temperature is chosen as $T^{\ast}=0.05$,
corresponding to a viscosity $\nu_{\rm s}^{\ast}\approx 0.0913\pm 0.0001$.
Furthermore, to observe magnetic field-induced changes in the velocity field
more clearly, we extend the range of concentrations to $n^{\ast}=5\times
10^{-3}$. Choosing again $\chi_{\rm L}=0$ eliminates contributions of the
demagnetization field. Simulations are performed for a total of $10^{5}$
integration steps and averages extracted after $5\times 10^{4}$ steps. We
calculate two-dimensional velocity and magnetization fields and their
derivatives using kernel smoothing methods as described in the appendix B.

Figure 5: Stationary velocity field for the flow around a square cylinder
(gray) for $h=0$ (top) and $h=4$ (bottom). For better visibility, only a part
of the velocity field is shown. The simulation parameters are chosen as
$T^{\ast}=0.05$, $n^{\ast}=5\times 10^{-3}$, $\chi_{\rm L}=0$, $v_{\rm
max}=0.1$.

Figure 6: The color-coded magnitude of the stationary velocity field is shown
($h=0$ in the top and $h=4$ in the bottom panel) for the same conditions as in
Fig. 5, but for a larger part of the flow field.

Figure 5 shows the stationary velocity fields with and without applied
magnetic field for a particular choice of model parameters
($T^{\ast}=0.05,n^{\ast}=5\times 10^{-3},v_{\rm max}=0.1$). We observe that
the applied magnetic field changes the flow field and causes a smaller
recirculation region. Qualitatively, such a change is expected due to the
field-induced increase in viscosity. Looking at a larger portion of the
velocity field, Fig. 6 shows the significant changes the external magnetic
field induces not only on the wake, but also on the larger-scale flow field.

Figure 7: (a) Length of the recirculation region $L_{\rm r}$ scaled with the
diameter of the cylinder $D$ as a function of the magnetic field strength $h$.
Diamonds and circles correspond to $n^{\ast}=2\times 10^{-3}$ and $5\times
10^{-3}$, respectively. Other parameters are chosen as $T^{\ast}=0.05,v_{\rm
max}=0.1$. (b) The data in panel (a) are shown as a function of effective
Reynolds number. In addition, open black squares show the result for $h=0$,
$n^{\ast}=5\times 10^{-3}$ but varying $v_{\rm max}$.

To make these observations more quantitative, we extract the length of the
recirculation zone, $L_{\rm r}$, from the velocity fields. In particular, we
determine $L_{\rm r}$ from the zero-crossing of the $x$-component of the
centerline velocity at the end of the wake. We use a parabolic fit to the
centerline velocity profile near the end of the wake to determine $L_{\rm r}$
and estimate error bars based on uncertainties in the fit parameters. Figure 7
shows the length of the recirculation region $L_{\rm r}$ in units of the
diameter $D$ of the cylinder as a function of the magnetic field and the
Reynolds number. Increasing the magnetic field strength leads to a decrease of
$L_{\rm r}$, the effect being more pronounced for larger concentrations. This
observation is consistent with the increased effective viscosity we found in
Sect. V. In the absence of an external magnetic field, by varying the maximum
inflow velocity $v_{\rm max}$, we recover earlier results showing a linear
increase of the length of the recirculation zone with Reynolds number,

$L_{\rm r}/D=c_{0}+c_{1}\ {\rm Re}$ (16)

with $c_{1}\approx 0.059\pm 0.007$, indicated as black dashed line in Fig.
7(b). Keeping instead $v_{\rm max}=0.1$ fixed but varying the magnetic field
strength $h$, we also find the linear relationship Eq. (16) when plotting the
data in Fig. 7(a) against the effective Reynolds number ${\rm Re}=v_{\rm
max}D/\nu(h)$, with $\nu(h)$ the field-enhanced effective kinematic viscosity
discussed in Sect. V. Interestingly, however, we observe that $L_{\rm r}$
varies much less strongly with ${\rm Re}$ due to a magnetic field compared to
the non-magnetic case. Indeed, we find there a slope $c_{1}\approx 0.021\pm
0.003$, less than half the value in the field-free case when varying $v_{\rm
max}$.

To further illustrate the effect of a magnetic field on ferrofluid flow, we
also study the drag coefficient $C_{\rm d}$. The drag coefficient is defined
in two dimensions as $C_{\rm d}=2F_{\|}/(Qv_{\rm max}^{2}D)$, where
$F_{\rm\|}$ denotes the force in flow direction exerted by the fluid on the
cylinder Breuer  _et al._ (2000); Lamura  _et al._ (2001). Within the MPC
scheme, the force $F_{\|}$ can readily be evaluated from the parallel
component of the total momentum the particles transfer to the cylinder wall
when bouncing back. Figure 8 shows the drag coefficient $C_{\rm d}$ as a
function of Reynolds number. In the absence of an externally applied field
(black open symbols), we observe the characteristic strong decrease of $C_{\rm
d}$ with ${\rm Re}$, approximately as an inverse power law, $C_{\rm d}\sim{\rm
Re}^{-z}$ with $z\approx 1.35\pm 0.15$, in agreement with finite-volume
simulations in Ref. Breuer  _et al._ (2000). However, when a magnetic field is
applied (filled colored symbols), we observe that the drag coefficient is
significantly smaller compared to the non-magnetic fluid at the same Reynolds
number. We observe this magnetic drag reduction in the low-Reynolds number
regime $5\lesssim{\rm Re}\lesssim 20$.

Figure 8: The dimensionless drag coefficient $C_{\rm d}$ is shown as a
function of Reynolds number for the same conditions and color coding as in
Fig. 7. Stars indicate the results of Ref. Breuer  _et al._ (2000) and dashed-
line a power-law fit.

##  VII Conclusions

We here present an extension of the MPC method to simulate dynamics and flow
of magnetic fluids, including fluctuation and backflow effects. In analogy to
recent extensions of MPC to nematic liquid crystals Lee and Mazza (2015);
Mandal and Mazza (2019), we equip each MPC particle with a magnetic moment and
include stochastic magnetization dynamics via additional rotational motion of
the individual particles. Fluid and magnetization dynamics are coupled to each
other via velocity gradients and the magnetic force density. We successfully
tested this hybrid scheme in a standard, two-dimensional channel geometry. For
Poiseuille flow, we reproduce the magnetoviscous effect and recover
quantitatively the theoretical result for the relative viscosity increase with
magnetic field strength. Using standard methods to implement no-slip
conditions in MPC simulations, we nevertheless observe small deviations of the
velocity gradient from the theoretical profile very close to the wall, leading
to corresponding deviations in the magnetization profile. This effect is
already known in the literature and more refined methods have been proposed to
better realize the no-slip condition Jonathan K Whitmer (2010). We leave this
more technical aspect for future research. We also illustrate the new MPC
method for magnetic fluids for the benchmark problem of flow around a square
cylinder. In the absence of an applied field, we verify known results for the
length of the recirculation region behind the cylinder increasing linearly
with Reynolds number and the decrease of the drag coefficient for non-magnetic
fluids. In addition, we study the dependence of these quantities on an
externally magnetic field for ferrofluids. We observe that the length of the
wake is increased and the drag is reduced by a magnetic field compared to the
non-magnetic case at the same effective Reynolds number. Being able to reduce
drag by a magnetic field in the low-Reynolds number regime might be of
interest in microfluidics applications. More generally, we find that these
quantities are not described by the effective Reynolds number alone, even when
accounting for the field-induced effective viscosity. This is likely due to
the anisotropy of the effective viscosity induced by the magnetic field
McTague (1969); Sreekumari and Ilg (2015), which leads to ferrofluids showing
different flow behavior than ordinary viscous fluids.

The method proposed here to simulate ferrofluid flow naturally includes
thermal fluctuations that are particularly relevant at small scales.
Furthermore, as a solver for fluctuating ferrohydrodynamics the method
inherits all the benefits of the MPC approach, which can easily be extended to
three spatial dimension and more complicated geometries. The MPC method is
also particularly well-suited as an efficient way to model solvent effects in
colloidal suspension Malevanets and Kapral (2000). Future studies may also
include the effect of demagnetization field that we have neglected here. In
addition, the method is very flexible and can straightforwardly accommodate
different magnetization equations, such as e.g. the so-called chain model
Zubarev and Iskakova (2000) to describe chain-forming ferrofluids or including
mean-field interactions Ilg and Odenbach (2008). Therefore, we expect this
method to be a useful tool for the simulation of flow and hydrodynamic effects
in magnetic fluids.

## Acknowledgments

Discussions with Anoop Varghese during a very early stage of this project are
gratefully acknowledged.

##  Appendix A Stochastic Heun algorithm for orientations

The stochastic Heun algorithm is a predictor-corrector scheme with second-
order weak convergence in the Stratonovich sense Garcia-Palacios (2000). For
rotational motion, we need to ensure that the orientations $\mathbf{u}_{i}(t)$
remain unit vectors for all times.

In the predictor step, new orientations $\breve{\mathbf{u}}_{i}$ are
calculated from an Euler scheme as indicated in Sect. III,

$\displaystyle\breve{\mathbf{u}}_{i}$ $\displaystyle=\frac{\mathbf{u}_{i}^{\rm
P}}{|\mathbf{u}_{i}^{\rm P}|},$ (17) $\displaystyle\mathbf{u}_{i}^{\rm P}$
$\displaystyle=\mathbf{u}_{i}(t)+\Delta\bm{\omega}_{i}(t)\times\mathbf{u}_{i}(t)$
(18)

and $\Delta\bm{\omega}_{i}(t)$ is given by Eq. (9) evaluated at time $t$.

These predicted orientations $\breve{\mathbf{u}}_{i}$ are now used to computed
new angular velocity changes $\Delta\breve{\bm{\omega}}_{i}$, where for the
latter, the right hand side of Eq. (9) is evaluated with $\mathbf{u}_{i}^{\rm
P}$ instead of $\mathbf{u}_{i}(t)$.

In the corrector step, the new orientations are obtained from

$\displaystyle\mathbf{u}_{i}(t+\Delta t)$
$\displaystyle=\frac{\mathbf{u}_{i}^{\rm C}}{|\mathbf{u}_{i}^{\rm C}|},$ (19)
$\displaystyle\mathbf{u}_{i}^{\rm C}$
$\displaystyle=\mathbf{u}_{i}(t)+\frac{1}{2}[\Delta\bm{\omega}_{i}(t)\times\mathbf{u}_{i}(t)+\Delta\breve{\bm{\omega}}_{i}\times\breve{\mathbf{u}}_{i}]$
(20)

##  Appendix B Kernel smoothing

Since evaluation of velocity and magnetostatic fields are important for the
model, we give here some details on their evaluation in the simulation.

We exemplify the method for the velocity component in flow direction, $v_{x}$.
The other fields are evaluated in the same manner. First, we use the Nadaraya-
Watson kernel regression estimator Loader (1999) to find the instantaneous
field as

$v_{x}(\mathbf{r};t)=\frac{\sum_{i=1}^{N}v_{x,i}(t)K\left(\frac{|\mathbf{r}_{i}(t)-\mathbf{r}|}{b}\right)}{\sum_{j=1}^{N}K\left(\frac{|\mathbf{r}_{j}(t)-\mathbf{r}|}{b}\right)}$
(21)

where $K(z)$ is known as the kernel and the parameter $b$ as bandwidth or
smoothing length. While the uniform kernel is frequently used, we here employ
the Epanechnikov kernel

$K(z)=\left\\{\begin{array}[]{cc}\frac{3}{4}(1-z^{2})&|z|\leq 1\\\ 0&{\rm
else}\end{array}\right.$ (22)

which is known to minimize the mean integrated square error. We found that
choosing $b$ equal to the linear size of the collision cells $a$ provides a
good compromise between smoothing and keeping local fluctuations.

Since there is a systematic attenuation bias due to smoothing, we resort to
finite-difference schemes to calculate spatial gradients. Having evaluated the
instantaneous velocity field $v_{x}$ at the centers of the collision cells
from Eq. (21), we use the central difference scheme to approximate the spatial
partial derivatives as

$\frac{\partial v_{x}(\mathbf{r};t)}{\partial
y}\approx\frac{v_{x}(\mathbf{r}+a\hat{\bf k};t)-v_{x}(\mathbf{r}-a\hat{\bf
k};t)}{2a}$ (23)

with the unit vector $\hat{\bf k}$ defined in Fig. 1. And correspondingly with
the unit vector $\hat{\bf i}$ in the $x$-direction for the partial derivative
with respect to $x$. Only near the channel walls we use a first-order
approximation instead. We found that the central difference scheme provides
good results, also in comparison to using higher-order, differentiable kernel
functions.

##  Appendix C Magnetostatics in channel geometry

For non-conducting fluids, Maxwell’s equations are given by Eqs. (3) and a
separate, medium-dependent magnetization equation. First, we consider the
exterior of the system to be non-magnetic,

$\mathbf{M}_{0}=0,\mathbf{B}_{0}=\mu_{0}\mathbf{H}_{0}$ (24)

and we assume the fields $\mathbf{H}_{0}$ (and consequently $\mathbf{B}_{0}$)
are spatially uniform. Thus, Maxwell’s equation (3) are identically satisfied
in the exterior, $\bm{\nabla}\times\mathbf{H}_{0}={\bf 0}$ and
$\bm{\nabla}\cdot\mathbf{B}_{0}=0$.

Let us denote the fields inside the fluid as $\mathbf{H}$ and $\mathbf{B}$
with

$\displaystyle\mathbf{B}$ $\displaystyle=\mu_{0}(\mathbf{H}+\mathbf{M})$ (25)

Define the contribution $\mathbf{H}_{\rm ff}$ of the magnetic fluid to the
internal field,

$\mathbf{H}=\mathbf{H}_{0}+\mathbf{H}_{\rm ff}$ (26)

Note that $\mathbf{H}_{\rm ff}$ is often denoted in terms of a
“demagnetization field”, $\mathbf{H}=\mathbf{H}_{0}-{\bf D}\cdot\mathbf{M}$,
with ${\bf D}$ the demagnetisation tensor, which depends on the shape of the
sample. Note that ${\bf D}$ is symmetric, has trace one, and all diagonal
elements are non-negative Moskowitz and Della Torre (1966).

We now specialize to the channel geometry sketched in Fig. 1. From the
continuity condition on the magnetic field, Eq. (10), we find that

$\displaystyle B_{{\rm}y}$ $\displaystyle=B_{0,y}$ (27) $\displaystyle
H_{{\rm}x}$ $\displaystyle=H_{0,x}\quad\Rightarrow H_{{\rm ff},x}=0$ (28)

Thus, we know that $\mathbf{H}_{\rm ff}=H_{\rm ff}\hat{\bf k}$ is oriented
normal to the wall. Next, to satisfy Maxwell’s equation
$\bm{\nabla}\times\mathbf{H}={\bf 0}$ we require
$\bm{\nabla}\times\mathbf{H}_{\rm ff}={\bf 0}$ and therefore conclude that
$H_{\rm ff}=H_{\rm ff}(y)$ must be independent of the horizontal position $x$
along the channel. Similarly, we can conclude from
$\bm{\nabla}\cdot\mathbf{B}=0$ and the continuity condition that $\partial
B_{{\rm}x}/\partial x=0$ and therefore $\partial M_{x}/\partial x=0$ since
$B_{{\rm}y}$ is constant. Therefore, the magnetization component $M_{x}$ is
also independent of the horizontal position along the channel. Finally, we can
also write $\bm{\nabla}\cdot\mathbf{B}=0$ as
$\bm{\nabla}\cdot(\mathbf{H}+\mathbf{M})=0$, which leads to $\partial(H_{\rm
ff}+M_{y})/\partial y=0$. Since we already concluded that $H_{\rm ff}$ is
independent of $x$, so must $M_{y}$ to satisfy this condition for all points
within the channel. Therefore, none of the magnetostatic fields depends on the
$x$-position within the channel.

## References

  * Rosensweig (1985) R. E. Rosensweig,  _Ferrohydrodynamics_ (Cambridge University Press, Cambridge, 1985).
  * Odenbach (2009) S. Odenbach, ed.,  _Colloidal Magnetic Fluids_ , Lecture Notes in Phys., Vol. 763 (Springer, Berlin, 2009).
  * McTague (1969) J. P. McTague, J. Chem. Phys. 51, 133 (1969).
  * Martsenyuk, Raikher, and Shliomis (1974) M. A. Martsenyuk, Y. L. Raikher, and M. I. Shliomis, Sov. Phys. JETP 38, 413 (1974).
  * Odenbach (2002) S. Odenbach,  _Magnetoviscous Effects in Ferrofluids_ , Lecture Notes in Phys., Vol. 71 (Springer, Berlin, 2002).
  * Kröger, Ilg, and Hess (2003) M. Kröger, P. Ilg, and S. Hess, J. Phys.: Condens. Matter 15, S1403 (2003).
  * Colombo  _et al._ (2012) M. Colombo, S. Carregal-Romero, M. F. Casula, L. Gutiérrez, M. P. Morales, I. B. Böhm, J. T. Heverhagen, D. Prosperi, and W. J. Parak, Chem. Soc. Rev. 41, 4306 (2012).
  * Ilg and Odenbach (2008) P. Ilg and S. Odenbach, in  _Colloidal Magnetic Fluids: Basics, Development and Applications of Ferrofluids_ , Lecture Notes in Phys., Vol. 763, edited by S. Odenbach (Springer, Berlin, 2008).
  * Felicia, Vinod, and Philip (2016) L. J. Felicia, S. Vinod, and J. Philip, Journal of Nanofluids 5, 1 (2016), publisher: American Scientific Publishers.
  * Rinaldi and Zahn (2002) C. Rinaldi and M. Zahn, Phys. Fluids 14, 2847 (2002).
  * Schumacher  _et al._ (2003) K. R. Schumacher, I. Sellien, G. S. Knoke, T. Cader, and B. A. Finlayson, Phys. Rev. E 67, 026308 (2003).
  * Hirabayashi, Chen, and Ohashi (2001) M. Hirabayashi, Y. Chen, and H. Ohashi, Phys. Rev. Lett. 87, 178301 (2001).
  * Sheikholeslami  _et al._ (2018) M. Sheikholeslami, M. Barzegar Gerdroodbary, S. V. Mousavi, D. D. Ganji, and R. Moradi, Journal of Magnetism and Magnetic Materials 460, 302 (2018).
  * Huang, Hädrich, and Michels (2019) L. Huang, T. Hädrich, and D. L. Michels, ACM Trans. Graph. 38, 93 (2019).
  * Kayal  _et al._ (2011) S. Kayal, D. Bandyopadhyay, T. K. Mandal, and R. V. Ramanujan, RSC Advances 1, 238 (2011).
  * Munaz, Shiddiky, and Nguyen (2018) A. Munaz, M. J. A. Shiddiky, and N.-T. Nguyen, Biomicrofluidics 12, 031501 (2018).
  * Malevanets and Kapral (1999) A. Malevanets and R. Kapral, Journal of Chemical Physics 110, 8605 (1999).
  * Ihle and Kroll (2001) T. Ihle and D. M. Kroll, Physical Review E 63, 020201 (2001).
  * Gompper  _et al._ (2009) G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler, Advanced Computer Simulation Approaches For Soft Matter Sciences III 221, 1 (2009).
  * Shendruk and Yeomans (2015) T. N. Shendruk and J. M. Yeomans, Soft Matter 11, 5101 (2015).
  * Lee and Mazza (2015) K.-W. Lee and M. G. Mazza, Journal of Chemical Physics 142, 164110 (2015).
  * Mandal and Mazza (2019) S. Mandal and M. G. Mazza, Physical Review E 99, 063319 (2019).
  * Rosensweig (2002) R. E. Rosensweig, in  _Ferrofluids. Magnetically Controllable Fluids and Their Applications_ , Lecture Notes in Physics No. 594, edited by S. Odenbach (Springer, Berlin, 2002) pp. 61–84.
  * Shliomis (2002) M. I. Shliomis, in  _Ferrofluids. Magnetically Controllable Fluids and Their Applications_ , Lecture Notes in Physics No. 594, edited by S. Odenbach (Springer, Berlin, 2002) pp. 85–111.
  * Leschhorn and Lücke (2006) A. Leschhorn and M. Lücke, Z. Phys. Chem. 220, 219 (2006).
  * Liu and Stierstadt (2008) M. Liu and K. Stierstadt, in  _Colloidal Magnetic Fluids: Basics, Development and Applications of Ferrofluids_ , Lecture Notes in Phys., Vol. 763, edited by S. Odenbach (Springer, Berlin, 2008).
  * Ilg, Kröger, and Hess (2002) P. Ilg, M. Kröger, and S. Hess, J. Chem. Phys. 116, 9078 (2002).
  * Soto-Aquino, Rosso, and Rinaldi (2011) D. Soto-Aquino, D. Rosso, and C. Rinaldi, Physical Review E 84, 056306 (2011), publisher: Cornell University Library.
  * (29) B. Kowalik and R. G. Winkler, Journal of Chemical Physics 138, 104903.
  * Jonathan K Whitmer (2010) J. K. Whitmer and E. Luijten, Journal of Physics Condensed Matter 22, 104106 (2010).
  * Osborn (1945) J. A. Osborn, Phys. Rev. 67, 351 (1945).
  * Lamura  _et al._ (2001) A. Lamura, G. Gompper, T. Ihle, and D. M. Kroll, Europhys. Lett. 56, 319 (2001).
  * Ripoll  _et al._ (2005) M. Ripoll, K. Mussawisade, R. Winkler, and G. Gompper, Phys. Rev. E 72, 016701 (2005).
  * Ilg, Kröger, and Hess (2005) P. Ilg, M. Kröger, and S. Hess, Phys. Rev. E 71, 031205 (2005).
  * Breuer  _et al._ (2000) M. Breuer, J. Bernsdorf, T. Zeiser, and F. Durst, International Journal of Heat and Fluid Flow 21, 186 (2000).
  * Sreekumari and Ilg (2015) A. Sreekumari and P. Ilg, Physical Review E 92, 012306 (2015).
  * Malevanets and Kapral (2000) A. Malevanets and R. Kapral, Journal of Chemical Physics 112, 7260 (2000).
  * Zubarev and Iskakova (2000) A. Y. Zubarev and L. Y. Iskakova, Phys. Rev. E 61, 5415 (2000).
  * Garcia-Palacios (2000) J. L. Garcia-Palacios,  _On the statics and dynamics of magneto-anisotropic nanoparticles_ (John Wiley & Sons, 2000).
  * Loader (1999) C. Loader,  _Local regression and likelihood_ , Statistics and Computing (Springer, 1999).
  * Moskowitz and Della Torre (1966) R. Moskowitz and E. Della Torre, IEEE Transactions on Magnetics 2, 739 (1966).