Modeling hydrodynamic interactions in soft materials with multiparticle collision dynamics

TL;DR

This paper reviews multiparticle collision dynamics (MPCD), a mesoscale simulation method for hydrodynamic interactions in soft materials, highlighting its applications, strengths, limitations, and future prospects.

Contribution

It provides a comprehensive overview of MPCD, summarizing recent applications, discussing current challenges, and outlining emerging developments to enhance its utility in soft matter research.

Findings

MPCD effectively simulates complex fluid flow and bacterial swimming.

Recent applications demonstrate MPCD's versatility in soft material systems.

Emerging developments aim to expand MPCD's applicability to technological systems.

Abstract

Multiparticle collision dynamics (MPCD) is a flexible and robust mesoscale computational technique for simulating solvent-mediated hydrodynamic interactions in soft materials. Here, we provide a critical overview of the MPCD method and summarize its current strengths and limitations. The capabilities of the method are highlighted by reviewing its recent applications to simulate diverse phenomena, ranging from the flow of complex fluids and thermo-osmotic transport to bacterial swimming and active particle self-assembly. We also discuss outstanding challenges and emerging methodological developments that are expected to greatly expand the applicability of MPCD to other systems of technological importance.

Equations

$${\mathbf{r}}_i(t+\Delta t)={\mathbf{r}}_i(t)+{\mathbf{v}}_i(t)\Delta t,$$

$${\mathbf{v}}_i(t+\Delta t)=\mathbf{u}(t)+\mathbf{R}\cdot \left({\mathbf{v}}_i(t)-\mathbf{u}(t)\right),$$

Full text

††journal: Current Opinion in Chemical Engineering

1 Introduction

Solvent-mediated hydrodynamic interactions (HI) play an important role in numerous applications, including self-assembly of active particles, mixing of (complex) liquids, and transport in microfluidic devices. Even at equilibrium, HI can have a dramatic effect on dynamic processes like diffusion. Computer simulations are powerful tools for studying these systems, providing microscopic insights into their static and dynamic properties. Simulations can also be used to explore large parameter spaces and to develop predictive theoretical models used in process design. Nonetheless, simulating soft materials in solution often presents a considerable computational challenge. Many soft materials consist of mesoscopic solutes such as colloids or polymer chains that can be orders of magnitude larger than solvent molecules and have correspondingly slower dynamics. Classical molecular dynamics (MD) simulations using an explicit, atomistically detailed solvent are often intractable for these systems because they require large numbers of atoms and many numerical integration steps to probe the length and time scales of interest.

Often, fine resolution of the solvent is unnecessary, and only the solvent’s effects on the solute need to be described to capture the relevant physics of the system. Such problems are amenable to a mesoscale simulation approach in which the microscopic details of the solvent are simplified by coarse-graining. Various methods have been developed for modeling HI in a coarse-grained fashion, including lattice-based models such as the Lattice-Boltzmann technique [1, 2]; implicit-solvent models like Brownian dynamics [3] and Stokesian dynamics [4]; and particle-based models such as dissipative particle dynamics (DPD) [5, 6] and multiparticle collision dynamics (MPCD) [7]. In this review, we focus on the MPCD method, which has proven to be a useful tool for studying both equilibrium and nonequilibrium behaviors of soft materials. We describe the MPCD algorithm in Section 2, highlight its recent applications to soft materials in Section 3, and then conclude by discussing outstanding challenges and future avenues for development.

2 Fundamentals

2.1 Algorithm

The MPCD method was first proposed by Malevanets and Kapral [7] and is sometimes referred to as stochastic rotation dynamics (SRD) based on their initial variant of the algorithm. The general spirit of MPCD, and other particle-based methods like DPD, is to adopt a computationally inexpensive, coarse grained solvent model that faithfully reproduces the solvent-mediated HI (Fig. 1). The MPCD solvent is a fluid of point particles, each having mass $m$, position $\mathbf{r}$, and velocity $\mathbf{v}$. The particle positions evolve according to Newton’s equations of motion. In the absence of external forces, the positions can be simply stepped forward over a time interval $\Delta t$,

[TABLE]

where the subscript $i$ denotes the particle index. Equation 1 is often referred to as the “streaming” step of the MPCD algorithm (Fig. 1(B)).

After streaming, the solvent undergoes a multiparticle collision that is the core of the MPCD algorithm (Fig. 1(C)). The collision can be performed according to various schemes [7, 10, 11, 12, 13], but all share a few common traits. First, they conserve linear momentum to ensure development of hydrodynamic correlations. Second, they are usually spatially localized and simple (e.g., avoiding pairwise force calculations) for computational efficiency. Finally, they are designed so that the MPCD model reproduces as many relevant transport coefficients and properties of the real solvent as possible, especially the kinematic viscosity and diffusivity, which are intimately connected with the propagation of HI.

In the original algorithm proposed by Malevanets and Kapral [7], commonly referred to as SRD, collisions are handled by sorting MPCD particles into cubic cells with edge length $a$, which sets the spatial resolution of HI [14]. The center-of-mass velocity of each cell, $\mathbf{u}=\sum m_{i}\mathbf{v}_{i}/\sum m_{i}$, is determined. All particles in the same cell then undergo a stochastic collision,

[TABLE]

in which the particle velocities relative to $\mathbf{u}$ are rotated by a norm-conserving rotation matrix $\mathbf{R}$ chosen randomly and independently for each cell. Most commonly, $\mathbf{R}$ rotates the velocities in a cell by a fixed angle $\alpha$ around an axis randomly selected from the unit sphere [10]. Equation 2 conserves both linear momentum and energy within a cell.

In principle, only eqs. 1 and 2 are needed to perform an MPCD simulation. Malevanets and Kapral demonstrated that this algorithm obeys the H theorem, meaning that the particle velocities relax to a Maxwell–Boltzmann distribution, and more importantly, that the algorithm yields the hydrodynamic equations for compressible flow [7]. Nonetheless, additional modifications to the MPCD algorithm are often needed. For instance, discretization of the MPCD particles into collision cells can lead to the loss of Galilean invariance when the particle mean-free path is smaller than the collision cell, resulting in spurious correlations between particle collisions. This issue is typically remedied by randomly shifting the origin of the MPCD cell grid before each collision step to restore Galilean invariance [15].

As another example, eq. 2 conserves linear momentum but not angular momentum [16]. Angular momentum conservation can be restored for eq. 2 [12]; the cost, however, is loss of energy conservation. The MPCD particles must then be coupled to a thermostat to maintain constant temperature, which is also required for nonequilibrium simulations. The most common thermostat for the SRD rule rescales $\mathbf{v}$ relative to $\mathbf{u}$ for each cell by a factor that enforces the correct temperature and ensemble [17]. Alternatively, temperature control can be incorporated directly into the collision rule, as in the Andersen thermostat (AT) scheme [10]. Additional technical details of the basic MPCD algorithm and its many variants can be found in Refs. [18, 19].

2.2 Mapping

One beneficial feature of MPCD is that the solvent transport coefficients can be estimated analytically [20, 21]. For the commonly used SRD parameters of rotation angle $\alpha=130^{\circ}$, collision time $\Delta t=0.1\,\tau$, and particle density $\rho=5\,m/a^{3}$, the analytic expressions predict solvent shear viscosity $\mu=3.96\,k_{\rm B}T\tau/a^{3}$ and self-diffusivity $D=0.064\,a^{2}/\tau$, where $\tau=\sqrt{ma^{2}/k_{\rm B}T}$ is the MPCD unit of time, $k_{\rm B}$ is Boltzmann’s constant, and $T$ is temperature. Analytic estimates of $\mu$ are typically in good agreement with direct calculations from simulation, whereas estimates for $D$ are smaller than the measured values due to hydrodynamic correlations not accounted for in the analytic expressions [22].

It is tempting to map the MPCD solvent model onto a real liquid like water. We set $T=298\,{\rm K}$ and the cell size $a=3\,{\rm\mbox{\AA}}$ to be comparable to the molecular diameter of water. If the MPCD solvent density is matched to that of water ($\rho=997\,{\rm kg}/{\rm m}^{3}$) through $m$, the viscosity of the MPCD solvent is $0.21\,{\rm mPa}\,{\rm s}$, which is lower than the viscosity of water ($0.89\,{\rm mPa}\,{\rm s}$) [23]. Repeating this mapping for hexane, another common solvent used in experiment, which has as a lower density ($\rho=655\,{\rm kg}/{\rm m}^{3}$) and viscosity ($0.30\,{\rm mPa}\,{\rm s}$) than water [23], gives a model viscosity ($0.17\,{\rm mPa}\,{\rm s}$) that is closer to the experimental value. The corresponding MPCD time step under both mappings is roughly 30 fs. This value is significantly larger than the ca. 1 – 2 fs time steps used in atomistic MD simulations, thereby allowing longer effective time scales to be probed.

This mapping has several drawbacks. First, the viscosities are not a perfect match, and the length scale $a$ (and time step) would need to be drastically reduced to obtain better agreement. Second, even taking $a=3\,{\rm\mbox{\AA}}$ may lead to prohibitively large MPCD simulations when the solute is big (tens of nanometers). An alternative mapping based on the solute, which sacrifices properties of the fluid, may be necessary in these cases [21]. Finally, the solvent diffusivity maps to $D\approx 16.8\,\mu{\rm m}^{2}/{\rm ms}$ for water, which is an order of magnitude larger than experimentally measured [24]. This mismatch is inherent to the chosen SRD parameters and mapping, as is apparent from the Schmidt number, ${\rm Sc}=\mu/\rho D\approx 12$. Although this value of Sc is usually considered liquid-like for MPCD [21, 22], it is an order of magnitude smaller than the values for typical solvents such as water. The implications of the MPCD solvent (and other models like DPD and the Lennard-Jones fluid [25]) having low Schmidt number are not fully understood and require further investigation.

2.3 Solute and boundary coupling

A notable advantage of MPCD compared to other techniques for modeling HI in soft materials is its compatibility with a variety of complex solutes and boundaries [25]. The MPCD solvent has been coupled to solutes such as macromolecules [26], rigid colloids [27, 28, 29, 25], and membranes [30]. These systems are typically simulated using a hybrid MD–MPCD scheme in which the solute is propagated using conventional MD techniques while also exchanging momentum with the MPCD solvent during the streaming and/or collision steps.

Macromolecules like polymers are coupled to the MPCD solvent during the collision step [26] as “heavy” particles [26] (Fig. 1(C)). The resulting dynamics of a linear polymer chain are consistent with the Zimm model [31], which includes pre-averaged HI. Slip or no-slip solid boundaries are incorporated during the streaming step using specular reflection (“bounce-back”) rules [32] (Fig. 1(D)) or by randomly redrawing the velocities of reflected particles [28, 29, 33, 25]. Reflection alone is not sufficient to guarantee that there is no slip at a surface [32] because collision cells sliced by the boundary are effectively underfilled with solvent; “virtual” solvent particles must be introduced inside the solid surface during the collision step [32] (Fig. 1(E)). Biasing the velocity of these particles helps enforce the no-slip condition in flow simulations [34].

Possibly the most difficult solute to reliably couple to the MPCD solvent is a moving solid boundary such as a spherical colloid [27, 28, 29, 25] or a deformable membrane [30]. Malevanets and Kapral used short-ranged pair interactions between solvent and solute [27], giving slip boundary conditions and requiring costly MD steps to propagate the solvent. Momentum exchange between the solute and solvent can alternatively be performed using bounce-back rules during streaming [7, 32, 33, 34]. Inoue et al. proposed another scheme [28] where the solvent is stochastically reflected from the solute surface [29, 25]. Virtual particles may be needed to achieve no-slip conditions [32], involving an additional transfer of momentum to the solute [33]. Poblete et al. proposed an alternative methodology where the colloid surface is represented by a penetrable set of point particles connected by springs that are included in the collision step [35].

It remains an open question as to which coupling scheme is optimal for a given problem. Padding and Louis identified spurious depletion forces between colloids when their interactions allowed overlap of volume excluded to the MPCD solvent by the coupling scheme [21]. Bolintineanu et al. found that the early- and long-time diffusion coefficients of colloids in suspension depended sensitively on the coupling scheme and collision rule [25]. Cerbelaud et al. measured the shear viscosity of colloidal suspensions and found that the different coupling schemes were largely equivalent in the dilute regime, but only stochastic bounce-back coupling reproduced the correct dynamics at higher volume fractions [36]. Recent simulations by Shakeri et al. [37] challenge whether MPCD with bounce-back rules is able to produce two-particle HI in the Stokes flow limit, but Poblete et al. [35] were able to reasonably reproduce the Stokes flow around a colloid with their scheme. Further study is required to resolve these outstanding issues and discrepancies. At present, we advocate carefully characterizing the colloid dynamics to identify an appropriate coupling scheme for the problem of interest.

2.4 Computational performance

The MPCD particle density (usually 5 or 10 particles per $a^{3}$) is higher than typical atomistic MD models (roughly 3 atoms per $a^{3}$ for water using the mapping in Section 2.2). Hence, MPCD requires more particles than MD to simulate the same system volume. The comparatively simple MPCD solvent interactions and equations of motions (eq. 1) are less computationally demanding, however, and are readily amenable to parallelization to obtain good performance even with large numbers of particles. Parallel implementations of MPCD have been developed for traditional CPUs [38, 39], but large simulations may require hundreds of CPU cores. Recently, the massive parallelism of graphics processing units (GPUs) has been leveraged to accelerate MPCD [40, 41]. Howard et al.’s open-source implementation of MPCD delivered the performance of nearly 3 CPU nodes (48 cores) using only a single GPU and scaled to run on hundreds of GPUs, enabling simulation volumes as large as $(400\,a)^{3}$ [41]. Moreover, the MPCD simulations of a polymer solution on a GPU were at most a factor of 1.4 slower than a Langevin dynamics simulations without HI, demonstrating the feasibility of routinely performing large-scale MPCD simulations [41]. Massively-parallel, open-source implementations of MPCD can be found in current versions of the HOOMD-blue (2.5.0) and LAMMPS (12 Dec 2018) simulation packages.

3 Applications

3.1 Complex fluids

The MPCD solvent described in Section 2 is a Newtonian fluid, i.e., the shear viscosity $\mu$ is independent of the shear rate. Many biologically and technologically relevant fluids, however, are non-Newtonian complex fluids. For instance, polymer solutions exhibit a distinct viscoelastic response to shear due to the elastic restoring forces counteracting the shear-induced deformation of the polymers. MPCD is well suited to study the dynamics of such complex fluids.

Viscoelasticity has been incorporated in MPCD by bonding pairs of solvent particles with harmonic springs to form dumbbells, resulting in a Maxwell fluid [42]. The MPCD streaming step is straightforward to implement for this model, but it does not reproduce certain nonequilibrium behaviors like shear-thinning due to the possibility of infinite bond extension. Ji et al. proposed to resolve this issue using finitely-extensible bonds [43], but this modification complicates the streaming step. Kowalik and Winkler instead constrained the bond length with a shear-rate-dependent spring constant [44]. Their model has a convenient streaming equation, reproduces the expected shear-thinning behavior, and yields a nonzero second normal stress coefficient.

Alternatively, viscoelasticity can be modeled by dispersing polymers in the MPCD solvent using the scheme in Section 2.3. The equilibrium and nonequilibrium dynamics of flexible polymer chains obtained for dilute [31] and semidilute [45] solutions were in excellent agreement with predictions of the Zimm model. Solutions of semiflexible polymers [46, 47] exhibited a significant increase in shear viscosity (and slowdown of chain dynamics) compared to the flexible polymers. Similar models have been used to study the equilibrium dynamics of other macromolecules, including star polymers [48, 49], ring polymers [50, 51], and microgels [52].

MPCD has also been applied to study the equilibrium dynamics of colloidal suspensions in Newtonian and non-Newtonian fluids. Several studies have investigated the shear viscosity and diffusivity of hard-sphere colloidal suspensions at volume fractions ranging from the dilute limit to $40\,\%$ [25, 36]. Recently, the diffusion of colloidal particles in semidilute polymer solutions was studied in the regime where the colloid diameter is comparable to the polymer mesh size [53, 54, 55] (Fig. 2(A)). These simulations revealed an intricate coupling between the colloid motion and the segmental relaxation and center-of-mass motion of the polymers, leading to subdiffusive dynamics for both species at intermediate timescales.

In addition to equilibrium dynamics, MPCD is well suited to study the nonequilibrium behaviors of complex fluids. Sedimenting hard-sphere colloids were found to exhibit a hydrodynamic (Rayleigh–Taylor-like) instability [61]. Sedimenting red blood cells [57] were shown to assume different shapes depending on the gravitational strength and the elastic moduli of their membranes (Fig. 2(B)), but their terminal velocities were mostly independent of the cell shape. For sedimenting dilute star polymers [58], their shape deformation became more pronounced at increased gravitational strength (Fig. 2(B)).

MPCD has been used extensively to study complex fluids under flow in microfluidic channels, providing insight into several surprising phenomena. Prohm et al. simulated the inertial focusing of spherical colloids in a tube [62] and found that the colloids formed an annulus around the centerline, an analog of the macroscopic Segré-Silberberg effect [63]. The annular distribution was smeared by thermal fluctuations, but sharpened for larger colloids and faster flows. Rigid colloids in dilute polymer solutions instead focused onto the channel center under specific conditions [56] (Fig. 2(A)). Dense colloidal suspensions in pressure-driven flow were found to form periodic velocity and density pulse trains, originating from jamming and release of the colloids [64]. Self-assembled micelles of amphiphilic Janus colloids grew in size and become more symmetric under intermediate shear flows before eventually breaking up under larger shear rates [65].

Recently, the influence of thermal gradients on flow behavior was also investigated with MPCD. Yang and Ripoll showed how microchannels with ratcheted walls can be used as effective microfluidic pumps without moving parts by locally heating the boundaries [59]. The temperature gradient induces thermo-osomotic flows (Fig. 2(C)) that could be exploited to create complex flow patterns [59] or rotate microgears [66]. Burelbach et al. studied the forces acting on a stationary colloid inside a temperature gradient [67]. They found that the thermophoretic forces were proportional to the temperature gradient, in agreement with theoretical predictions, and that the magnitude of this force depended on the hydrodynamic boundary conditions of the colloid. Thermo-osmotic flow from a light-adsorbing colloid led to long-ranged interactions that affected the colloid’s dynamics near a wall [60] (Fig. 2(C)). In addition, the motion of an optically trapped colloid in a light-adsorbing fluid was studied, with the resulting thermo-osmotic flow exerting a torque on the trapped particle.

3.2 Active soft matter

An emerging area of application of MPCD is modeling the hydrodynamics of “active” systems. Unlike the examples above, where solute dynamics arise from equilibrium thermal fluctuations or external perturbations, active particles self-propel through solutions by converting available energy into work. Motile microorganisms such as protozoa and bacteria are the canonical biological example of active particles. They use appendages known as flagella to swim through fluids [68]. Because of their small size and slow swimming speeds, their HI are described by the linear, inertialess Stokes equations. In this regime, swimming motions that are symmetric upon time reversal result in zero net displacement and hence cannot be used to move through the fluid. Consequently, flagellated microorganisms evolved to swim using motions that break time-reversal symmetry [68].

Stark and coworkers [69] used MPCD simulations to investigate the swimming motions of the protozoan parasite Trypanosoma brucei. The trypanosome swims by beating an attached flagellum that is wrapped around its spindle-shaped body in a left-handed helical half-turn (Fig. 3(A)) [69, 70]. The attachment of the flagellum distorts the trypanosome’s flexible cell body into an asymmetric chiral shape, which results in corkscrew-like swimming trajectories that are asymmetric under time reversal. The cell body and flagellum were modeled as an elastic network of small beads connected via harmonic springs. A time-dependent and spatially varying potential was applied along the flagellum to drive its sinusoidal beating motion. The model accurately described the trypanosome’s corkscrew-like forward swimming mode, and it was found that the location of the flagellum along the trypanosome body maximizes its swimming speed. Similar MPCD-based models have been used to elucidate features for swimming of other unicellular flagellates including sperm cells [71] and the bacterium Escherichia coli [72] (Fig. 3(A)).

MPCD simulations have also been used to investigate collective behaviors of active solutes such as clustering and gas–liquid phase separation. In passive systems, these phenomena reflect metastable or equilibrium phases arising from attractive interactions between particles. In active systems, by contrast, they are nonequilibrium steady states that are thought to arise from particle motility, steric repulsions, and HI [76, 74]. Understanding motility-induced clustering (MIC) and motility-induced phase separation (MIPS) in active systems is key to controlling processes such as the formation of the surface-attached bacterial aggregates responsible for biofouling and to develop new routes to self-assemble particles into novel materials [76].

Theers et al. [74] used MPCD to study the effect of particle shape and hydrodynamics on the collective behavior of rigid, prolate spheroidal squirmers in the quasi-two-dimensional confinement of a narrow slit (Fig. 3(B)). Squirmers are generic models for self-propelled solutes. They differ in the type of active stress applied to the surrounding fluid to mimic the flow fields generated by the appendages (flagellum or cilia) of microswimmer classes, including pushers (e.g., Escherichia coli), pullers (e.g., Chlamydomonas), and neutral swimmers (e.g., Paramecium) (Fig. 3(B)). In the simulations, the swimming velocity and active stress were specified by imposing a spatially varying slip velocity across the squirmer surface. Comparison with complementary simulations of non-hydrodynamic active Brownian particles showed that HI suppress and enhance MIPS in systems of neutral spherical and elongated squirmers, respectively, indicating a remarkable sensitivity to swimmer shape (Fig. 3(B)). Further, whereas pullers exhibited an increased tendency to phase separate over neutral swimmers, MIPS was suppressed in systems of pushers. Hence, these findings reveal that collective behavior of active solutes can depend on the complex interplay between HI, particle geometry, and propulsion mechanism.

The collective behavior of self-propelled chemical motors has also been simulated using the MPCD framework. Kapral and coworkers [75] developed a model for Janus motors, which are spherical colloids with asymmetric catalytic activity on their surface (Fig. 3(C)). Reactions occurring on the catalytically active side of the surface generate a local concentration gradient of chemical species and a diffusiophoretic force that propels the particle. The MPCD solvent was modeled as a binary A–B mixture of solvent species that interact with the Janus particle through hard reflective collisions during the MPCD streaming step. The two species were identical, but component A was given a slightly larger collision radius. Particles of species A that collided with the catalytic cap (C) were converted to species B to model the irreversible reaction $\rm{A}+\rm{C}\rightarrow\rm{B}+\rm{C}$. Steady-state conditions were maintained by simulating the reaction $\rm{B}\rightarrow{A}$ in the bulk using a reactive MPCD scheme [77]. Janus motors with small catalytic caps were found to exhibit an increased propensity to cluster at high particle volume fractions (Fig. 3(C)). Particles with large caps, by contrast, tended to exhibit MIC at low densities. This behavior suggests that catalytic cap size and particle volume fraction are key design parameters in controlling the self-assembly and collective dynamics of Janus motors.

4 Conclusions and outlook

Over the past twenty years, MPCD has emerged as a powerful mesoscale technique for simulating HI in soft matter systems such as flowing complex fluids, sedimenting colloidal suspensions, swimming microorganisms, and self-propelling chemical motors. Massively parallel implementations of the algorithm in open-source molecular simulation packages have enabled larger, longer MPCD simulations and lowered the barrier to application of the method. Looking forward, the MPCD method may have promising applications in adaptive resolution schemes where the level of coarse-graining varies in space [78], or in hybrid particle-continuum simulations that utilize microscopic simulations for generating on-the-fly constitutive models for continuum computational fluid dynamics [79].

Despite extensive prior research, there still remain several aspects of the MPCD method in need of further investigation. First, it is challenging to directly map the transport coefficients of the MPCD solvent onto a real liquid (Section 2.2), partially due to the thermal fluctuations in the solvent that are inherent to the collision rule and are hence difficult to tune. The result is a fluid with a high self-diffusivity and low Schmidt number; there is an open question as to whether a solvent with these properties adequately reproduces the relevant HI in the Stokes flow regime [37]. Second, it has been suggested that the anisotropy of the cells underlying the MPCD collision may introduce undesirable artifacts. The recently proposed isotropic SRD collision rule [13] circumvents this issue, but additional work is needed to understand the strengths and limitations of such schemes. Finally, MPCD has had only limited application to multiphase fluids [19, 80], a technologically important class of complex fluids. Future advances in these areas will greatly expand the capabilities of MPCD and solidify its status as one of the most flexible and robust methods for modeling HI in soft materials.

5 Conflicts of interest

The authors declare no conflicts of interest.

6 Acknowledgements

We thank Jacinta Conrad for comments on the manuscript and discussions on bacterial motility. Work by M.P.H. was supported as part of the Center for Materials for Water and Energy Systems, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award #DE-SC0019272. A.N. received support from the German Research Foundation (DFG) under project number NI 1487/2-1. J.C.P gratefully acknowledges support from the Welch Foundation (Grant E-1882) and the National Science Foundation (CBET-1705968).

7 References and recommended reading