---
abstract: 'Over the last few decades, classical density-functional theory (DFT) and its dynamic extensions (DDFTs) have become powerful tools in the study of colloidal fluids. Recently, previous DDFTs for spherically-symmetric particles have been generalised to take into account both inertia and hydrodynamic interactions, two effects which strongly influence non-equilibrium properties. The present work further generalises this framework to systems of anisotropic particles. Starting from the Liouville equation and utilising Zwanzig’s projection-operator techniques, we derive the kinetic equation for the Brownian particle distribution function, and by averaging over all but one particle, a DDFT equation is obtained. Whilst this equation has some similarities with DDFTs for spherically-symmetric colloids, it involves a translational-rotational coupling which affects the diffusivity of the (asymmetric) particles. We further show that, in the overdamped (high friction) limit, the DDFT is considerably simplified and is in agreement with a previous DDFT for colloids with arbitrary shape particles.'
author:
- 'Miguel A.'
- 'Benjamin D.'
- Serafim Kalliadasis
title: Dynamical density functional theory for orientable colloids including inertia and hydrodynamic interactions
---

Introduction {#intro}
============

The study of colloidal fluids, which typically involve particles of micrometer size suspended in a simple fluid bath (where the particles are the atoms-molecules themselves, hence of nanometer size), goes back to the 19th Century with the work of @brown_brief_1828. Since then, numerous studies have been devoted to the description of the dynamics of such systems. A particular challenge is the multiscale nature of the dynamics due to the mass separation of the bath and colloidal particles which in turn implies that the velocities of the bath particles are much higher than those of the colloidal particles. While the time evolution of colloidal fluids can be formally studied by modelling the full systems, i.e. considering the Newtonian dynamics of both bath and colloidal particles, such an approach is computationally intractable, due to both the number of particles and the very different timescales one has to consider[@goddard_general_2012; @goddard_unification_2013], and thus there is a need for coarse graining leading to reduced models. One such reduced model consists of the stochastic Langevin Equations (LEs) for colloidal positions and momenta or equivalently the corresponding Fokker-Planck equation (FPE) for the $N$-body probability distribution function. The LEs were originally proposed heuristically, but can be formally justified as shown in Appendix \[app:projection-operator-tech\] (where for the first time a detailed microscopic derivation of the rotational part of the FPE is offered starting from the full system of bath and colloidal particles). A consequence of this coarse graining down to the dynamics of only the colloidal-particle degrees of freedom is that we no longer deal with deterministic but with random variables. The advantage is that it removes the enormous number of variables related to the dynamical variables of the bath, traditionally via projection-operator techniques which effectively allow to average out the bath[@kirkwood_statistical_1946; @murphy_brownian_1972; @deutch_molecular_1971; @wilemski_derivation_1976; @ermak_brownian_1978; @grabert_microdynamics_1980; @darve_computing_2009]. Despite such a simplification, the resultant equations are still computationally prohibitive because of the large number of particles that need to be taken into consideration to describe the behavior of colloidal fluids. A theoretical way out is to obtain mean field approximations which are independent of the number of particles. For an $N$-particle system, a standard procedure in the statistical mechanics of classical fluids is to average over $N-n$ (with $n < N$) colloidal particles, resulting in a time-evolution equation for the $n$-body reduced distribution function. However, to obtain a closed set of equations for a given $n$ requires knowledge of the relationships between the full and the reduced distributions. Dynamical density functional theories (DDFTs) perform this procedure for $n=1$, i.e. for the one-body density, and typically result in a continuity equation for the density and a time-evolution equation for the current, that is now a functional of the density. The functional that relates the current and the density is generally unknown, but, for a practical implementation of DDFTs it needs to be approximated. A standard approach, for the overdamped dynamics, is to take it as the free-energy functional of a system with the same density at equilibrium. This functional has been well-studied in the statistical mechanics of fluids[@evans_nature_1979; @lutsko_recent_2010]. Also, approximating the functional relating the current and the density with the free energy at equilibrium also ensures that the DDFTs reduce to the corresponding equilibrium DFTs[@goddard_general_2012; @goddard_unification_2013].

The first derivations of a DDFT go back to the early works of @evans_nature_1979 and @dieterich_nonlinear_1990. These initial formulations were mostly phenomenological and subsequent efforts focused on deriving rigorously the same results and also DDFTs that include additional effects such as inertia and hydrodynamic interactions (HIs)[@marconi_dynamic_2000; @archer_dynamical_2004; @archer_dynamical_2009; @espanol_derivation_2009; @goddard_overdamped_2012; @goddard_general_2012; @goddard_unification_2013; @donev_dynamic_2014].

There have also been several efforts to test the derived DDFTs and their predictive capabilities. In fact there have been successful applications to a wide spectrum of physical settings, from hard rod and hard sphere systems[@goddard_unification_2013], mixtures[@goddard_multi-species_2013], thermodynamic phase transitions such as nucleation and spinodal decomposition[@van_teeffelen_colloidal_2008; @lutsko_dynamical_2012] to the calculation of the van Hove distribution function for Brownian hard spheres[@hopkins_van_2010] and crystallization[@neuhaus_density_2014].

DDFTs have been derived for both the overdamped and inertial regimes and although they were initially applied to spherically-symmetric colloids, they have also been generalised to systems where both orientational and translational motions of particles are taken into account, either under weak[@rex_dynamical_2007] or strong[@wittkowski_dynamical_2011] coupling between the two motions. However, these generalisations are restricted to the overdamped regime and neglect HIs. The main goal of the present study is precisely to formulate a DDFT for orientable particles to include both inertia and HIs. Our methodology follows closely the momentum-moments approach adopted by @archer_dynamical_2009 to obtain a DDFT for spherically-symmetric colloidal particles with inertia and also by @goddard_general_2012[@goddard_general_2012; @goddard_unification_2013] to include both inertia and HIs in this DDFT.

In Section \[sec:evolution-equations\] we introduce the equations of motion for a system of $N$ interacting orientable colloidal particles immersed in a thermal bath. Analogous to the spherically-symmetric case, these equations are the LEs for the phase-space coordinates: positions, Eulerian angles, linear and angular (also referred to as rotational) momenta. These equations have already been derived and examined thoroughly before[@dickinson_brownian_1985; @wolynes_dynamical_1977; @hernandez-contreras_brownian_1996], but for the sake of clarity and completeness they are reviewed in Appendices \[app:projection-operator-tech\] and \[app:neglecting-bath-inertia\], in particular highlighting their range of applicability. The derived LEs naturally contain the full rotational-translational coupling which is crucial in the DDFT derivation that follows. In Section \[sec:continuity-current-equation\], we follow[@goddard_general_2012; @goddard_unification_2013] and integrate the Kramer’s equation over all but one particle, yielding a continuity equation for the density distribution and a time-evolution equation for the current. At this point we require closure approximations to obtain a final equation dependent on the density distribution only. These approximations, discussed in detail in Section \[sec:extended-unified-DDFT\], are: (i) the local-equilibrium approximation for the one-body distribution function; (ii) the adiabatic approximation involving the equilibrium free-energy functional; (iii) separation of the $N$-body distribution into local equilibrium and non-equilibrium contributions; and (iv) an approximation for the two-body interactions. We then arrive at the central result of this work, i.e. a generalised DDFT equation to describe systems of orientable particles taking account of both inertia and HIs. Finally, we compare our derived DDFT with previous DDFTs. As expected, if the rotational degrees of freedom are ignored, our DDFT reduces to the one obtained for point-like or spherical particles in the presence of both inertia and HIs[@goddard_general_2012; @goddard_unification_2013], and hence also to several other DDFTs[@marconi_dynamic_2000; @archer_dynamical_2009; @archer_dynamical_2004; @espanol_derivation_2009; @goddard_overdamped_2012] which have neglected inertia or HIs or both. For orientable particles, neglecting inertial effects and HIs allows us to connect our DDFT with previous ones, recovering in particular the DDFTs derived by @rex_dynamical_2007 and @wittkowski_dynamical_2011. Concluding remarks along with a discussion of open problems are offered in Section \[sec:conclusions\].

Evolution equations {#sec:evolution-equations}
===================

Consider a closed system of $N$ identical, asymmetric colloidal particles with mass $m$ immersed in a fluid of $n\gg N$ identical bath particles with mass $m_b \ll m$. This mass separation plays a vital role in deriving the Brownian equations of motion governing colloidal dynamics via projection-operator techniques[@kirkwood_statistical_1946; @murphy_brownian_1972; @deutch_molecular_1971; @wilemski_derivation_1976; @ermak_brownian_1978; @grabert_microdynamics_1980; @darve_computing_2009] (see also Appendix \[app:projection-operator-tech\]). We start by describing the configuration state of the colloidal particles by the position vectors $\mathbf{r}^N=\{\mathbf{r}_1,\dots,\mathbf{r}_N\}$ for their centres of mass and the Eulerian angles ${\boldsymbol{\alpha}}^N=\{{\boldsymbol{\alpha}}_1,\dots,{\boldsymbol{\alpha}}_N\}$, with ${\boldsymbol{\alpha}}_i=(\theta_i,\phi_i,\chi_i)$ defined in the same manner as @condiff_brownian_1966 [@goldstein_classical_2002] or @jose_classical_2013. These angles characterise the orientation of principal-axes frame of the particles, $\mathfrak{B}$, relative to the space-fixed frame, $\mathfrak{S}$. We also denote the translational and angular momenta as ${\mathbf{p}}^N=\{{\mathbf{p}}_1,\dots,{\mathbf{p}}_N\}$ and ${\mathbf{L}}^N=\{{\mathbf{L}}_1,\dots,{\mathbf{L}}_N\}$, respectively. The angular velocity of the $i$th colloidal particle, $\boldsymbol{\omega}_i$, is related to its angular momentum by the equation ${\mathbf{L}}_i=\mathbb{I}_i\boldsymbol{\omega}_i$ with $\mathbb{I}_i$ the inertia tensor. In the principal-axes frames $\mathfrak{B}$, such a tensor is diagonal and constant, $\mathbb{I}_i\equiv\mathbb{I}=\text{diag}(I_{1},I_{2},I_{3})$. The transformation between the frames $\mathfrak{B}$ and $\mathfrak{S}$ is given by the rotation operator $\mathcal{R}_i$,[@jose_classical_2013] such that $\mathbf{x}=\mathcal{R}_i\,\mathbf{x}'$ where $\mathbf{x}\in\mathfrak{B}$ and $\mathbf{x}'\in\mathfrak{S}$. Hence $\boldsymbol{\omega}_i=\mathcal{R}_i\boldsymbol{\omega}_i'$, with $\boldsymbol{\omega}_i'$ the angular velocity in $\mathfrak{S}$, and $\mathbb{I}_i'=\mathcal{R}_i^\top\,\mathbb{I}\,\mathcal{R}_i$, which is neither diagonal nor constant. The angular velocity ${\boldsymbol{\omega}}_i$ can be related to the Eulerian velocities, $\dot{{\boldsymbol{\alpha}}}_i$, via the transformation ${\boldsymbol{\omega}}_i({\boldsymbol{\alpha}}_i) \equiv
\dot{{\boldsymbol{\Phi}}_i}({\boldsymbol{\alpha}}_i) = {\boldsymbol{\Lambda}}_i^\top \dot{{\boldsymbol{\alpha}}}_i$, where the dot denotes time derivative and ${\boldsymbol{\Phi}}_i$ highlights the fact that ${\boldsymbol{\omega}}_i$ is the vector accounting for the rate of change of angular displacement over the Cartesian frame $\mathfrak{B}$. The exact definition of ${\boldsymbol{\Lambda}}^\top$ can be found in Appendix \[app:projection-operator-tech\].

For simplicity, we will consider angular quantities with respect to the body-fixed frames. As such, the dynamics of colloidal particles are described by the Langevin equations[@dickinson_brownian_1985; @dickinson_brownian_1985-1] (see detailed derivation in Appendix \[app:projection-operator-tech\]) $$\begin{aligned}
 \dot{{\mathbf{r}}}_j(t)=&\,\frac{1}{m}{\mathbf{p}}_j,\quad
 \dot{{\boldsymbol{\Phi}}}_j(t)=\,\mathbb{I}^{-1}{\mathbf{L}}_j,\label{eq:1}\\
 \dot{{\mathbf{p}}}_j(t)=&\, -\frac{\partial}{\partial {\mathbf{r}}_j}V({\mathbf{r}}^N,{\boldsymbol{\alpha}}^N)-\sum_{k=1}^N\left({\boldsymbol{\Gamma}}_{jk}^{TT}{\mathbf{p}}_k+{\boldsymbol{\Gamma}}_{jk}^{TR}{\mathbf{L}}_k\right)+\sum_{k=1}^N \left( \mathbf{A}_{jk}^{TT}\mathbf{f}_k(t)+\mathbf{A}_{jk}^{TR}\mathbf{t}_k(t) \right),
 \notag
 \\
 \dot{{\mathbf{L}}}_j(t)=&
 -\frac{\partial}{\partial {\boldsymbol{\Phi}}_j}{V}({\mathbf{r}}^N,{\boldsymbol{\alpha}}^N)-(\mathbb{I}^{-1}{\mathbf{L}}_j)\times{\mathbf{L}}_j
 \,
 -\sum_{k=1}^N\left({\boldsymbol{\Gamma}}_{jk}^{RT}{\mathbf{p}}_k+{\boldsymbol{\Gamma}}_{jk}^{RR}{\mathbf{L}}_k\right)
 +\sum_{k=1}^N \left( \mathbf{A}_{jk}^{RT}\mathbf{f}_k(t)+\mathbf{A}_{jk}^{RR}\mathbf{t}_k(t) \right) \notag,\end{aligned}$$ or, more compactly, $$\begin{aligned}
 \dot{{\boldsymbol{\mathfrak{r}}}}_j(t)=&\ {\boldsymbol{\mathfrak{m}}}_j^{-1}{\boldsymbol{\mathfrak{p}}}_j(t)\label{eq:2},\\
 \dot{{\boldsymbol{\mathfrak{p}}}}_j(t)=&\ {\boldsymbol{\mathfrak{F}}}_j-\sum_{k=1}^N{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N){\boldsymbol{\mathfrak{p}}}_k + \sum_{k=1}^N\mathbf{A}_{jk}\boldsymbol{\xi}_k(t)\notag,\end{aligned}$$ with ${\boldsymbol{\mathfrak{m}}}_j = \text{diag}(m\mathbf{1},\mathbb{I})$, $\boldsymbol{\mathfrak{r}}_j=(\mathbf{r}_j,\boldsymbol{\Phi}_j)$, $\boldsymbol{\mathfrak{p}}_j=(\mathbf{p}_j,\mathbf{L}_j)$ and ${\boldsymbol{\mathfrak{F}}}_j=(-\frac{\partial V}{\partial {\mathbf{r}}_j},-\frac{\partial V}{\partial {\boldsymbol{\Phi}}_j}-(\mathbb{I}^{-1}{\mathbf{L}}_j)\times{\mathbf{L}}_j)$, where $V({\boldsymbol{\mathfrak{r}}}^N)$ is the solvent-averaged interaction-potential of mean force[@snook_langevin_2006]. The gradient operators are defined by $\boldsymbol{\nabla}_{\mathfrak{r}_j}=(\frac{\partial}{\partial
\mathbf{r}_j},\frac{\partial }{\partial \boldsymbol{\Phi}_j})^\top$, with $\frac{\partial}{\partial {\boldsymbol{\Phi}}}=\mathbf{e}_x\frac{\partial}{\partial
\Phi_x}+\mathbf{e}_y\frac{\partial}{\partial
\Phi_y}+\mathbf{e}_z\frac{\partial}{\partial \Phi_z}$ the angular gradient \[see Equations (\[eq:a-17\])–(\[eq:a-18\])\], and $\boldsymbol{\nabla}_{\boldsymbol{\mathfrak{p}}_j}= (\frac{\partial
}{\partial \mathbf{p}_j},\frac{\partial }{\partial \mathbf{L}_j})^\top$. Here ${\boldsymbol{\xi}}_j=(\mathbf{f}_j,\mathbf{t}_j)^\top$ is a 6-dimensional Gaussian white noise vector where $\mathbf{f}_j$ and $\mathbf{t}_j$ are random forces and torques acting upon the $j$-th particle, respectively, such that $\langle
\xi_j^a(t)\rangle=0$ and $\langle \xi_j^a(t)\xi_k^b(t')\rangle =
2\delta_{jk}\delta^{ab}\delta(t-t')$. As for the spherical case, the motion of the colloidal particles induces flows in the bath which in turn result in forces on the colloids, what we have already referred to in Section \[intro\] as HIs. These HIs are represented by the $6\times6$ friction tensors, $
 {\boldsymbol{\Gamma}}_{jk}=({\boldsymbol{\Gamma}}_{jk}^{\mu\nu})
$ where the superscripts $\mu,\nu\in\{T,R\}$ denote the coupling between the translational and angular momenta. Furthermore, the interaction between bath and colloidal particles has a random component representing unpredictable collisions of the bath particles with the colloidal particles. The strength of these random forces is given by the $6\times 6$ tensors $\mathbf{A}_{jk}=(\mathbf{A}_{jk}^{\mu\nu})$ which are coupled to the friction tensors via the fluctuation-dissipation relation,[@wolynes_dynamical_1977; @dickinson_brownian_1985] $$\begin{aligned}
 {\boldsymbol{\mathfrak{m}}}\,{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N) = \beta\sum_{l=1}^N\mathbf{A}_{jl}({\boldsymbol{\mathfrak{r}}}^N)\mathbf{A}_{kl}({\boldsymbol{\mathfrak{r}}}^N)\label{eq:3}.\end{aligned}$$ Here $\beta=1/k_BT$ with $k_B$ the Boltzmann constant and $T$ the temperature, is assumed to be constant. We define the $6N\times6N$ total-mass tensor ${\boldsymbol{\mathfrak{m}}}_{_N}=\text{diag}({\boldsymbol{\mathfrak{m}}}_j)$, along with the $6N\times 6N$ tensors ${\boldsymbol{\Gamma}}=({\boldsymbol{\Gamma}}_{jk})$, and $\mathbf{A}=(\mathbf{A}_{jk})$, transforming Equation (\[eq:3\]) to ${\boldsymbol{\mathfrak{m}}}_{_N}{\boldsymbol{\Gamma}}= \beta\mathbf{A}\mathbf{A}^\top$. Due to the asymmetry of the particles, the friction and noise matrices are not generally symmetric. However, as discussed by @brenner_stokes_1964-2 and @condiff_brownian_1966, the asymmetric part of ${\boldsymbol{\Gamma}}_{ij}$ is exclusively due to nondissipative gyrostatic forces and torques acting on the particles. If such effects can be neglected[@condiff_brownian_1966], then one can assume ${\boldsymbol{\Gamma}}_{jk}^{TR}=({\boldsymbol{\Gamma}}_{jk}^{RT})^\top$. The fluctuation-dissipation relation is then $\mathbf{A} = \sqrt{k_BT {\boldsymbol{\mathfrak{m}}}_{_N}
{\boldsymbol{\Gamma}}}$, analogous to that for spherical particles.

Starting from the Smoluchowski equation, a fluctuating DDFT equation was derived recently[@donev_dynamic_2014]. The authors of this study argued that any equation of motion that accounts for inertial effects of the colloidal particles must include the inertia of the fluid bath. Moreover, they claim that there are only two consistent ways of dealing with this problem: either *fluctuating hydrodynamics*[@hauge_fluctuating_1973] or the overdamped LEs based on the work of @hinch_application_1975 and @roux_brownian_1992. They justify this claim using the estimate obtained by @roux_brownian_1992, who also refers to @hinch_application_1975, of the typical relaxation time of the bath, eventually invalidating the inertial LEs. Nevertheless, while the inertial equations (\[eq:1\]) are neither formally exact nor completely general, it can be shown that there exist very sensible regimes in which neglecting inertia in the bath is appropriate, namely when the colloidal-particle density, $\rho_B$, is much larger than that of the bath particles, $\rho_b$. This is done in Appendix \[app:neglecting-bath-inertia\] using dimensional analysis showing that there is indeed a mass-induced time-scale separation when $\rho_b/\rho_B\ll1$. Under such circumstances the LEs and the associated FPE must be applicable, as explicitly pointed out also by @roux_brownian_1992.

From the system of stochastic differential equations (SDEs) in Eq. (\[eq:2\]), one can obtain[@risken_fokker-planck_1996] (see Appendix \[app:projection-operator-tech\]) the time-evolution equation of the phase-space probability density function (PDF), ${f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t)$, to find each particle $j$ with positions ${\mathbf{r}}_j,\ {\boldsymbol{\alpha}}_j$ and momenta ${\mathbf{p}}_j,\ {\mathbf{L}}_j$ at time $t$, $$\begin{aligned}
\partial_t {f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t) +&
\sum_{j=1}^N{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\cdot \big({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_j{f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t) \big)+\sum_{j=1}^N{\boldsymbol{\mathfrak{F}}}_j\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}{f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t) \label{eq:4}\\
&= \sum_{j,k=1}^N \,{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}\cdot{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N)\left({\boldsymbol{\mathfrak{p}}}_k+k_BT{\boldsymbol{\mathfrak{m}}}\,{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_k}\right){f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t)\notag.\end{aligned}$$ It is well-known that Equations (\[eq:2\]) or (\[eq:4\]) are computationally tractable only for small systems. This problem can be circumvented to an extent by assuming a very simple structure for the friction tensor, namely ${\boldsymbol{\Gamma}}= \text{diag}(\gamma_{TT}\mathbf{1},\gamma_{RR}\mathbf{1})$. While such an approximation certainly reduces the numerical complexity \[since determining $\mathbf{A}$ is now $O(N)$ rather than $O(N^3)$\], it physically corresponds to ignoring HIs, thus neglecting crucial physical information. Another approach to tackle the problem, which is the one we adopt, is to derive a coarse-grained/mean-field model by averaging over the degrees of freedom of all but a few particles[@wu_density-functional_2007; @chan_time-dependent_2005; @goddard_unification_2013; @goddard_multi-species_2013]. This method yields a lower-dimensional system but requires knowledge of the functional relation between ${f^{(N)}}$ and the reduced distribution function. Therefore, convenient and accurate closures for the dynamics of the reduced distribution function must be found.

Continuity and flux equations \[sec:continuity-current-equation\]
=================================================================

Given (\[eq:4\]), we can reduce to a lower dimensional problem by averaging over all but a few particles. Our first assumption involves the solvent-averaged interaction potential $V({\boldsymbol{\mathfrak{r}}}^N)$. Inspired by the decomposition proposed for spherical particles[@goddard_unification_2013] we assume that this potential can be split into linear combinations of $n$-particle interactions, i.e. $$\begin{aligned}
 V({\boldsymbol{\mathfrak{r}}}^N;t)= \sum_{j=1}^NV_1({\boldsymbol{\mathfrak{r}}}_j;t)+\sum_{n=2}^N\frac{1}{n!}\sum_{k_1\neq\dots\neq k_n=1}^n V_n({\boldsymbol{\mathfrak{r}}}_{{k}_1},\dots,{\boldsymbol{\mathfrak{r}}}_{{k}_n})\label{eq:5}.\end{aligned}$$ For ease of notation we also adopt the decomposition ${\boldsymbol{\Gamma}}_{ij}^{\mu\nu}=\gamma_{\mu\nu}\delta_{ij}\mathbf{1}+\gamma_{\mu\nu}\widetilde{{\boldsymbol{\Gamma}}}_{ij}^{\mu\nu}$, recalling that we are not making use of Einstein’s summation convention. Thus, $$\begin{aligned}
 {\boldsymbol{\Gamma}}_{ij}=\delta_{ij}{\boldsymbol{\Upsilon}}+\widetilde{{\boldsymbol{\Upsilon}}}_{ij} ,\end{aligned}$$ with the $6\times 6$ tensor ${\boldsymbol{\Upsilon}}= (\gamma_{\mu\nu}\mathbf{1})$, where $\gamma_{\mu\nu}$ is the constant friction coefficient for an isolated particle and $\widetilde{{\boldsymbol{\Upsilon}}}_{ij}=(\gamma_{\mu\nu}\widetilde{{\boldsymbol{\Gamma}}}_{ij}^{\mu\nu})$, where $\widetilde{{\boldsymbol{\Gamma}}}_{ij}^{\mu\nu}$ are the translational-rotational HI tensors. In the special case when $\gamma_{\mu\nu}=0,\ \forall\mu,\nu$, the FPE (\[eq:4\]) becomes the Liouville equation initially derived by @condiff_brownian_1966[@condiff_brownian_1966; @condiff_transport_1969; @dahler_kinetic_1963] for molecular fluids consisting of non-spherical particles. In this case one usually invokes a Kirkwood-type time smoothing technique[@condiff_brownian_1966; @kirkwood_statistical_1946], rather than a local-equilibrium approximation involving a free-energy functional (see Section \[sec:extended-unified-DDFT\]).

We now define the reduced distributions obtained by integration of the phase-space PDF, $$\begin{aligned}
 {f^{(n)}}({\boldsymbol{\mathfrak{r}}}^n,{\boldsymbol{\mathfrak{p}}}^n;t)=\frac{N!}{(N-n)!} \int d{\boldsymbol{\mathfrak{r}}}^{N-n}d{\boldsymbol{\mathfrak{p}}}^{N-n}\,{f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t)
\label{eq:7}\end{aligned}$$ where $d{\boldsymbol{\mathfrak{r}}}^{N-n}\equiv d{\boldsymbol{\mathfrak{r}}}_{n+1}\dots d{\boldsymbol{\mathfrak{r}}}_N$, with an analogous expression for $d{\boldsymbol{\mathfrak{p}}}$, and the corresponding configuration-space distribution, $$\begin{aligned}
{\rho^{(n)}}({\boldsymbol{\mathfrak{r}}}^n;t) = \int d{\boldsymbol{\mathfrak{p}}}^n\,{f^{(n)}}({\boldsymbol{\mathfrak{r}}}^n,{\boldsymbol{\mathfrak{p}}}^n;t).
 \label{eq:8}\end{aligned}$$ We denote the domains of integration by $\mathfrak{V}=\mathcal{V}\times\mathcal{S}$, with $\mathcal{V}=\mathbb{R}^3$ and $\mathcal{S}=[0,2\pi)\times[0,\pi)\times[0,2\pi)$, and so $\int d{\boldsymbol{\mathfrak{r}}}\equiv \int_{\mathcal{V}}d\mathcal{V}\int_{\mathcal{S}}d\mathcal{S}$ and $\int d{\boldsymbol{\mathfrak{p}}}\equiv \int_{\mathcal{V}}d\mathcal{V}\int_{\mathcal{V}}d\mathcal{V}$, where $$\begin{aligned}
 \int_{\mathcal{V}}d\mathcal{V}=& \int_{0}^\infty dx \int_{0}^\infty dy \int_{0}^\infty dz \label{eq:9}\\
 \int_{\mathcal{S}}d\mathcal{S}=&\int_{0}^{2\pi}d\theta \int_{0}^{\pi}d\phi\, \sin(\phi)\int_{0}^{2\pi} d\chi.
 \notag\end{aligned}$$ For ease of notation $\rho^{(1)}$ is replaced with $\rho$, also known as the one-body density. We now take successive momentum moments of Equation (\[eq:4\]) starting at zeroth order and integrating over all but one-particle position and momentum, i.e. $N \int \text{(\ref{eq:4}) } d{\boldsymbol{\mathfrak{r}}}^{N-1} d{\boldsymbol{\mathfrak{p}}}^{N}$, resulting in the continuity equation $$\begin{aligned}
 \partial_t\rho({\boldsymbol{\mathfrak{r}}}_1;t) + {\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot{\boldsymbol{\mathfrak{j}}}({\boldsymbol{\mathfrak{r}}}_1;t) = 0,
 \label{eq:10}\end{aligned}$$ with the flux defined as ${\boldsymbol{\mathfrak{j}}}({\boldsymbol{\mathfrak{r}}}_1;t)=\int d{\boldsymbol{\mathfrak{p}}}\,({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1)f^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1;t)$.

To obtain the first-momentum moment equation we use $N\int ({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1)\otimes \text{(\ref{eq:4}) }d{\boldsymbol{\mathfrak{r}}}^{N-1}d{\boldsymbol{\mathfrak{p}}}^N$ giving $$\begin{aligned}
 \partial_t{\boldsymbol{\mathfrak{j}}}({\boldsymbol{\mathfrak{r}}}_1;t)&+{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot\int d{\boldsymbol{\mathfrak{p}}}_1({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1\otimes{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1)f^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1;t)\label{eq:11}\\
 &+{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}V_1({\boldsymbol{\mathfrak{r}}}_1;t)\rho({\boldsymbol{\mathfrak{r}}}_1;t)
 +{\boldsymbol{\mathfrak{m}}}^{-1}\sum_{n=2}^N\int d{\boldsymbol{\mathfrak{r}}}_2\dots d{\boldsymbol{\mathfrak{r}}}_n{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}V_n({\boldsymbol{\mathfrak{r}}}^n;t)\rho^{(n)}({\boldsymbol{\mathfrak{r}}}^n;t)\notag\\
 &+{\boldsymbol{\Upsilon}}_{\parallel}{\boldsymbol{\mathfrak{j}}}({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\mathfrak{j}}}_{\bot}({\boldsymbol{\mathfrak{r}}}_1;t)+N{\boldsymbol{\mathfrak{m}}}^{-1}\sum_{j=1}^{N}\int d{\boldsymbol{\mathfrak{r}}}^{N-1}d{\boldsymbol{\mathfrak{p}}}^N\,\widetilde{{\boldsymbol{\Upsilon}}}_{1j}({\boldsymbol{\mathfrak{r}}}^N){\boldsymbol{\mathfrak{p}}}_j{f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t) = 0\notag\end{aligned}$$ where $\otimes$ is the dyadic product. Here we have defined ${\boldsymbol{\Upsilon}}_{\parallel}=\text{diag}(\gamma_{TT}\mathbf{1},\gamma_{RR}\mathbf{1})$, ${\boldsymbol{\Upsilon}}_{\bot}=\text{diag}(\gamma_{TR}\mathbf{1},\gamma_{RT}\mathbf{1})$ and ${\boldsymbol{\mathfrak{j}}}_{\bot}=\int
d{\boldsymbol{\mathfrak{p}}}_1({\boldsymbol{\mathfrak{m}}}^{-1}_{\bot}{\boldsymbol{\mathfrak{p}}}_1)f^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1;t)$, with ${\boldsymbol{\mathfrak{m}}}^{-1}_{\bot}=\text{adiag}(m^{-1}\mathbf{1},\mathbb{I}^{-1})$ the anti-diagonal block matrix with submatrices, $({\boldsymbol{\mathfrak{m}}}_{\bot}^{-1})_{1,2}=
m^{-1}\mathbf{1}$, $({\boldsymbol{\mathfrak{m}}}_\bot^{-1})_{2,1}=\mathbb{I}^{-1}$ and $({\boldsymbol{\mathfrak{m}}}_\bot^{-1})_{i,j}=\mathbf{0}$ otherwise. For particles where coupling cannot be neglected, such as screw-like colloids[@wittkowski_dynamical_2011], this first-momentum moment equation includes an additional term compared to that for spherical particles. This term, involving ${\boldsymbol{\mathfrak{j}}}_\bot({\boldsymbol{\mathfrak{r}}}_1;t)$, reflects the coupled translational and rotational forces. When such a coupling vanishes, as in the case of orthotropic-like particles[@wittkowski_dynamical_2011], the rotational and translational motions are completely decoupled and Equation (\[eq:11\]) can be split into two independent equations.

Equations (\[eq:10\])–(\[eq:11\]), are the basis of our DDFT. In the next section we introduce the closures which link the terms $f^{(1)}$, $\rho^{(n)}$ and ${f^{(N)}}$ of Eq (\[eq:11\]) to the one-body density, following the same reasoning as that applied for the spherical-particle case.

Extending the unified DDFT equation to non-spherical particles\[sec:extended-unified-DDFT\]
===========================================================================================

The closure procedure will be performed in four steps: (i) the term involving the dyadic product will be rewritten using the definition of the kinetic stress, ${\overline{\overline{\boldsymbol{\Pi}}}}({\boldsymbol{\mathfrak{r}}}_1;t)$, and the equilibrium distribution function, $f_{eq}^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1)$; (ii) an adiabatic approximation will be used to substitute the $n$-body density terms with an expression involving the free-energy functional, $\mathcal{F}[\rho]$; (iii) the $N$-body distribution function will be expressed as a sum of a local-equilibrium and non-equilibrium distributions, yielding an extension of the result obtained for spherical-colloidal fluids, and finally; (iv) we restrict our attention to the case of two-body HIs and use an approximation.

Equilibrium distribution, kinetic stress and adiabatic approximation \[subsec:equilibrium-kinetic-stress-adiabatic-approximation\]
----------------------------------------------------------------------------------------------------------------------------------

At equilibrium the one-body distribution function has a Maxwell-Boltzmann structure[@condiff_brownian_1966] $$\begin{aligned}
 f_{eq}^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1)=\frac{\sqrt{|\mathbb{I}^{-1}|}}{m^{3/2}(2\pi k_BT)^3}\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)\exp\left(-\frac{({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1)\cdot{\boldsymbol{\mathfrak{p}}}_1}{2k_BT}\right)
 \label{eq:12},\end{aligned}$$ with $|\mathbb{I}^{-1}|=\prod_{i}(I_{ii})^{-1}$. Hence, one can easily check that $
 \int d{\boldsymbol{\mathfrak{p}}}_1 ({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1\otimes {\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1)f_{eq}^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1)=  k_BT{\boldsymbol{\mathfrak{m}}}^{-1}\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)$. By adding and subtracting ${\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot \left(
k_BT{\boldsymbol{\mathfrak{m}}}^{-1}\rho({\boldsymbol{\mathfrak{r}}}_1)\right)$ on the left-hand side of Equation (\[eq:11\]) we can write, $$\begin{aligned}
 k_BT{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\rho({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot{\overline{\overline{\boldsymbol{\Pi}}}}({\boldsymbol{\mathfrak{r}}}_1;t)&=\,{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot\int d{\boldsymbol{\mathfrak{p}}}_1({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1\otimes{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1)f^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1;t) \label{eq:13}\end{aligned}$$ with the kinetic-stress tensor defined as[@goddard_unification_2013] $$\begin{aligned}
 {\overline{\overline{\boldsymbol{\Pi}}}}({\boldsymbol{\mathfrak{r}}}_1;t)&=\,\int d{\boldsymbol{\mathfrak{p}}}_1\left(({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1\otimes{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_1)-k_BT{\boldsymbol{\mathfrak{m}}}^{-1}\right)f^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1;t)\label{eq:14}.\end{aligned}$$ According to this definition, the term involving the one-body distribution function and the dyadic product in Equation (\[eq:11\]) can be identified as the remainder of the spatial and angular divergence of the kinetic stress. Although at this stage Equation (\[eq:13\]) only seems to be a more convenient way of expressing the second term on the left-hand side of Equation (\[eq:11\]), the advantage will become clear in Section \[subsec:local-equilibrium\] where we discuss how to deal with out-of-equilibrium distributions.

But before this, we will make use of the adiabatic approximation. This approximation identifies, at any time $t$, the one-body density, $\rho({\boldsymbol{\mathfrak{r}}}_1;t)$, with the equilibrium one, $\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)$. We now consider the Helmholtz free-energy functional[@evans_nature_1979; @lutsko_recent_2010] $$\begin{aligned}
 \mathcal{F}[\rho_{eq}]={\textcolor{black}{k_BT}}\int_{\mathfrak{V}}d{\boldsymbol{\mathfrak{r}}}_1\,\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)\left(\ln\left(\Lambda^{3}\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)-1\right)\right)
 +\mathcal{F}_{exc}[\rho_{eq}]+\int_{\mathfrak{V}}d{\boldsymbol{\mathfrak{r}}}_1\,\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)V_1({\boldsymbol{\mathfrak{r}}}_1)
 \label{eq:15}\end{aligned}$$ where the first term on the right-hand side accounts for the ideal-rotator gas component, the last term contains the external potential contribution to the free energy, $\mathcal{F}_{exc}$ is an excess over ideal gas contribution and $\Lambda$ is the de Broglie wavelength (which will turn out to be irrelevant). From classical DFT[@lutsko_recent_2010] we know that there is a unique external potential, $\mathfrak{W}({\boldsymbol{\mathfrak{r}}}_1)=w({\boldsymbol{\mathfrak{r}}}_1)+V_1({\boldsymbol{\mathfrak{r}}}_1)$, for which $\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)$ is the equilibrium density. $${\textcolor{black}{
 \frac{k_BT{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\rho({\boldsymbol{\mathfrak{r}}}_1)}{\rho({\boldsymbol{\mathfrak{r}}}_1)}+{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\mathfrak{W}({\boldsymbol{\mathfrak{r}}}_1) =
 -\left.{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}_{exc}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}\right|_{\rho=\rho_{eq}}}}
 \label{eq:16-a}$$ $$ {\textcolor{black}{{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}w({\boldsymbol{\mathfrak{r}}}_1)=-\left.{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}\right|_{\rho=\rho_{eq}}\equiv
 -{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\mu({\boldsymbol{\mathfrak{r}}}_1)
 }}
 \label{eq:16}$$

Furthermore, under the same circumstances Equation (\[eq:11\]) becomes the first equation of Yvon-Born-Green (YBG) hierarchy[@barrat_basic_2003; @liboff_kinetic_2003] $$k_BT{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)+{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\mathfrak{W}({\boldsymbol{\mathfrak{r}}}_1;t)\rho_{eq}({\boldsymbol{\mathfrak{r}}}_1)=
 {\boldsymbol{\mathfrak{m}}}^{-1}\sum_{n=2}^N\int d{\boldsymbol{\mathfrak{r}}}_2\dots d{\boldsymbol{\mathfrak{r}}}_n{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}V_n({\boldsymbol{\mathfrak{r}}}^n)\rho^{(n)}_{eq}({\boldsymbol{\mathfrak{r}}}^n)
 \label{eq:17}$$ Making the approximation that Equations (\[eq:15\])–(\[eq:17\]) hold out of equilibrium, i.e. substitute $\rho_{eq}^{(n)}$ by $\rho^{(n)}(t)$ for all $n$, yields, $$\begin{aligned}
 \rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}=&\,
 \rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}V_1({\boldsymbol{\mathfrak{r}}}_1;t)+k_BT{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\rho({\boldsymbol{\mathfrak{r}}}_1;t)\label{eq:18}\\
 &+\sum_{n=2}^N\int d{\boldsymbol{\mathfrak{r}}}_2\dots d{\boldsymbol{\mathfrak{r}}}_n{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}V_n({\boldsymbol{\mathfrak{r}}}^n)\rho^{(n)}({\boldsymbol{\mathfrak{r}}}^n;t). \notag\end{aligned}$$ Substituting Equations (\[eq:13\]) and (\[eq:18\]) into Equation (\[eq:11\]) produces, $$\begin{aligned}
 \partial_t{\boldsymbol{\mathfrak{j}}}({\boldsymbol{\mathfrak{r}}}_1;t)+&{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot{\overline{\overline{\boldsymbol{\Pi}}}}({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\mathfrak{m}}}^{-1}\rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}
 +{\boldsymbol{\Upsilon}}_{\parallel}{\boldsymbol{\mathfrak{j}}}({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\mathfrak{j}}}_{\bot}({\boldsymbol{\mathfrak{r}}}_1;t)\label{eq:19}\\
 +&N{\boldsymbol{\mathfrak{m}}}^{-1}\sum_{j=1}^{N}\int d{\boldsymbol{\mathfrak{r}}}^{N-1}d{\boldsymbol{\mathfrak{p}}}^N\,\widetilde{{\boldsymbol{\Upsilon}}}_{1j}({\boldsymbol{\mathfrak{r}}}^N){\boldsymbol{\mathfrak{p}}}_j{f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t) = 0.\notag\end{aligned}$$

Local-equilibrium approximation and beyond\[subsec:local-equilibrium\]
----------------------------------------------------------------------

The next step of the derivation is the introduction of a widely used approximation in statistical physics[@pottier_nonequilibrium_2014; @liboff_kinetic_2003; @goddard_unification_2013], namely the local-equilibrium approximation. It allows rewriting the kinetic stress in terms of the one-body density. With the definitions of local density, $\rho$, velocity, ${\boldsymbol{\mathfrak{v}}}$, and temperature, $T$, given by $$\begin{aligned}
    \int d{\boldsymbol{\mathfrak{p}}}\left\{\begin{array}{c}
             1\\
             {\boldsymbol{\mathfrak{p}}}\\
             ({\boldsymbol{\mathfrak{m}}}^{-1}\mathbf{c})\cdot\mathbf{c}
            \end{array}\right\}
          f^{(1)}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)=
          \left\{\begin{array}{c}
           1\\
           {\boldsymbol{\mathfrak{m}}}{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\\
           6k_BT
          \end{array}\right\}\rho({\boldsymbol{\mathfrak{r}}};t)
          \label{eq:20},\end{aligned}$$ with[@liboff_kinetic_2003] $\mathbf{c}\equiv\mathbf{c}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)={\boldsymbol{\mathfrak{p}}}-{\boldsymbol{\mathfrak{m}}}{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)$ the deviation-from-the-mean microscopic momentum, the local-equilibrium approximation to the true probability distribution, $$\begin{aligned}
 f^{(1)}_{leq}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)=&\, \frac{|\mathbb{I}^{-1}|^{1/2}\rho({\boldsymbol{\mathfrak{r}}};t)}{m^{3/2}(2\pi k_BT)^{3}}\label{eq:21}
  \exp\left(-\frac{({\boldsymbol{\mathfrak{m}}}^{-1}\mathbf{c})\cdot\mathbf{c}}{2k_BT}\right)\end{aligned}$$ fulfils the normalisation conditions given in Equation (\[eq:20\]). While approximating $f^{(1)}$ with Equation (\[eq:21\]) can be reasonable as a first-order approach, there is in fact no actual reason to impose such a restriction. We instead introduce the expansion $f^{(1)}=f^{(1)}_{leq}\,\left[1+\varPsi \right]\equiv
f^{(1)}_{leq}+f^{(1)}_{neq}$, which takes into account both local-equilibrium and non-equilibrium effects. The non-equilibrium term must satisfy $$\begin{aligned}
 \int d{\boldsymbol{\mathfrak{p}}}\left\{\begin{array}{c}
             1\\
             {\boldsymbol{\mathfrak{p}}}\\
             ({\boldsymbol{\mathfrak{m}}}^{-1}\mathbf{c})\cdot\mathbf{c}
            \end{array}\right\}
          f^{(1)}_{neq}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)=
          \left\{\begin{array}{c}
           0\\
           \mathbf{0}\\
           0
          \end{array}\right\}.\label{eq:22}\end{aligned}$$ Hence, $f_{neq}^{(1)}$ only contributes to higher moments such as the stress, ${\overline{\overline{\boldsymbol{\Pi}}}}$, and the heat flux, $\boldsymbol{\mathfrak{Q}}$. Using Equations (\[eq:20\])–(\[eq:22\]) it is readily found that ${\boldsymbol{\mathfrak{j}}}({\boldsymbol{\mathfrak{r}}};t)={\boldsymbol{\mathfrak{j}}}_{leq}({\boldsymbol{\mathfrak{r}}};t)=\rho({\boldsymbol{\mathfrak{r}}};t){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)$. Thus, recalling the definition given in Equation (\[eq:11\]), the kinetic-stress tensor satisfies the decomposition, $$\begin{aligned}
 {\overline{\overline{\boldsymbol{\Pi}}}}({\boldsymbol{\mathfrak{r}}};t)=& ({\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\otimes{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t))\rho({\boldsymbol{\mathfrak{r}}};t)+
 \int d {\boldsymbol{\mathfrak{p}}}\left({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}\otimes {\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}\right)f_{neq}^{(1)}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)\label{eq:23} ={\overline{\overline{\boldsymbol{\Pi}}}}_{leq}({\boldsymbol{\mathfrak{r}}};t) + {\overline{\overline{\boldsymbol{\Pi}}}}_{neq}({\boldsymbol{\mathfrak{r}}};t),\end{aligned}$$ where the first term stands for the local-equilibrium stress, which embodies the kinetic contribution[@condiff_transport_1969] due to translational, rotational and coupled velocities. On the other hand, the heat flux only shows a non-equilibrium contribution and, given the symmetry of $f_{leq}^{(1)}$, $$\begin{aligned}
 \boldsymbol{\mathfrak{Q}}({\boldsymbol{\mathfrak{r}}};t)=\frac{1}{2}\int d{\boldsymbol{\mathfrak{p}}}\,|\mathbf{c}|^2\,{\boldsymbol{\mathfrak{m}}}^{-1}\mathbf{c}\,f^{(1)}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)
=\frac{1}{2}\int d{\boldsymbol{\mathfrak{p}}}\,|\mathbf{c}|^2\,{\boldsymbol{\mathfrak{m}}}^{-1}\mathbf{c}\,f_{neq}^{(1)}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)\label{eq:24}.
 \end{aligned}$$ Now, making use of the identity $$\begin{aligned}
{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}}\cdot\left[({\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\otimes{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t))\rho({\boldsymbol{\mathfrak{r}}};t)\right] =
{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}}\cdot\left({\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\rho({\boldsymbol{\mathfrak{r}}};t)\right)+\rho({\boldsymbol{\mathfrak{r}}};t)\left({\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}}\right){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t),\label{eq:25}\end{aligned}$$ in Equations (\[eq:23\]) and (\[eq:19\]), we arrive at the time-evolution equation $$\begin{aligned}
 \rho({\boldsymbol{\mathfrak{r}}}_1;t)\mathcal{D}_t{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)\label{eq:26}+&{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot {\overline{\overline{\boldsymbol{\Pi}}}}_{neq}({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\mathfrak{m}}}^{-1}\rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}
 +({\boldsymbol{\Upsilon}}_{\parallel}+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\epsilon}})\rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)\\
 +&N{\boldsymbol{\mathfrak{m}}}^{-1}\sum_{j=1}^{N}\int d{\boldsymbol{\mathfrak{r}}}^{N-1}d{\boldsymbol{\mathfrak{p}}}^N\,\widetilde{{\boldsymbol{\Upsilon}}}_{1j}({\boldsymbol{\mathfrak{r}}}^N){\boldsymbol{\mathfrak{p}}}_j{f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t) = 0\notag\end{aligned}$$ where we define the material derivative as $\mathcal{D}_t=(\partial_t+{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}})$ and ${\boldsymbol{\epsilon}}=\text{adiag}(m^{-1}\mathbb{I},m\mathbb{I}^{-1})$. At this stage, there are two terms which still have to be treated, namely the non-equilibrium and the friction tensor terms. The former will be analysed in what follows and the latter in Section \[subsec:enskog-closure\].

To treat the non-equilibrium term, we follow Grad’s method[@grad_kinetic_1949; @liboff_kinetic_2003]. Specifically, the one-body distribution is rewritten as $f^{(1)}=f^{(1)}_{leq}\,\left[1+\varPsi
\right]$ and $\Psi({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)$ is expanded in a basis of orthogonal polynomials of the $\mathbf{p}\oplus\mathbf{L}$ space. In particular, we choose a basis of the multiple generalised Hermite polynomials[@dunkl_classical_2001] on $\mathbb{R}^d$, $$\begin{aligned}
P^{[\alpha]}_{n,\boldsymbol{\ell}}(x_1,\dots,x_d)=H_{\ell_1}^{[\alpha]}(x_1)\dots H_{\ell_d}^{[\alpha]}(x_d)\label{eq:27},
\quad \sum\ell_k=n\end{aligned}$$ with $n\in\mathbb{N}_0=\mathbb{Z}^+\cup\{0\}$, $\boldsymbol{\ell}\in\mathbb{N}^d_0$ and $H_{k}^{[\alpha]}(x)\equiv
H_{k}(x/\sqrt{\alpha})/\sqrt{k!}$ the generalised Hermite polynomials, which are an orthogonal basis for $L^2\left(\mathbb{R}\right)$ with respect to the weight $(2\pi\alpha)^{-1/2} \exp( -x^2/(2\alpha) )$. Denoting $\mathbf{e}_j=(\delta_{jk})$, the first $P^{[\alpha]}_{n,\boldsymbol{\ell}}$ are given by $$\begin{aligned}
 P^{[\alpha]}_{0,\boldsymbol{0}}(\mathbf{x})=\,1,\label{eq:28}\quad P^{[\alpha]}_{1,\mathbf{e}_{j}}(\mathbf{x})=\,\frac{x_j}{\sqrt{\alpha}},\quad
 P^{[\alpha]}_{2,\mathbf{e}_{j}+\mathbf{e}_{k}}(\mathbf{x})=\,\frac{x_jx_k}{\alpha}-\delta_{jk}.\end{aligned}$$ Setting now $\alpha= \sqrt{k_BT}$, the one-body distribution function can be expanded as $$\begin{aligned}
 f^{(1)}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)=f^{(1)}_{leq}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)\left(1+\sum_{n,\boldsymbol{\ell}}B_{n,\boldsymbol{\ell}}({\boldsymbol{\mathfrak{r}}};t)\,
 P^{[\alpha]}_{n,\boldsymbol{\ell}}\left({\boldsymbol{\mathfrak{m}}}^{-1/2}\mathbf{c}\right)\right).\label{eq:29}\end{aligned}$$ with $B_{n,\boldsymbol{\ell}}({\boldsymbol{\mathfrak{r}}};t)$ the moments of the distribution $f_{neq}^{(1)}$. From Equation (\[eq:22\]), we then find that $B_{0,\mathbf{0}}=B_{1,\mathbf{e}_j}=0$ and $\sum_{j}{B}_{2,2\mathbf{e}_j}=0$. Finally, the non-equilibrium term in Equation (\[eq:26\]) can be expressed as $$\begin{aligned}
\boldsymbol{\Pi}_{jk}^{neq}({\boldsymbol{\mathfrak{r}}};t) = k_B T \rho({\boldsymbol{\mathfrak{r}}};t) B_{2,\mathbf{e}_j+\mathbf{e}_k}({\boldsymbol{\mathfrak{r}}};t),
\label{eq:30}\end{aligned}$$ which represents an extension of the analogous relationship found for the spherical case (Section 3.2 of Ref. ). This relationship does not necessarily imply the non-local kinetic pressure to be small, since the functions $B_{2,\mathbf{e}_j+\mathbf{e}_k}({\boldsymbol{\mathfrak{r}}};t)$ do not depend on $\rho,\, T$ or ${\boldsymbol{\mathfrak{v}}}$. However, at this point we make the assumption that $f_{neq}^{(1)}$ can be either approximated or neglected, inspired by previous rigorous analyses for spherical particles[@goddard_overdamped_2012; @goddard_unification_2013]. Therefore, to obtain a DDFT, we need to deal with the last term in Equation (\[eq:26\]).

Two-body interactions and closure equation\[subsec:enskog-closure\]
-------------------------------------------------------------------

With the aim of closing Equation (\[eq:26\]) in $\rho$, we restrict our attention to systems where the HIs are given by linear combinations of two-body interactions. In particular, we assume the general structure[@goddard_unification_2013; @goddard_overdamped_2012] $$\begin{aligned}
 \widetilde{{\boldsymbol{\Upsilon}}}_{jk}({\boldsymbol{\mathfrak{r}}}^N)=\delta_{jk}\sum_{r\neq j}\mathbf{Z}_1({\boldsymbol{\mathfrak{r}}}_j,{\boldsymbol{\mathfrak{r}}}_r)+(1-\delta_{jk})\mathbf{Z}_2({\boldsymbol{\mathfrak{r}}}_j,{\boldsymbol{\mathfrak{r}}}_k)\label{eq:31}.\end{aligned}$$ The integral of Equation (\[eq:26\]) depends only on $f^{(2)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2,{\boldsymbol{\mathfrak{p}}}_1,{\boldsymbol{\mathfrak{p}}}_2;t)$ and we now need a closure relation between the one- and two-body distributions. From a previous rigorous derivation of DDFT[@chan_time-dependent_2005], it is known that the one-body density determines the $N$-body distribution and, therefore, all the time-dependent properties of the system can be expressed as functionals of $\rho({\boldsymbol{\mathfrak{r}}};t)$. Thus, we make the approximation $$\begin{aligned}
f^{(2)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2,{\boldsymbol{\mathfrak{p}}}_1,{\boldsymbol{\mathfrak{p}}}_2;t) = f^{(1)}({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{p}}}_1;t)f^{(1)}({\boldsymbol{\mathfrak{r}}}_2,{\boldsymbol{\mathfrak{p}}}_2;t)g({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2;[\rho]),
\label{eq:32}\end{aligned}$$ although the exact form of the functional $g({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2;[\rho])$ is unknown, accurate approximations are known at equilibrium. However, the problem of obtaining an expression for $g$ is beyond the scope of the present study. From now on, we assume there exists at least a good approximation.

Substituting Equations (\[eq:31\]) and (\[eq:32\]) into (\[eq:30\]) yields $$\begin{aligned}
 \rho({\boldsymbol{\mathfrak{r}}}_1;t)\mathcal{D}_t{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)+&{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot {\overline{\overline{\boldsymbol{\Pi}}}}_{neq}({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\mathfrak{m}}}^{-1}\rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}
 +({\boldsymbol{\Upsilon}}_{\parallel}+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\epsilon}})\rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)\label{eq:33}\\
+&{\boldsymbol{\mathfrak{m}}}^{-1}\rho({\boldsymbol{\mathfrak{r}}}_1;t)\int d{\boldsymbol{\mathfrak{r}}}_2\left[
 \mathbf{Z}_1({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2){\boldsymbol{\mathfrak{m}}}{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)+\mathbf{Z}_2({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2){\boldsymbol{\mathfrak{m}}}{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_2;t)
 \right]\,\rho({\boldsymbol{\mathfrak{r}}}_2;t)\,g({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2;[\rho]) = 0. \notag\end{aligned}$$ Neglecting the non-local equilibrium term as discussed above, and dividing by the one-body density, we finally obtain our equation for the evolution of the flux, $$\begin{aligned}
 \mathcal{D}_t{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)+&{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}\label{eq:34}
 +({\boldsymbol{\Upsilon}}_{\parallel}+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\epsilon}}){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)\\
 +&{\boldsymbol{\mathfrak{m}}}^{-1}\int d{\boldsymbol{\mathfrak{r}}}_2\left[
 \mathbf{Z}_1({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2){\boldsymbol{\mathfrak{m}}}{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)+\mathbf{Z}_2({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2){\boldsymbol{\mathfrak{m}}}{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_2;t)
 \right]\,\rho({\boldsymbol{\mathfrak{r}}}_2;t)\,g({\boldsymbol{\mathfrak{r}}}_1,{\boldsymbol{\mathfrak{r}}}_2;[\rho]) = 0\notag.\end{aligned}$$ Combining this equation with the conservation law resulting from Equations (\[eq:20\]) and (\[eq:22\]) in (\[eq:10\]), namely $$\begin{aligned}
 \partial_t\rho({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot(\rho({\boldsymbol{\mathfrak{r}}};t){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t))=0\label{eq:35},\end{aligned}$$ we finally have the generalised DDFT for systems of orientable particles.

We can readily check that for the special cases of point-like and smooth spherical particles (with perfect slip on the surface[@dickinson_brownian_1985]), these equations reduce to the unified DDFT recently derived by @goddard_unification_2013. Section \[subsec:wl-eq\] will be devoted to showing how Equation (\[eq:34\]) turns into the equation derived by @rex_dynamical_2007 and @wittkowski_dynamical_2011 for systems of arbitrary-shape colloids by setting $\mathbf{Z}_1=\mathbf{Z}_2=0$, and showing that the generalised Einstein relationships[@dickinson_brownian_1985; @peters_smoluchowski_2000] can be recovered in such a case. Finally, we note that the overall effect of the last two terms in Equation (\[eq:34\]) will be a friction-like retardation of the translational and angular velocities ${\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1,t)$, akin to what happened for spherical particles[@goddard_unification_2013].

Connection with previous DDFTs\[subsec:wl-eq\]
----------------------------------------------

Here we connect our DDFT with the DDFTs obtained by @rex_dynamical_2007 and @wittkowski_dynamical_2011 for systems of arbitrary-shape colloids with overdamped dynamics, neglecting HIs. Within the context of our study, this is equivalent to setting $\mathbf{Z}_1=\mathbf{Z}_2=0$. We then have $$\begin{aligned}
  ({\boldsymbol{\Upsilon}}_{\parallel}+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\epsilon}}){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t) = -\mathcal{D}_t{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)-&{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}\label{eq:36}.\end{aligned}$$ Using this result in the continuity equation (\[eq:35\]) gives $$\begin{aligned}
 \partial_t\rho({\boldsymbol{\mathfrak{r}}}_1;t)={\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot\left(\rho({\boldsymbol{\mathfrak{r}}}_1;t)\boldsymbol{\zeta}^{-1}\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}+{\boldsymbol{\mathfrak{m}}}\mathcal{D}_t{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)\right)\right)\label{eq:37}\end{aligned}$$ where we defined $$\begin{aligned}
({\boldsymbol{\Upsilon}}_{\parallel}+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\epsilon}}){\boldsymbol{\mathfrak{m}}}=\boldsymbol{\zeta}\equiv\,\begin{pmatrix}
                                                           \boldsymbol{\zeta}_{\text{TT}} & \boldsymbol{\zeta}_{\text{TR}}\\
                                                           \boldsymbol{\zeta}_{\text{RT}} & \boldsymbol{\zeta}_{\text{RR}}
                                                          \end{pmatrix}\label{eq:38}
 \end{aligned}$$ along with $\boldsymbol{\zeta}_{\text{TT}} =\, m\gamma_{TT}\mathbf{1},\ \boldsymbol{\zeta}_{\text{TR}} =\, m^{-1}\gamma_{TR}\mathbb{I}^2,\ \boldsymbol{\zeta}_{\text{RT}} =\, m^2\gamma_{RT}\mathbb{I}^{-1},$ and $\boldsymbol{\zeta}_{\text{RR}} =\, \gamma_{RR}\mathbb{I}$. In the overdamped regime, inertial forces and unsteady acceleration are neglected, i.e. ${\boldsymbol{\mathfrak{m}}}\,\mathcal{D}_t{\boldsymbol{\mathfrak{v}}}\rightarrow0$. Therefore, Equation (\[eq:38\]) becomes $$\begin{aligned}
\partial_t\rho({\boldsymbol{\mathfrak{r}}}_1;t)={\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot\left(\rho({\boldsymbol{\mathfrak{r}}}_1;t)\beta\mathbf{D}({\boldsymbol{\mathfrak{r}}}_1;t)\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}\right)\right)\label{eq:39}\end{aligned}$$ where $\beta\mathbf{D}=\boldsymbol{\zeta}^{-1}$ is the diffusion tensor. Making use of the block-matrix inversion formula \[e.g. Equation (2.8.18) of Ref. \] we can check that the diffusion tensor fulfils the generalised Einstein relationships[@dickinson_brownian_1985; @peters_smoluchowski_2000] $$\begin{aligned}
  \beta\mathbf{D}_{\text{TT}} =& \left(\boldsymbol{\zeta}_{\text{TT}}-\boldsymbol{\zeta}_{\text{TR}}\boldsymbol{\zeta}_{\text{RR}}^{-1}\boldsymbol{\zeta}_{\text{RT}}\right)^{-1}\label{eq:40}\\
  \beta\mathbf{D}_{\text{TR}} =& -\left(\boldsymbol{\zeta}_{\text{TT}}-\boldsymbol{\zeta}_{\text{TR}}\boldsymbol{\zeta}_{\text{RR}}^{-1}\boldsymbol{\zeta}_{\text{RT}}\right)^{-1}\boldsymbol{\zeta}_{\text{TR}}\boldsymbol{\zeta}_{\text{RR}}^{-1}\equiv
  \left(\boldsymbol{\zeta}_{\text{RT}}-\boldsymbol{\zeta}_{\text{TT}}\boldsymbol{\zeta}_{\text{TR}}^{-1}\boldsymbol{\zeta}_{\text{RR}}\right)^{-1}\notag\\
  \beta\mathbf{D}_{\text{RT}} =& -\left(\boldsymbol{\zeta}_{\text{TT}}-\boldsymbol{\zeta}_{\text{TR}}\boldsymbol{\zeta}_{\text{RR}}^{-1}\boldsymbol{\zeta}_{\text{RT}}\right)^{-1}\boldsymbol{\zeta}_{\text{RT}}\boldsymbol{\zeta}_{\text{TT}}^{-1}
  \equiv\left(\boldsymbol{\zeta}_{\text{TR}}-\boldsymbol{\zeta}_{\text{RR}}\boldsymbol{\zeta}_{\text{RT}}^{-1}
  \boldsymbol{\zeta}_{\text{TT}}\right)^{-1}\notag\\
  \beta\mathbf{D}_{\text{RR}} =& \left(\boldsymbol{\zeta}_{\text{RR}}-\boldsymbol{\zeta}_{\text{RT}}\boldsymbol{\zeta}_{\text{TT}}^{-1}\boldsymbol{\zeta}_{\text{TR}}\right)^{-1}\notag
 \end{aligned}$$ In the special case when the translational and rotational motions are fully decoupled, $\mathbf{D}_{\text{TR}}=\mathbf{D}_{\text{RT}}\equiv 0$, we recover the equation derived by @rex_dynamical_2007. This is indeed the case for spherically isotropic colloids \[e.g. Equation (37) of Ref. \] far from a wall. Furthermore, if the colloidal particles are smooth spheres we recover the DDFT derived by @marconi_dynamic_2000, as the stress tensor is then uniquely determined by the translational component. However, the strength of the hydrodynamic translation-rotation coupling increases rapidly as the distance of the particles from a wall decreases[@dickinson_brownian_1985; @beenakker_many-sphere_1984]. The coupling is then unavoidable near the walls even for the simple spherical case. In this regard, Equation (\[eq:39\]) can be seen as a generalisation of the uncoupled DDFT obtained by @rex_dynamical_2007.

Of course, our overdamped DDFT (\[eq:39\]) does not include the self-propulsion forces present in the DDFT of @wittkowski_dynamical_2011 as these forces were not present in our original equations to begin with (see Appendix \[app:projection-operator-tech\]). Nevertheless, they can readily be included as an additional term, ${\boldsymbol{\mathfrak{F}}}_{i}^{A}({\boldsymbol{\mathfrak{r}}}_i;t)$, inside the solvent-averaged interaction forces and torques ${\boldsymbol{\mathfrak{F}}}$. This results in Equation (\[eq:11\]) having the extra term, $-{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{F}}}_{1}^{A}({\boldsymbol{\mathfrak{r}}}_1;t)\rho({\boldsymbol{\mathfrak{r}}}_1;t)$. Following exactly the same steps as in (\[eq:12\])–(\[eq:16\]), the same first equation of the YBG hierarchy (\[eq:17\]) can be obtained, by setting ${\boldsymbol{\mathfrak{F}}}_1^A=0$ at equilibrium[@wittkowski_dynamical_2011]. Therefore, the additional term $-{\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{F}}}_{1}^{A}({\boldsymbol{\mathfrak{r}}}_1;t)\rho({\boldsymbol{\mathfrak{r}}}_1;t)$ would survive throughout the whole derivation so that the left-hand side of our DDFT equation (\[eq:34\]) would contain it. Eventually, this leads us to $$\begin{aligned}
\partial_t\rho({\boldsymbol{\mathfrak{r}}}_1;t)={\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot\left(\rho({\boldsymbol{\mathfrak{r}}}_1;t)\beta\mathbf{D}({\boldsymbol{\mathfrak{r}}}_1;t)\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}-{\boldsymbol{\mathfrak{F}}}_{1}^{A}({\boldsymbol{\mathfrak{r}}}_1;t)\right)\right)\label{eq:41}\end{aligned}$$ which is now in complete agreement with the equation previously derived within the overdamped regime[@wittkowski_dynamical_2011], with the diffusion tensor still satisfying the generalised Einstein relationships (\[eq:40\]).

Conclusions {#sec:conclusions}
===========

In this work we have formulated a DDFT for orientable particles with inertia and HIs. There are numerous examples where orientation of particles plays a key role. From the fundamental study of perfectly-rough spheres[@condiff_brownian_1966], loaded spherocylinders[@condiff_brownian_1966; @dahler_kinetic_1963], nematic solutions[@hartel_towing_2010; @zhang_isotropic-nematic_2006] or liquid-crystal nucleation[@schilling_self-poisoning_2004; @zhang_isotropic-nematic_2006; @miller_hierarchical_2009], to the study of biological processes such as[@dickinson_brownian_1985] protein adsorption and trapping, antibody-antigen interaction, biochemical assembly by monomer aggregation or polymerization, bone formation[@cantaert_nanoscale_2013; @gomez-morales_progress_2013] or *in vivo* protein crystallization[@bechtel_electron_1976; @koopmann_vivo_2012], to name but a few. We believe that all these problems could be tackled with the aid of the DDFT developed here.

The present study also addresses the problem posed by @wittkowski_dynamical_2011: generalising their work to include HIs. But our study goes a step further and also includes the inertia of the colloidal particles. We are currently in the process of implementing our DDFT equations in the pseudospectral scheme, based on Chebychev collocation[@boyd_chebyshev_2001], we have developed for the numerical treatment of DDFT for spherical particles and which has been successfully applied to a wide spectrum of physical settings; from fluids in a confining potential[@goddard_unification_2013] and mixtures in such a potential[@goddard_multi-species_2013] to adsorbed films on a substrate[@yatsyshin_spectral_2012], fluids in confined geometry[@yatsyshin_geometry-induced_2013; @yatsyshin_wetting_2015; @yatsyshin_density_2015] and contact lines[@nold_fluid_2014; @nold_nanoscale_2015]. While the theory presented here is valid for three-dimensional orientable particles, from a computational point of view systems of two-dimensional orientable particles would be more tractable. Finally, interesting extensions of the framework developed here would be to include torsional degrees of freedom to describe systems of flexible-chain molecules[@evans_momentum_1980], consider multiple-particle species, thus extending the previously developed DDFT for mixtures of spherical particles[@goddard_multi-species_2013], and confined geometry[@yatsyshin_geometry-induced_2013; @yatsyshin_wetting_2015; @yatsyshin_density_2015]. We shall examine these and related questions in future studies.

We are grateful to the anonymous referees for useful comments and suggestions and to Andreas Nold for stimulating discussions that led to the scaling arguments in Appendix B. We acknowledge financial support from the European Research Council via Advanced Grant No. 247031 and from EPSRC via grants No. EP/L020564 and EP/L025159.

Derivation of generalised Langevin equations {#app:projection-operator-tech}
============================================

The Fokker-Planck Equation
--------------------------

Here we outline the derivation of the time-evolution equation for the probability distribution function of a system of arbitrary-shape particles. This derivation can be thought as an application of the works of @kirkwood_statistical_1946, @murphy_brownian_1972, @deutch_molecular_1971 and @wilemski_derivation_1976 in conjunction with the results of @dahler_kinetic_1963, @condiff_brownian_1966, @evans_cumulant_1976 and @condiff_transport_1969, to consider arbitrary-shape particles. .

As in Section \[sec:evolution-equations\], consider $N$ identical, asymmetric colloidal particles with mass $m$ immersed in a fluid of $n\gg N$ identical bath particles with mass $m_b$. Throughout this section upper case letters will refer to colloidal particles while lower case ones indicate bath particles. The phase-space coordinates of the $i$-th bath particle are $\mathbf{x}_j\doteq(\mathbf{r}_j,\,\mathbf{p}_j)$, with “$\doteq$” denoting “by definition", $\mathbf{r}_j$ being the position vector and $\mathbf{p}_j=m_b\dot{\mathbf{r}}_j$ the canonical momentum. In addition, the dynamical state of the $j$-th colloidal particle is determined by $\mathbf{X}_j\doteq(\mathbf{R}_j,\mathbf{P}_j)$, where $\mathbf{R}_j$ is its centre-mass position vector and $\mathbf{P}_j=m\dot{\mathbf{R}}_j$ is its conjugate momentum, along with the pair $\boldsymbol{\Omega}_j\doteq(\boldsymbol{\alpha}_j,\boldsymbol{\pi}_j)$ which comprises the rotational degrees of freedom, with $\boldsymbol{\alpha}_j$ being the Eulerian angles as in Section \[sec:evolution-equations\] and $\boldsymbol{\pi}_j$ their conjugate momenta.[@dahler_kinetic_1963] The angular velocity of a particle, $\boldsymbol{\omega}_j$, is determined by the time derivative of the Euler angles. In the principal-axes frame $\mathfrak{B}$, we have the relation $\boldsymbol{\omega}_j=\boldsymbol{\Lambda}_j^\top\dot{\boldsymbol{\alpha}}_j$, where[@jose_classical_2013; @condiff_brownian_1966] $$\boldsymbol{\Lambda}_j^\top = \begin{pmatrix}
  \cos\chi_j & \sin\theta_j\sin\chi_j & 0\\
  -\sin\chi_j & \sin\theta_j\cos\chi_j & 0\\
  0 & \cos\theta_j & 1
 \end{pmatrix}\ \Leftrightarrow
 \boldsymbol{\Lambda}_j^{-1} =
 \begin{pmatrix}
  \cos\chi_j & \csc\theta_j\cos\chi_j & -\cot\theta_j\sin\chi_j\\
  -\sin\chi_j & \csc\theta_j\sin\chi_j & -\cot\theta_j\cos\chi_j\\
  0 & 0 & 1
 \end{pmatrix}.
 \label{eq:a-1}$$ Accordingly, $\boldsymbol{\omega}_j'=\mathcal{R}_j^\top\boldsymbol{\omega}_j=
\boldsymbol{\Xi}_j\,\dot{\boldsymbol{\alpha}}_j$ under the space-fixed frame $\mathfrak{S}$, with $\boldsymbol{\Xi}_j\doteq\mathcal{R}_j^\top\boldsymbol{\Lambda}_j^\top$. These can be related with their corresponding angular momenta by, $\mathbf{L}_j=\mathbb{I}\,\boldsymbol{\omega}_j$ and $\mathbf{L}'_j=\mathbb{I}'_j\boldsymbol{\omega}'_j$ respectively.

Thus, the dynamical state of the system at any given instant represents a single point in a $6(N+n)-$dimensional space, $\Gamma$,[@liboff_kinetic_2003] $$\mathfrak{s}(t)\doteq(\mathbf{x}_1(t),\dots,\mathbf{x}_n(t),\mathbf{X}_1(t),\dots,\mathbf{X}_N(t),\boldsymbol{\Omega}_1(t),\dots,\boldsymbol{\Omega}_N(t))\equiv
 (\mathbf{x}^n(t),\mathbf{X}^N(t),\boldsymbol{\Omega}^N(t))\in\Gamma
 \label{eq:a-2}$$ where we made use of the notation $\mathbf{x}^n= \mathbf{x}_1\dots
\mathbf{x}_n$, $\mathbf{X}^N=\mathbf{X}_1\dots\mathbf{X}_N$ and $\boldsymbol{\Omega}^N=\boldsymbol{\Omega}_1,\dots,\boldsymbol{\Omega}_N$. From classical mechanics, the evolution of the system is completely determined by the initial conditions for positions and momenta of all particles. This time evolution, which prescribes a unique trajectory $\mathfrak{s}(t;t_0)$, is fully described by the Lagrangian and the Hamiltonian of the system. In the following we use the former to get the relation between $\boldsymbol{\pi}_j$ and $\dot{\boldsymbol{\alpha}}_j$ in order to construct the Hamiltonian. Then Hamilton’s equations along with Liouville’s theorem will be employed to get the time-evolution equation for the PDF of the system so that its phases lie in a differential region of $\Gamma$ with centre placed at $\mathfrak{s}(t)$. First, the Lagrangian of the system can be written as,[^1] $$\begin{aligned}
 \mathcal{L}=&\mathcal{L}_b+\mathcal{L}_B,\label{eq:a-3}\\
 \mathcal{L}_b=&\sum_{i=1}^n\frac{1}{2}m_b\dot{\mathbf{r}}_j\cdot\dot{\mathbf{r}}_j-\mathcal{U}(\mathbf{r}^N,\mathbf{R}^n,\boldsymbol{\alpha}^n),\notag\\
 \mathcal{L}_B=&\sum_{i=1}^N\frac{1}{2}
 \left(m\dot{\mathbf{R}}_j\cdot\dot{\mathbf{R}}_j+\dot{\boldsymbol{\alpha}}_j\cdot(\boldsymbol{\Xi}_j^\top\mathbb{I}'_j\boldsymbol{\Xi}_j) \dot{\boldsymbol{\alpha}}_j\right)-V(\mathbf{R}^N,\boldsymbol{\alpha}^N)\notag,\end{aligned}$$ with $V(\mathbf{R}^N,\boldsymbol{\alpha}^N)$ the potential energy due to short-range interactions exclusively between colloidal particles, and $$\mathcal{U}(\mathbf{r}^n,\mathbf{R}^N,\boldsymbol{\alpha}^N) = U(\mathbf{r}^n)+\sum_{\mu=1}^N\mathfrak{u}_\mu(\mathbf{r}^N,\mathbf{R}_\mu,\boldsymbol{\alpha}_\mu),\label{eq:a-4}$$ the short-range intermolecular potential energy coming from the interaction between bath particles, $U(\mathbf{r}^n)$, and the interaction of each colloidal particle with the whole bath, $\mathfrak{u}_\mu(\mathbf{r}^N,\mathbf{R}_\mu), \forall \mu=1,\dots,N$. From these equations we obtain $$\begin{aligned}
 \mathbf{p}_j=\frac{\partial \mathcal{L}}{\partial \dot{\mathbf{r}}_j}=m_b\dot{\mathbf{r}_j},\quad
 \mathbf{P}_j=\frac{\partial \mathcal{L}}{\partial \dot{\mathbf{R}}_j}=m\dot{\mathbf{R}_j},\quad
 \boldsymbol{\pi}_j=\frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol{\alpha}}_j}= (\boldsymbol{\Xi}_j^\top\mathbb{I}'_j\boldsymbol{\Xi}_j)\dot{\boldsymbol{\alpha}_j},
 \label{eq:a-5}\end{aligned}$$ Therefore, $\boldsymbol{\pi}_j$ can be easily related with the angular momentum $\mathbf{L}_j$ via[@condiff_brownian_1966; @condiff_transport_1969] $$\boldsymbol{\pi}_j = (\boldsymbol{\Xi}_j^\top\mathbb{I}'_j\boldsymbol{\Xi}_j)\dot{\boldsymbol{\alpha}}_j
\equiv
\boldsymbol{\Lambda}_j{\mathbf{L}}_j,
 \label{eq:a-6}$$ which will be used to perform the appropriate change of variables later on. We now write the Hamiltonian function of the system as $$\begin{aligned}
 \mathcal{H} =& \mathcal{H}_b + \mathcal{H}_B \label{eq:a-7}\\
 \mathcal{H}_b =& \sum_{i=1}^n\frac{\mathbf{p}_j\cdot\mathbf{p}_j}{2m_b}+\mathcal{U}(\mathbf{r}^N,\mathbf{R}^n,\boldsymbol{\alpha}^n)\notag\\
 \mathcal{H}_B =&\sum_{i=1}^N\left(\frac{\mathbf{P}_j\cdot\mathbf{P}_j}{2m}+\frac{1}{2}\boldsymbol{\pi}_j\cdot(\boldsymbol{\Xi}_j^\top\mathbb{I}'_j\boldsymbol{\Xi}_j)^{-1}\boldsymbol{\pi}_j
 \right)+V(\mathbf{R}^N,\boldsymbol{\alpha}^n).\notag\end{aligned}$$ According to Liouville’s theorem, the $n+N$ particle distribution function, $\varrho^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N;t)$, will evolve according to $$\partial_t \varrho^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\Omega}^N;t)+\left(i\mathfrak{L}_b+i\mathfrak{L}_B^T+i\mathfrak{L}_B^R\right)\,\varrho^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N;t)=0
 \label{eq:a-8}$$ with the Liouvillian, $\mathfrak{L}\doteq \mathfrak{L}_b+\mathfrak{L}_B^T+\mathfrak{L}_B^R$, $$\begin{aligned}
i\mathfrak{L}_b =& \sum_{j=1}^n\left(\frac{\mathbf{p}_j}{m_b}\cdot\frac{\partial }{\partial \mathbf{r}_j}+
\mathbf{f}_j\cdot\frac{\partial }{\partial \mathbf{p}_j}\right)\label{eq:a-9}\\
i\mathfrak{L}_B^T =& \sum_{j=1}^N\left(\frac{\mathbf{P}_j}{m}\cdot\frac{\partial }{\partial \mathbf{R}_j}+
\mathbf{F}_j\cdot\frac{\partial }{\partial \mathbf{P}_j}\right)\notag\\
i\mathfrak{L}_B^R =& \sum_{j=1}^N\left(\dot{\boldsymbol{\alpha}}_j\cdot\frac{\partial }{\partial \boldsymbol{\alpha}_j}+
\dot{\boldsymbol{\pi}}_j\cdot\frac{\partial }{\partial \boldsymbol{\pi}_j}\right)\notag\end{aligned}$$ where $\mathbf{f}_j=-\frac{\partial}{\partial\mathbf{r}_j}(U+\sum_{\mu=1}^N\mathfrak{u}_\mu)$ and $\mathbf{F}_j=-\frac{\partial}{\partial \mathbf{R}_j}(V+\mathfrak{u}_j)$ are the instantaneous forces acting on bath and colloidal particles, respectively. From a conceptual point of view,[@dahler_kinetic_1963] the quantities ${\boldsymbol{\alpha}}_j$ and ${\boldsymbol{\pi}}_j$ are less convenient than the angular velocities and momenta, $\boldsymbol{\omega}_j$ and $\mathbf{L}_j$. This can be easily checked by considering the rotational kinetic energy $$\begin{aligned}
K_{R}= \sum_{j=1}^N \frac{\dot{\boldsymbol{\alpha}}_j}{2}\cdot(\boldsymbol{\Xi}_j^\top\mathbb{I}'_j\boldsymbol{\Xi}_j) \dot{\boldsymbol{\alpha}}_j \equiv &\ \sum_{j=1}^N \frac{\boldsymbol{\omega}_j}{2}\cdot \mathbb{I}\boldsymbol{\omega}_j\label{a-10}\\
K_{R}= \sum_{j=1}^N \frac{\boldsymbol{\pi}_j}{2}\cdot(\boldsymbol{\Xi}_j^\top\mathbb{I}'_j\boldsymbol{\Xi}_j)^{-1}\boldsymbol{\pi}_j
\equiv&\ \sum_{i=1}^N \frac{\mathbf{L}_j}{2}\cdot(\mathbb{I}^{-1}\mathbf{L}_j).\notag\end{aligned}$$ Transforming now Equation (\[eq:a-8\]) into an equivalent equation for $$\begin{aligned}
 F^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\alpha}^N,\mathbf{L}^N;t)=&\left|
 \frac{\partial(\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\alpha}^N,\boldsymbol{\pi}^N)}{\partial (\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\alpha}^N,\mathbf{L}^N)}
 \right|\varrho^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\Omega}^N;t)\label{eq:a-11}\\
 =&\prod_{j}\sin\theta_j\,\varrho^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\Omega}^N;t)\notag\end{aligned}$$ gives $$\begin{aligned}
 \partial_t F^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\Omega}^N;t)+\left(i\mathfrak{L}_b+i\mathfrak{L}_B^T+i\widetilde{\mathfrak{L}}_B^R\right)\,F^{(n+N)}(\mathbf{x}^n,\mathbf{X}^N,\boldsymbol{\Omega}^N;t)=0
 \label{eq:a-12}\end{aligned}$$ with[@condiff_brownian_1966; @condiff_transport_1969; @dahler_kinetic_1963] $$i\widetilde{\mathfrak{L}}_B^R = \sum_{j=1}^N\dot{\boldsymbol{\alpha}}_j \cdot \hat{\partial}_{\boldsymbol{\alpha}_j}+\mathbf{T}_j\cdot\frac{\partial }{\partial \mathbf{L}_j}
\label{eq:a-13}$$ where the operator[@condiff_brownian_1966] $\hat{\partial}_{\boldsymbol{\alpha_j}}\doteq
(\csc\theta_j\frac{\partial}{\partial \theta_j}\sin\theta_j,
\frac{\partial}{\partial \phi_j}, \frac{\partial}{\partial\chi_j})^\top$, deduced by using the chain rule[@risken_fokker-planck_1996], has been introduced. In the latter equation, $\mathbf{T}_j = \mathbf{N}_j -
\boldsymbol{\omega}_j\times\mathbf{L}_j$ denotes the net torque acting on the $j$-th colloidal particle, with[@condiff_brownian_1966] $$\begin{aligned}
 \mathbf{N}_j &= - \boldsymbol{\Lambda}_j^{-1}\frac{\partial}{\partial \boldsymbol{\alpha}_j}\left[V(\mathbf{R}^N,\boldsymbol{\alpha}^n)+\mathfrak{u}_j(\mathbf{r}^N,\mathbf{R}_j,\boldsymbol{\alpha}_j)\right]
 \label{eq:a-14}\end{aligned}$$ the torque due to intermolecular interactions along the principal axes of inertia. At this point, it is necessary to introduce the angular-gradient operator (also known as orientational gradient[@peters_smoluchowski_2000]) $$\begin{aligned}
 \frac{\partial }{\partial \boldsymbol{\Phi}_j}&\doteq
 \boldsymbol{\Lambda}_j^{-1}\frac{\partial }{\partial \boldsymbol{\alpha}_j}\equiv \hat{\partial}_{\boldsymbol{\alpha}_j}\cdot(\boldsymbol{\Lambda}^\top)^{-1}\label{eq:a-15}\end{aligned}$$ such that $\frac{\partial }{\partial {\boldsymbol{\Phi}}}= \mathbf{e}_x\frac{\partial
}{\partial \Phi_{x}}+\mathbf{e}_y\frac{\partial }{\partial
\Phi_{y}}+\mathbf{e}_z\frac{\partial }{\partial \Phi_{z}}$, where $\mathbf{e}_i$ is the unitary vector along axis $i\in\{x,y,z\}$ of the Cartesian frame $\mathfrak{B}$, and $$\begin{aligned}
 \frac{\partial }{\partial \Phi_{x}^j}&= \cos\chi_j\frac{\partial}{\partial \theta_j}
 +\csc\theta_j\cos\chi_j\frac{\partial }{\partial \phi_j}
 -\cot\theta_j\sin\chi_j\frac{\partial }{\partial \chi_j}\label{eq:a-16}\\
 \frac{\partial }{\partial \Phi_y^j}&=
 -\sin\chi_j\frac{\partial}{\partial \theta_j}
 +\csc\theta_j\sin\chi_j\frac{\partial }{\partial \phi_j}
 -\cot\theta_j\cos\chi_j\frac{\partial }{\partial \chi_j}
 \notag\\
 \frac{\partial }{\partial \Phi_z^j}&=\frac{\partial }{\partial\chi_j}.\notag\end{aligned}$$ It is worth mentioning here that the derivative operators $(\frac{\partial
}{\partial \Phi_{x}},\frac{\partial }{\partial \Phi_{y}},\frac{\partial
}{\partial \Phi_{z}})$ are the generators of rotations of a rigid body about the body-fixed Cartesian frame[@gray_theory_1984]. This results in $$\begin{aligned}
i\widetilde{\mathfrak{L}}_B^R=\sum_{j=1}^N \frac{\partial }{\partial \boldsymbol{\Phi}_j}\cdot\boldsymbol{\omega}_j +  \left(\mathbf{N}_j+\mathbf{L}_j\times\boldsymbol{\omega}_j\right)\cdot\frac{\partial }{\partial\mathbf{L}_j}
\label{eq:a-17}.\end{aligned}$$ Now we can define ${\boldsymbol{\mathfrak{m}}}\doteq \text{diag}(m\mathbf{1},\mathbb{I})$ along with the vectors $\boldsymbol{\mathfrak{r}}_j\doteq
(\mathbf{R}_j,\boldsymbol{\Phi}_j)$, $\boldsymbol{\mathfrak{p}}_j\doteq(\mathbf{P}_j,\mathbf{L}_j)$ and $\boldsymbol{\mathfrak{f}}_j\doteq(\mathbf{F}_j,\mathbf{T}_j)$, and the operators $\boldsymbol{\nabla}_{\mathfrak{r}_j}\doteq(\frac{\partial}{\partial
\mathbf{R}_j},\frac{\partial }{\partial \boldsymbol{\Phi}_j})^\top$ and $\boldsymbol{\nabla}_{\boldsymbol{\mathfrak{p}}_j}\doteq (\frac{\partial
}{\partial \mathbf{P}_j},\frac{\partial }{\partial \mathbf{L}_j})^\top$, enabling us to rewrite Equation (\[eq:a-12\]) in a more compact and convenient way, $$\begin{aligned}
 \partial_t F^{(n+N)}(t)+
 \sum_{j=1}^n\left(\frac{\mathbf{p}_j}{m_b}\cdot\frac{\partial }{\partial \mathbf{r}_j}+
\mathbf{f}_j\cdot\frac{\partial }{\partial \mathbf{p}_j}\right){F^{(n+N)}}(t)
 + \sum_{j=1}^N\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\cdot {\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_j + {\boldsymbol{\mathfrak{f}}}_j\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}
 \right){F^{(n+N)}}(t)=0\label{eq:a-18},\end{aligned}$$ where explicit dependence on the phase-space coordinates was omitted, but recalled through the superscript $(n+N)$. Following @murphy_brownian_1972, we introduce the scaling quantity, ${\boldsymbol{\lambda}}\doteq {\boldsymbol{\mathfrak{m}}}^{-1/2}$, so that $\widetilde{{\boldsymbol{\mathfrak{p}}}}_j={\boldsymbol{\lambda}}{\boldsymbol{\mathfrak{p}}}_j$ and hence, the last term of the previous equation becomes $$\begin{aligned}
 \sum_{j=1}^N{\boldsymbol{\lambda}}\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\cdot {\widetilde{\boldsymbol{\mathfrak{p}}}}_j + {\boldsymbol{\mathfrak{f}}}_j\cdot{\boldsymbol{\nabla}}_{{\widetilde{\boldsymbol{\mathfrak{p}}}}_j}
 \right){F^{(n+N)}}(t). \label{eq:a-19}\end{aligned}$$ Substitution of Equation (\[eq:a-19\]) into Equation (\[eq:a-18\]) the results in an equation resembling Liouville’s equation for spherical colloidal particles. Such a result provides a description of the time evolution of the full system. However, our interest rests exclusively on colloidal particles. Thus, our aim is to get the time-evolution equation for the $N$-particle distribution, $$\begin{aligned}
 {f^{(N)}}(t) \doteq \int d{\mathbf{x}}^n  {F^{(n+N)}}(t).\label{eq:a-20}\end{aligned}$$ For this purpose, Zwanzig’s projection technique can be applied as in the case of spherical colloidal particles[@deutch_molecular_1971; @murphy_brownian_1972; @wilemski_derivation_1976; @ermak_brownian_1978]. Following the work of @murphy_brownian_1972, for the arbitrary initial state at $t=-t_I$ we choose one where the bath particles are in equilibrium with the instantaneous positions of the colloidal particles. This means that $${F^{(n+N)}}({\mathbf{x}}^n,{\mathbf{X}}^N,{\boldsymbol{\Omega}}^N;-t_I)=\rho_n^\dag({\mathbf{x}}^n)\,{f^{(N)}}({\mathbf{X}}^N,{\boldsymbol{\Omega}}^N;-t_I)\label{eq:a-21},$$ with $\rho_n^\dag$ the canonical distribution of the $n$ bath particles in the instantaneous potential created by the colloidal particles. The last step before integrating Equation (\[eq:a-18\]) to remove the dependence on fast variables involves the definition of the projection operator, $$\hat{\mathcal{P}}\doteq \rho_n^\dag({\mathbf{x}}^n)\int d{\mathbf{x}}^n\,,$$ and its complementary operator $\hat{\mathcal{Q}}=1-\hat{\mathcal{P}}$. Thus, the integration of Equation (\[eq:a-18\]) over the fast variables ${\mathbf{x}}^n$ is equivalent to applying $(\rho_n^\dag)^{-1}\hat{\mathcal{P}}$ on both sides of such an equation. Although considerable algebraic manipulations are required, we can follow the work of @murphy_brownian_1972 and @lebowitz_microscopic_1965 step by step to finally reach the desired time-evolution equation for the projected distribution function ${f^{(N)}}$, $$\begin{aligned}
 \partial_t{f^{(N)}}(t)+\sum_{j=1}^N&{\boldsymbol{\lambda}}\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\cdot {\widetilde{\boldsymbol{\mathfrak{p}}}}_j+{\boldsymbol{\mathfrak{F}}}_j\cdot{\boldsymbol{\nabla}}_{{\widetilde{\boldsymbol{\mathfrak{p}}}}_j}\right){f^{(N)}}(t) \label{eq:a-23}\\
 =& \sum_{j,k=1}^N {\boldsymbol{\lambda}}^2\,{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}\cdot\int_{-t_I}^{t}dt'{\boldsymbol{\gamma}}_{jk}({\boldsymbol{\lambda}};t,t')\left({\boldsymbol{\nabla}}_{{\widetilde{\boldsymbol{\mathfrak{p}}}}_k}+\beta\,{\widetilde{\boldsymbol{\mathfrak{p}}}}_k\right){f^{(N)}}(t')\notag\end{aligned}$$ where $\beta=1/k_BT$, where $k_B$ is the Boltzmann constant, $T$ is the temperature imposed by the bath, and ${\boldsymbol{\mathfrak{F}}}$ is the equilibrium average force and torque, i.e. $$\begin{aligned}
  {\boldsymbol{\mathfrak{F}}}_j &\equiv
   \left\langle
  \begin{array}{l}
  {\mathbf{F}}_j\\
  {\mathbf{T}}_j
 \end{array}
  \right\rangle^\dag
  =\int d{\mathbf{x}}^n \rho_n^\dag({\mathbf{x}}^n)\,{\boldsymbol{\mathfrak{f}}}_j({\mathbf{r}}^n,{\mathbf{R}}^N,{\boldsymbol{\alpha}}^N)
 \label{eq:a-24}\\
 &=-\left[
 \begin{array}{l}
  -\frac{\partial }{\partial {\mathbf{R}}_j} V({\mathbf{R}}^N,{\boldsymbol{\alpha}}^N)\\
  -\frac{\partial }{\partial {\boldsymbol{\Phi}}_j}V({\mathbf{R}}^N,{\boldsymbol{\alpha}}^N)- {\boldsymbol{\omega}}_j\times{\mathbf{L}}_j
 \end{array}
 \right]+\left\langle{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\mathfrak{u}({\mathbf{r}}^n,{\mathbf{R}}_j,{\boldsymbol{\alpha}}_j)\right\rangle^\dag =
 \left[
 \begin{array}{l}
  -\frac{\partial }{\partial {\mathbf{R}}_j} (V+\psi)\\
  -\frac{\partial }{\partial {\boldsymbol{\Phi}}_j}(V+\psi)- {\boldsymbol{\omega}}_j\times{\mathbf{L}}_j
 \end{array} \right],
 \notag\end{aligned}$$ with $\langle.\rangle^\dag$ the equilibrium average over the fast variables. Equation (\[eq:a-24\]) also includes the definition of the potential of mean force $\psi\doteq
\langle\mathfrak{u}\rangle^\dag$, i.e. the potential which gives rise to the average (over all configurations of the $n$ bath molecules) force and torque acting on the $j$th colloidal particle at any given configuration keeping all the colloidal particles frozen. Thus, the fluid-equilibrium average force and torque includes the contribution of a postulated vacuum colloid-colloid interaction potential, $V$, and the solvent contribution to the total force and torque, $\psi$. This combination in turn results in a solvent-averaged potential of mean force, $\widetilde{V}\doteq V+\psi$, which can be obtained from a given physical model, e.g. the DLVO theory for the interaction of charged colloidal particles[@snook_langevin_2006] or the ten Wolde-Frenkel potential[@wolde_enhancement_1997]. Finally, the tensor ${\boldsymbol{\gamma}}_{jk}$ is given by $$\begin{aligned}
 {\boldsymbol{\gamma}}_{jk}({\boldsymbol{\lambda}};t,t')&=\left\langle {\boldsymbol{\mathfrak{f}}}_j(t)\otimes e^{i(t'-t)\hat{\mathcal{Q}}\mathfrak{L}}({\boldsymbol{\mathfrak{f}}}_k(t)-{\boldsymbol{\mathfrak{F}}}_k)\right\rangle^\dag
 \label{eq:a-25}\\
 &\equiv\left\langle
 \begin{array}{ccc}
  {\mathbf{F}}_j(t)\otimes e^{i(t'-t)\hat{\mathcal{Q}}\mathfrak{L}}({\mathbf{F}}_k(t)-\boldsymbol{\mathcal{F}}_k)
&& {\mathbf{F}}_j(t)\otimes e^{i(t'-t)\hat{\mathcal{Q}}\mathfrak{L}}({\mathbf{T}}_k(t)-\boldsymbol{\mathcal{T}}_k)\\
  {\mathbf{T}}_j(t)\otimes e^{i(t'-t)\hat{\mathcal{Q}}\mathfrak{L}}({\mathbf{F}}_k(t)-\boldsymbol{\mathcal{F}}_k)
&& {\mathbf{T}}_j(t)\otimes e^{i(t'-t)\hat{\mathcal{Q}}\mathfrak{L}}({\mathbf{T}}_k(t)-\boldsymbol{\mathcal{T}}_k)
  \end{array}
 \right\rangle^\dag\notag.\end{aligned}$$ where $\boldsymbol{\mathcal{F}}\doteq\langle{\mathbf{F}}_k\rangle^\dag$ and $\boldsymbol{\mathcal{T}}\doteq\langle{\mathbf{T}}\rangle^\dag$. The behaviour of Equation (\[eq:a-23\]) as the product $m_b\,{\boldsymbol{\mathfrak{m}}}^{-1}$ vanishes can be obtained by simply letting ${\boldsymbol{\lambda}}\rightarrow \mathbf{0}$. In such a limit, the friction tensor ${\boldsymbol{\gamma}}$ can be approximated by the first term of a multi-power series in ${\boldsymbol{\lambda}}$,[@murphy_brownian_1972] $$\begin{aligned}
 {\boldsymbol{\gamma}}_{jk}({\boldsymbol{\lambda}};t,t')\sim \left\langle
 {\boldsymbol{\mathfrak{f}}}_j(t)\otimes ({\boldsymbol{\mathfrak{f}}}_k(t')-{\boldsymbol{\mathfrak{F}}}_k)
 \right\rangle\label{eq:a-26}\end{aligned}$$ which yields $$\begin{aligned}
 \partial_t{f^{(N)}}(t)+{\boldsymbol{\lambda}}\sum_{j=1}^N&\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\cdot {\widetilde{\boldsymbol{\mathfrak{p}}}}_j+{\boldsymbol{\mathfrak{F}}}_k\cdot{\boldsymbol{\nabla}}_{{\widetilde{\boldsymbol{\mathfrak{p}}}}_j}\right){f^{(N)}}(t) \label{eq:a-27}\\
 =& \sum_{j,k=1}^N {\boldsymbol{\lambda}}^2\,{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}\cdot\int_{-t_I}^{t}dt'\left\langle
 {\boldsymbol{\mathfrak{f}}}_j(t)\otimes ({\boldsymbol{\mathfrak{f}}}_k(t')-{\boldsymbol{\mathfrak{F}}}_k)
 \right\rangle\left({\boldsymbol{\nabla}}_{{\widetilde{\boldsymbol{\mathfrak{p}}}}_k}+\beta\,{\widetilde{\boldsymbol{\mathfrak{p}}}}_k\right){f^{(N)}}(t'). \notag\end{aligned}$$ The last step, the “Markovianization” of this time-evolution equation with memory, is the most controversial one[@snook_langevin_2006; @mazo_theory_1969; @mazur_molecular_1970; @bocquet_microscopic_1997] as it requires the assumption that ${f^{(N)}}(t')$ is very slowly varying compared to the correlation $\langle {\boldsymbol{\mathfrak{f}}}_j(t)\otimes{\boldsymbol{\mathfrak{f}}}_k(t) \rangle^\dag$, so that[@murphy_brownian_1972] $$\int_{-t_I}^{t}dt'\left\langle
 {\boldsymbol{\mathfrak{f}}}_j(t)\otimes ({\boldsymbol{\mathfrak{f}}}_k(t')-{\boldsymbol{\mathfrak{F}}}_k)
 \right\rangle\left({\boldsymbol{\nabla}}_{{\widetilde{\boldsymbol{\mathfrak{p}}}}_k}+\beta\,{\widetilde{\boldsymbol{\mathfrak{p}}}}_k\right){f^{(N)}}(t')\sim \int_{-\infty}^{t}dt'\left\langle
 {\boldsymbol{\mathfrak{f}}}_j(t)\otimes ({\boldsymbol{\mathfrak{f}}}_k(t')-{\boldsymbol{\mathfrak{F}}}_k)
 \right\rangle\left({\boldsymbol{\nabla}}_{{\widetilde{\boldsymbol{\mathfrak{p}}}}_k}+\beta\,{\widetilde{\boldsymbol{\mathfrak{p}}}}_k\right){f^{(N)}}(t). \label{eq:a-28}$$ Such an approximation is widely known to produce an incorrect description of the velocity correlation function if the colloidal particles have a similar density to that of the bath particles[@snook_langevin_2006; @mazo_theory_1969; @bocquet_microscopic_1997; @mazur_molecular_1970]. Equation (\[eq:a-28\]) gives an exponential decay, $\sim e^{-t}$, for the velocity correlation at large times, while with memory, such a decay is algebraic[@snook_langevin_2006], $\sim t^{-3/2}$. Nevertheless, it has also been pointed out that such long-time tails are very small compared to the exponential component predicted under the Markovianized theory[@snook_langevin_2006]. Although these could be reasons to avoid this critical step, the advantages of getting an FPE are significant. In contrast, the approximation is completely valid when both $m_b/m$ and $N_{\text{Kn}}\doteq r_0/R_0$ (the Knudsen number, with $r_0$ being a characteristic length scale for fluid intermolecular interactions and $R_0$ is a characteristic colloidal particle length scale) are considered very small[@peters_fokker-planck_1999]. Nevertheless, it was argued by @bocquet_microscopic_1997 that under such circumstances sedimentation of colloidal particles could occur. This seems to restrict the applicability of the resultant theory to microgravity scenarios. However, it is not clear at all what the actual significance of these algebraic long-time tails is and each case should be judged on its own merits[@mazo_theory_1969]. Thus Equation (\[eq:a-28\]) comprises, on the one hand, an uncontrolled approximation. On the other hand, it has been extensively used in statistical mechanics[@curtiss_kinetic_1957; @wolynes_dynamical_1977; @peters_fokker-planck_1999; @hernandez-contreras_brownian_1996]. For instance, the same hypothesis underlies a recent derivation of a unified DDFT to include inertia and HIs[@goddard_unification_2013], and is also involved in modern theories describing nucleation of colloidal systems and macromolecules[@lutsko_dynamical_2012; @lutsko_classical_2013; @duran-olivencia_mesoscopic_2015; @lutsko_two-parameter_2015]. Moreoever, the results obtained are in perfect agreement with experiments and simulations, corroborating the smallness of the error related to Equation (\[eq:a-28\]) when it comes to describing systems of interacting colloidal particles. More akin to the problem at hand, the hypothesis is tacitly assumed within the seminal work of @dickinson_brownian_1985, where a generalised algorithm to simulate protein diffusional problems is proposed. Therefore, while we cannot really justify such an assumption in a rigorous manner it does, nevertheless, represents the state-of-the-art in modelling colloidal systems. With this in mind, we can finally obtain the FPE related to Equation (\[eq:a-27\]) when Equation (\[eq:a-28\]) is taken into consideration, $$\begin{aligned}
 \partial_t{f^{(N)}}(t)+\sum_{j=1}^N\left({\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\cdot {\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_j+{\boldsymbol{\mathfrak{F}}}_j\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}\right){f^{(N)}}(t) = \sum_{j,k=1}^N \,{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}\cdot{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N)\left({\boldsymbol{\mathfrak{p}}}_k+k_BT{\boldsymbol{\mathfrak{m}}}\,{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_k}\right){f^{(N)}}(t),
 \label{eq:a-29}\end{aligned}$$ with the translational, rotational and coupled translational-rotational friction tensors, $${\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N)
\equiv
 \left(
 \begin{array}{cc}
  {\boldsymbol{\Gamma}}_{jk}^{TT} & {\boldsymbol{\Gamma}}_{jk}^{TR}\\
  {\boldsymbol{\Gamma}}_{jk}^{RT} & {\boldsymbol{\Gamma}}_{jk}^{RR}\\
 \end{array}
 \right)
\doteq \beta {\boldsymbol{\mathfrak{m}}}^{-1}\int_{0}^{\infty}ds\,\left\langle
{\boldsymbol{\mathfrak{f}}}_j(t)\otimes ({\boldsymbol{\mathfrak{f}}}_k(t-s)-{\boldsymbol{\mathfrak{F}}}_k)
 \right\rangle. \label{eq:a-30}$$

Equations of motion
-------------------

The FPE previously derived can be rewritten in the less compact but more explicit form, $$\begin{aligned}
\partial_t {f^{(N)}}(t)+&\sum_{j=1}^N
 {\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_j}\cdot
 \left({\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_j{f^{(N)}}(t)\right)
 -
  {\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j}
 \cdot \left[\left({\boldsymbol{\mathfrak{F}}}_j-\sum_{k=1}^N{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N){\boldsymbol{\mathfrak{p}}}_k
 \right){f^{(N)}}(t)\right]\label{eq:a-31}\\
 &= \sum_{j,k=1}^N
\left[
{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_j} \otimes
{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{p}}}_k}
\right]:\left(k_BT{\boldsymbol{\mathfrak{m}}}\,{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N){f^{(N)}}(t)\right)\notag\end{aligned}$$ which is equivalent to the system of SDEs[@risken_fokker-planck_1996; @kampen_stochastic_2011; @berendsen_simulating_2007] $$\begin{aligned}
 \dot{{\boldsymbol{\mathfrak{r}}}}_j(t)=&\ {\boldsymbol{\mathfrak{m}}}^{-1}{\boldsymbol{\mathfrak{p}}}_j(t)\label{eq:a-32}\\
 \dot{{\boldsymbol{\mathfrak{p}}}}_j(t)=&\ {\boldsymbol{\mathfrak{F}}}_j-\sum_{k=1}^N{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N){\boldsymbol{\mathfrak{p}}}_k + \sum_{k=1}^N\mathbf{A}_{jk}\boldsymbol{\xi}_k(t)\notag\end{aligned}$$ where ${\boldsymbol{\xi}}_j=(\mathbf{f}_j,\mathbf{t}_j)^\top$ is a 6-dimensional Gaussian white noise representing the random forces, $\mathbf{f}_j$, and torques, $\mathbf{t}_j$, acting upon the $j$-th particle, such that $\langle
\xi_j^a(t)\rangle=0$ and $\langle \xi_j^a(t)\xi_k^b(t')\rangle  =
2\delta_{jk}\delta^{ab}\delta(t-t')$, where $\langle.\rangle$ refers to the average over an ensemble of the white-noise realisations. The strength of these random forces and torques is given by the tensor $\mathbf{A}_{jk}$ which obeys the fluctuation-dissipation relation, $$k_BT {\boldsymbol{\mathfrak{m}}}\,{\boldsymbol{\Gamma}}_{jk}({\boldsymbol{\mathfrak{r}}}^N) = \sum_{l=1}^N\mathbf{A}_{jl}({\boldsymbol{\mathfrak{r}}}^N)\mathbf{A}_{kl}({\boldsymbol{\mathfrak{r}}}^N)
\label{eq:a-33},$$ and $$\begin{aligned}
 \mathbf{A}_{jk}{\boldsymbol{\xi}}_k\equiv
 \begin{pmatrix}
\mathbf{A}_{jk}^{TT}&\mathbf{A}_{jk}^{TR}\\
\mathbf{A}_{jk}^{RT}&\mathbf{A}_{jk}^{RR}
 \end{pmatrix}
 \begin{pmatrix}
  \mathbf{f}_{k}\\
  \mathbf{t}_{k}
 \end{pmatrix}.
 \label{eq:a-34}\end{aligned}$$ Coming back to the expanded notation, the system of equations (\[eq:a-32\]) becomes $$\begin{aligned}
 \frac{d{\mathbf{R}}_j}{dt}=&\,\frac{1}{m}{\mathbf{P}}_j,\label{eq:a-35}\\
 \frac{d{\boldsymbol{\Phi}}_j}{dt}=&\,\mathbb{I}^{-1}{\mathbf{L}}_j\label{eq:a-36},\\
 \frac{d{\mathbf{P}}_j}{dt}=&\, -\frac{\partial}{\partial {\mathbf{R}}_j}\widetilde{V}({\mathbf{R}}^N,{\boldsymbol{\alpha}}^N)-\sum_{k=1}^N\left({\boldsymbol{\Gamma}}_{jk}^{TT}{\mathbf{P}}_k+{\boldsymbol{\Gamma}}_{jk}^{TR}{\mathbf{L}}_k\right)+\sum_{k=1}^N\mathbf{A}_{jk}^{TT}\mathbf{f}_k(t)+\mathbf{A}_{jk}^{TR}\mathbf{t}_k(t)
 \label{eq:a-37}\\
 \frac{d{\mathbf{L}}_j}{dt}\equiv\mathbb{I}\frac{d{\boldsymbol{\omega}}_j}{dt}=&
 -\frac{\partial}{\partial {\boldsymbol{\Phi}}_j}\widetilde{V}({\mathbf{R}}^N,{\boldsymbol{\alpha}}^N)-{\boldsymbol{\omega}}_j\times{\mathbf{L}}_j
 \,
 -\sum_{k=1}^N\left({\boldsymbol{\Gamma}}_{jk}^{RT}{\mathbf{P}}_k+{\boldsymbol{\Gamma}}_{jk}^{RR}{\mathbf{L}}_k\right)+\sum_{k=1}^N\mathbf{A}_{jk}^{RT}\mathbf{f}_k(t)+\mathbf{A}_{jk}^{RR}\mathbf{t}_k(t)
 \label{eq:a-38}\end{aligned}$$ with Equation (\[eq:a-36\]) equivalent to the relation ${\boldsymbol{\omega}}_j={\boldsymbol{\Lambda}}_j^\top\dot{{\boldsymbol{\alpha}}_j}$, as can be verified by using Equation (\[eq:a-15\]). Equations (\[eq:a-35\]-\[eq:a-38\]), a much less convenient representation of the rotational-translational Langevin equation (\[eq:a-32\]), have been introduced in the studies by @wolynes_dynamical_1977, @dickinson_brownian_1985[@dickinson_brownian_1985; @dickinson_brownian_1985-1] and @hernandez-contreras_brownian_1996 but these authors started from postulated equations instead of the detailed microscopic derivation from the full system of bath and colloidal particles offered here.

On the rapid relaxation of the fluid: neglecting inertia in the bath {#app:neglecting-bath-inertia}
====================================================================

Here we address the question of when the inertia of the fluid bath can be neglected while having finite viscous forces. With this aim we will make use of some of the results derived by @peters_fokker-planck_1999 in his study on the FPE for coupled rotational and translational motions of structured Brownian particles. The main conclusion, the rapid relaxation of the fluid bath, is ultimately connected with the assumption of negligible inertial effects in the fluid bath whilst considering inertial effects of the colloids.

In a very detailed study, @peters_fokker-planck_1999 applied the multiple time-scale expansion to the derivation of the FPE for arbitrary-shape colloids. With this method @peters_fokker-planck_1999 showed that Equation (\[eq:a-29\]) is the formal time-evolution equation for the distribution function, up to $(m/M)^3$, when both $m/M$ and $N_{\text{Kn}}$ are considered small. It was also argued that the rapid relaxation of the fluid depends upon $N_{\text{Kn}}$ for a system with $m/M$ small. If we now make use of the fact that $m \sim r_0^3\rho_b$ and $M\sim R_0^3\rho_B$, the conditions for the FPE would to be a good description of the colloidal system can be reduced to requiring that $\rho_b/\rho_B$ must be small, which is indeed the condition pointed out by @bocquet_microscopic_1997 and many others[@roux_brownian_1992; @hauge_fluctuating_1973; @hinch_application_1975; @masters_time-scale_1986; @michaels_long-time_1975]. In the following we analyse whether or not this is possible while the ratio between inertia in the fluid bath and viscous forces is small. That is, whether is possible to have a low Reynolds number for the fluid bath along with the condition on the ratio between bath and colloidal densities. In such a case, it would be justified to neglect inertial effects in the bath.

To this end, we consider the case of solid and spherical colloids (far from walls) so that the coupling components of the friction and diffusion tensors vanish, and neglect HIs for the moment. If the radius of the colloidal particles is denoted by $R_0$, the friction tensor takes the simple form[@dickinson_brownian_1985], $$\begin{aligned}
 {\boldsymbol{\Gamma}}_{jk}^{TT} = \gamma_{T}\delta_{jk}\,\mathbf{1},\quad {\boldsymbol{\Gamma}}_{jk}^{RR}=\gamma_{R}\delta_{jk}\,\mathbf{1},\quad {\boldsymbol{\Gamma}}_{jk}^{TR}={\boldsymbol{\Gamma}}_{jk}^{RT}=0,\label{eq:b-1}\end{aligned}$$ with $\gamma_{T}=6\pi\eta\, R_0/M$ and $\gamma_{R}=8\pi\eta\,R_0^3/I=20\pi\eta\,R_0/M$, and $\eta$ the dynamic viscosity, satisfying the Stokes-Einstein formula[@dickinson_brownian_1985]. For the sake of generality, the friction components will be considered equally important, i.e. $\gamma_{T}/\gamma_{R}\sim\mathcal{O}(1)$. For this reason, we can define the following two time scales $$\begin{aligned}
 t_0^T = \left(\frac{MR_0^2\,\gamma_{TT}}{k_BT}\right),\quad t_0^R=\left(\frac{MR_0^2\,\gamma_{RR}}{k_BT}\right)
 \label{eq:b-2},\end{aligned}$$ which are indeed of the same order-of-magnitude. Considering the natural physical scales of the system, the following dimensionless translational and rotational variables (denoted by an asterisk) are pertinent $$\begin{aligned}
\begin{array}{ccccc}
 t = t_0^T\,t*,\ &\ {\mathbf{R}}_j = R_0\,{\mathbf{R}}^*_j,\ &\ {\mathbf{P}}_j =  \frac{k_BT}{\gamma_TR_0}\,{\mathbf{P}}_j^*, \ &\ \mathbf{F}_j = \frac{k_BT}{R_0}\,\mathbf{F}_j^*,\ &\ \mathbf{F}_{j,\text{noise}}= \frac{k_BT}{R_0}\,\mathbf{F}_{j,\text{noise}}^*,\ \\
 t = t_0^R\,t*,\ &\ {\boldsymbol{\omega}}_j = \frac{k_BT}{MR_0^2\gamma_R}{\boldsymbol{\omega}}^*_j,\ &\ {\mathbf{L}}_j =  \frac{k_BT}{\gamma_R}\,{\mathbf{L}}_j^*, \ &\ \mathbf{T}_j = k_BT\,\mathbf{T}_j^*,\ &\ \mathbf{T}_{j,\text{noise}}= k_BT\,\mathbf{T}_{j,\text{noise}}^*,\ \\
\end{array}\label{eq:b-3}\end{aligned}$$ which, when applied in Equations (\[eq:a-35\])-(\[eq:a-38\]), yield the dimensionless equations of motion, $$\begin{aligned}
 \frac{d{\mathbf{R}}_j^*}{dt^*}=&\,{\mathbf{P}}_j^*,\label{eq:b-4}\quad  \frac{d{\boldsymbol{\Phi}}_j^*}{dt^*}=\,{\boldsymbol{\omega}}_j^*\\
 \frac{d{\mathbf{P}}_j^*}{dt^*}=&\,\gamma_T^{*^2}\left(\mathbf{F}_j^*-{\mathbf{P}}_j^*+\mathbf{F}_{j,\text{noise}}^*(t)\right)\notag\\
 \frac{d{\mathbf{L}}_j^*}{dt^*}=&\,\gamma_R^{*^2}\left(\mathbf{T}_j^*-{\mathbf{L}}_j^*+\mathbf{T}_{j,\text{noise}}^*(t)\right)\notag\end{aligned}$$ along with the definitions $$\begin{aligned}
 \gamma_T^{*}=\gamma_T R_0\sqrt{\frac{M}{k_BT}},\quad \gamma_R^*=\gamma_{R} R_0\sqrt{\frac{M}{k_BT}}.
 \label{eq:b-5}\end{aligned}$$ What we wish to test is whether or not it is possible to have finite viscous forces, i.e. $\gamma_T^{*}\sim\mathcal{O}(1)$, such that at the same time inertia forces in the fluid bath are negligible. For this purpose both translational and rotational Reynolds numbers must be small. They can be defined as[@happel_low_1981], $$\begin{aligned}
\text{Re}^T=\frac{U\,R_0\,\rho_b}{\eta},\quad \text{Re}^R=\frac{\Omega\,R_0^2\,\rho_b}{\eta},
\label{eq:b-6}\end{aligned}$$ respectively. The quantities $U$ and $\Omega$ represent typical linear and angular velocity scales, which can be obtained from the nondimensionalised momenta, $$U=\frac{k_BT}{MR_0\gamma_{T}},\quad \Omega=\frac{k_BT}{MR_0^2\gamma_{R}}.\label{eq:b-7}$$ Making use of (\[eq:b-2\]), (\[eq:b-3\]), (\[eq:b-5\]) and (\[eq:b-7\]) into (\[eq:b-6\]), we finally reach $$\text{Re}^T=\frac{9}{2\gamma_{T}^{*^2}}\frac{\rho_b}{\rho_B},\quad \text{Re}^R=\frac{15}{\gamma_{R}^{*^2}}\frac{\rho_b}{\rho_B}\label{eq:b-8}.$$ As already noted, the regime we consider involves, $\gamma_T^{*}\sim\gamma_{R}^*\sim\mathcal{O}(1)$ and $\text{Re}^T\sim\text{Re}^R\ll1$, which is undoubtedly satisfied when $\rho_b/\rho_B\ll 1$. Then, neglecting inertia forces in the fluid bath is consistent with the separation of time scales we already assumed to obtain the FPE. Further details of the physical interpretation and consequences of this limiting condition were given in Section \[app:projection-operator-tech\].

References {#references .unnumbered}
==========

[86]{} \[1\][\#1]{} \[1\][\#1]{} urlstyle \[1\][DOI \#1]{}

Archer, A.J.: Dynamical density functional theory for molecular and colloidal fluids: [A]{} microscopic approach to fluid mechanics. J. Chem. Phys. **130**(1), 014509 (2009)

Archer, A.J., Evans, R.: Dynamical density functional theory and its application to spinodal decomposition. J. Chem. Phys. **121**(9), 4246–4254 (2004)

Barrat, J., Hansen, J.: Basic [Concepts]{} for [Simple]{} and [Complex]{} [Liquids]{}. Cambridge University Press (2003)

Bechtel, D.B., Bulla, L.A.: Electron [Microscope]{} [Study]{} of [Sporulation]{} and [Parasporal]{} [Crystal]{} [Formation]{} in [Bacillus]{} thuringiensis. J. Bacteriol. **127**(3), 1472–1481 (1976)

Beenakker, C.W.J., Saarloos, W.v., Mazur, P.: Many-sphere hydrodynamic interactions. Physica A **127**(3), 451 – 472 (1984)

Berendsen, H.J.C.: Simulating the [Physical]{} [World]{}: [Hierarchical]{} [Modeling]{} from [Quantum]{} [Mechanics]{} to [Fluid]{} [Dynamics]{}. Cambridge University Press (2007)

Bernstein, D.S.: Matrix [Mathematics]{}: [Theory]{}, [Facts]{}, and [Formulas]{} with [Application]{} to [Linear]{} [Systems]{} [Theory]{}. Princeton University Press (2005)

Bocquet, L., Piasecki, J.: Microscopic derivation of non-[Markovian]{} thermalization of a [Brownian]{} particle. J. Stat. Phys. **87**(5-6), 1005–1035 (1997)

Boyd, J.P.: Chebyshev and [Fourier]{} [Spectral]{} [Methods]{}: [Second]{} [Revised]{} [Edition]{}. Courier Corporation (2001)

Brenner, H.: The [Stokes]{} resistance of an arbitrary particle[II]{}. Chem. Eng. Sci. **19**(9), 599 – 629 (1964)

Brown, R.: A brief account of microscopical observations made in the months of [June]{}, [July]{} and [August]{}, 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies. Phil. Mag. **4**, 161–173 (1828)

Cantaert, B., Beniash, E., Meldrum, F.C.: Nanoscale [Confinement]{} [Controls]{} the [Crystallization]{} of [Calcium]{} [Phosphate]{}: [Relevance]{} to [Bone]{} [Formation]{}. Chem.-Eur. J. **19**(44), 14918–14924 (2013)

Chan, G.K.L., Finken, R.: Time-[Dependent]{} [Density]{} [Functional]{} [Theory]{} of [Classical]{} [Fluids]{}. Phys. Rev. Lett. **94**(18), 183001 (2005)

Condiff, D.W., Brenner, H.: Transport [Mechanics]{} in [Systems]{} of [Orientable]{} [Particles]{}. Phys. Fluids **12**(3), 539–551 (1969)

Condiff, D.W., Dahler, J.S.: Brownian [Motion]{} of [Polyatomic]{} [Molecules]{}: [The]{} [Coupling]{} of [Rotational]{} and [Translational]{} [Motions]{}. J. Chem. Phys. **44**(10), 3988–4004 (1966)

Curtiss, C.F., Muckenfuss, C.: Kinetic [Theory]{} of [Nonspherical]{} [Molecules]{}. [II]{}. J. Chem. Phys. **26**(6), 1619–1636 (1957)

Dahler, J.S., Sather, N.F.: Kinetic [Theory]{} of [Loaded]{} [Spheres]{}. [I]{}. J. Chem. Phys. **38**(10), 2363–2382 (1963)

Darve, E., Solomon, J., Kia, A.: Computing generalized [Langevin]{} equations and generalized [Fokker]{}[Planck]{} equations. Proceedings of the National Academy of Sciences **106**(27), 10,884–10,889 (2009)

Deutch, J.M., Oppenheim, I.: Molecular [Theory]{} of [Brownian]{} [Motion]{} for [Several]{} [Particles]{}. J. Chem. Phys. **54**(8), 3547–3555 (1971)

Dickinson, E.: Brownian dynamic with hydrodynamic interactions: the application to protein diffusional problems. Chem. Soc. Rev. **14**(4), 421–455 (1985)

Dickinson, E., Allison, S.A., McCammon, J.A.: Brownian dynamics with rotationtranslation coupling. J. Chem. Soc., Faraday Trans. 2 **81**(4), 591–601 (1985)

Dieterich, W., Frisch, H.L., Majhofer, A.: Nonlinear diffusion and density functional theory. Z. Phys. B Cond. Mat. **78**(2), 317–323 (1990)

Donev, A., Vanden-Eijnden, E.: Dynamic density functional theory with hydrodynamic interactions and fluctuations. J. Chem. Phys. **140**(23), 234115 (2014)

Dunkl, C.F., Xu, Y.: Classical and [Generalized]{} [Classical]{} [Orthogonal]{} [Polynomials]{}. In: Orthogonal [Polynomials]{} of [Several]{} [Variables]{}, Encyclopedia of [Mathematics]{} and its [Applications]{}. Cambridge University Press (2001)

Dur[á]{}n-Olivencia, M.A., Lutsko, J.F.: Mesoscopic nucleation theory for confined systems: [A]{} one-parameter model. Phys. Rev. E **91**(2), 022402 (2015)

Ermak, D.L., McCammon, J.A.: Brownian dynamics with hydrodynamic interactions. J. Chem. Phys. **69**(4), 1352–1360 (1978)

Espa[ñ]{}ol, P., L[ö]{}wen, H.: Derivation of dynamical density functional theory using the projection operator technique. J. Chem. Phys. **131**(24), 244101 (2009)

Evans, G.T.: Cumulant expansion of a [Fokker]{}[Planck]{} equation: [Rotational]{} and translational motion in dense fluids. J. Chem. Phys. **65**(8), 3030–3039 (1976)

Evans, G.T.: Momentum space diffusion equations for chain molecules. J. Chem. Phys. **72**(7), 3849–3858 (1980)

Evans, R.: The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids. Adv. Phys. **28**(2), 143–200 (1979)

Goddard, B.D., Pavliotis, G.A., Kalliadasis, S.: The [Overdamped]{} [Limit]{} of [Dynamic]{} [Density]{} [Functional]{} [Theory]{}: [Rigorous]{} [Results]{}. Multiscale Model. Simul. **10**(2), 633–663 (2012)

Goddard, B.D., Nold, A., Kalliadasis, S.: Multi-species dynamical density functional theory. J. Chem. Phys. **138**(14), 144904 (2013)

Goddard, B.D., Nold, A., Savva, N., Pavliotis, G.A., Kalliadasis, S.: General [Dynamical]{} [Density]{} [Functional]{} [Theory]{} for [Classical]{} [Fluids]{}. Phys. Rev. Lett. **109**(12), 120603 (2012)

Goddard, B.D., Nold, A., Savva, N., Yatsyshin, P., Kalliadasis, S.: Unification of dynamic density functional theory for colloidal fluids to include inertia and hydrodynamic interactions: derivation and numerical experiments. J. Phys.: Condens. Matter **25**(3), 035101 (2013)

Goldstein, H., Poole, C.P., Safko, J.L.: Classical [Mechanics]{}. Addison Wesley (2002)

G[ó]{}mez-Morales, J., Iafisco, M., Delgado-L[ó]{}pez, J.M., Sarda, S., Drouet, C.: Progress on the preparation of nanocrystalline apatites and surface characterization: [Overview]{} of fundamental and applied aspects. Prog. Cryst. Growth Charact. Mater. **59**(1), 1–46 (2013)

Grabert, H., H[ä]{}nggi, P., Talkner, P.: Microdynamics and nonlinear stochastic processes of gross variables. J. Stat. Phys. **22**(5), 537–552 (1980)

Grad, H.: On the kinetic theory of rarefied gases. Commun. Pure Appl. Math. **2**(4), 331–407 (1949)

Gray, C.G., Gubbins, K.E.: Theory of [Molecular]{} [Fluids]{}: [I]{}: [Fundamentals]{}. OUP Oxford (1984)

Happel, J., Brenner, H.: Low [Reynolds]{} number hydrodynamics, *Mechanics of fluids and transport processes*, vol. 1. Springer Netherlands, Dordrecht (1981)

H[ä]{}rtel, A., Blaak, R., L[ö]{}wen, H.: Towing, breathing, splitting, and overtaking in driven colloidal liquid crystals. Phys. Rev. E **81**(5), 051703 (2010)

Hauge, E.H., Martin-L[ö]{}f, A.: Fluctuating hydrodynamics and [Brownian]{} motion. J. Stat. Phys. **7**(3), 259–281 (1973)

Hern[á]{}ndez-Contreras, M., Medina-Noyola, M.: Brownian motion of interacting nonspherical tracer particles: [General]{} theory. Phys. Rev. E **54**(6), 6573–6585 (1996)

Hinch, E.J.: Application of the [Langevin]{} equation to fluid suspensions. J. Fluid Mech. **72**(03), 499–511 (1975)

Hopkins, P., Fortini, A., Archer, A.J., Schmidt, M.: The van [Hove]{} distribution function for [Brownian]{} hard spheres: [Dynamical]{} test particle theory and computer simulations for bulk dynamics. J. Chem. Phys. **133**(22), 224505 (2010)

Jos[é]{}, J.V., Saletan, E.J.: Classical [Dynamics]{}: [A]{} [Contemporary]{} [Approach]{}. Cambridge University Press Textbooks (2013)

Kampen, N.G.V.: Stochastic [Processes]{} in [Physics]{} and [Chemistry]{}. Elsevier (2011)

Kirkwood, J.G.: The [Statistical]{} [Mechanical]{} [Theory]{} of [Transport]{} [Processes]{} [I]{}. [General]{} [Theory]{}. J. Chem. Phys. **14**(3), 180–201 (1946)

Koopmann, R., Cupelli, K., Redecke, L., Nass, K., DePonte, D.P., White, T.A., Stellato, F., Rehders, D., Liang, M., Andreasson, J., Aquila, A., Bajt, S., Barthelmess, M., Barty, A., Bogan, M.J., Bostedt, C., Boutet, S., Bozek, J.D., Caleman, C., Coppola, N., Davidsson, J., Doak, R.B., Ekeberg, T., Epp, S.W., Erk, B., Fleckenstein, H., Foucar, L., Graafsma, H., Gumprecht, L., Hajdu, J., Hampton, C.Y., Hartmann, A., Hartmann, R., Hauser, G., Hirsemann, H., Holl, P., Hunter, M.S., Kassemeyer, S., Kirian, R.A., Lomb, L., Maia, F.R.N.C., Kimmel, N., Martin, A.V., Messerschmidt, M., Reich, C., Rolles, D., Rudek, B., Rudenko, A., Schlichting, I., Schulz, J., Seibert, M.M., Shoeman, R.L., Sierra, R.G., Soltau, H., Stern, S., Str[ü]{}der, L., Timneanu, N., Ullrich, J., Wang, X., Weidenspointner, G., Weierstall, U., Williams, G.J., Wunderer, C.B., Fromme, P., Spence, J.C.H., Stehle, T., Chapman, H.N., Betzel, C., Duszenko, M.: In vivo protein crystallization opens new routes in structural biology. Nat. Methods **9**(3), 259–262 (2012)

Lebowitz, J.L., R[é]{}sibois, P.: Microscopic [Theory]{} of [Brownian]{} [Motion]{} in an [Oscillating]{} [Field]{}; [Connection]{} with [Macroscopic]{} [Theory]{}. Phys. Rev. **139**(4A), A1101–A1111 (1965)

Liboff, R.: Kinetic [Theory]{} - [Classical]{}, [Quantum]{}, and [Relativistic]{} [Descriptions]{}. Graduate [Texts]{} in [Contemporary]{} [Physics]{}, third edn. Springer-Verlag (2003)

Lutsko, J.F.: Recent [Developments]{} in [Classical]{} [Density]{} [Functional]{} [Theory]{}. In: Advances in [Chemical]{} [Physics]{}, pp. 1–92. John Wiley & Sons, Inc. (2010)

Lutsko, J.F.: A dynamical theory of nucleation for colloids and macromolecules. J. Chem. Phys. **136**(3), 034509 (2012)

Lutsko, J.F., Dur[á]{}n-Olivencia, M.A.: Classical nucleation theory from a dynamical approach to nucleation. J. Chem. Phys. **138**(24), 244908 (2013)

Lutsko, J.F., Dur[á]{}n-Olivencia, M.A.: A two-parameter extension of classical nucleation theory. J. Phys.: Condens. Matter **27**(23), 235101 (2015)

Marconi, U.M.B., Tarazona, P.: Dynamic density functional theory of fluids. J. Phys.: Condens. Matter **12**(8A), A413 (2000)

Masters, A.J.: Time-scale separations and the validity of the [Smoluchowski]{}, [Fokker]{}-[Planck]{} and [Langevin]{} equations as applied to concentrated particle suspensions. Mol. Phys. **57**(2), 303–317 (1986)

Mazo, R.M.: On the theory of [Brownian]{} motion. [I]{}. [Interaction]{} between [Brownian]{} particles. J. Stat. Phys. **1**(1), 89–99 (1969)

Mazur, P., Oppenheim, I.: Molecular theory of [Brownian]{} motion. Physica **50**(2), 241–258 (1970)

Michaels, I.A., Oppenheim, I.: Long-time tails and brownian motion. Physica A **81**(2), 221 – 240 (1975)

Miller, W.L., Cacciuto, A.: Hierarchical self-assembly of asymmetric amphiphatic spherical colloidal particles. Phys. Rev. E **80**(2), 021404 (2009)

Murphy, T.J., Aguirre, J.L.: Brownian [Motion]{} of [N]{} [Interacting]{} [Particles]{}. [I]{}. [Extension]{} of the [Einstein]{} [Diffusion]{} [Relation]{} to the [N]{}-[Particle]{} [Case]{}. J. Chem. Phys. **57**(5), 2098–2104 (1972)

Neuhaus, T., H[ä]{}rtel, A., Marechal, M., Schmiedeberg, M., L[ö]{}wen, H.: Density functional theory of heterogeneous crystallization. Eur. Phys. J. Spec. Top. **223**(3), 373–387 (2014)

Nold, A., Sibley, D.N., Goddard, B.D., Kalliadasis, S.: Fluid structure in the immediate vicinity of an equilibrium three-phase contact line and assessment of disjoining pressure models using density functional theory. Phys. Fluids **26**(7), 072001 (2014)

Nold, A., Sibley, D.N., Goddard, B.D., Kalliadasis, S.: Nanoscale [Fluid]{} [Structure]{} of [Liquid]{}-solid-vapour [Contact]{} [Lines]{} for a [Wide]{} [Range]{} of [Contact]{} [Angles]{}. Math. Model. Nat. Phenom. **10**(4), 111–125 (2015)

Peters, M.H.: Fokker-[Planck]{} equation and the grand molecular friction tensor for coupled rotational and translational motions of structured [Brownian]{} particles near structured surfaces. J. Chem. Phys. **110**(1), 528–538 (1999)

Peters, M.H.: The [Smoluchowski]{} diffusion equation for structured macromolecules near structured surfaces. J. Chem. Phys. **112**(12), 5488–5498 (2000)

Pottier, N.: Nonequilibrium [Statistical]{} [Physics]{}: [Linear]{} [Irreversible]{} [Processes]{}. Oxford University Press (2014)

Rex, M., Wensink, H.H., L[ö]{}wen, H.: Dynamical density functional theory for anisotropic colloidal particles. Phys. Rev. E **76**(2), 021403 (2007)

Risken, H.: The [Fokker]{}-[Planck]{} [Equation]{}: [Methods]{} of [Solutions]{} and [Applications]{}. 2nd ed. edn. Springer (1996)

Roux, J.N.: Brownian particles at different times scales: a new derivation of the [Smoluchowski]{} equation. Physica A **188**(4), 526–552 (1992)

Schilling, T., Frenkel, D.: Self-poisoning of crystal nuclei in hard-rod liquids. J. Phys.: Condens. Matter **16**(19), S2029 (2004)

Snook, I.: The [Langevin]{} and [Generalised]{} [Langevin]{} [Approach]{} to the [Dynamics]{} of [Atomic]{}, [Polymeric]{} and [Colloidal]{} [Systems]{}. Elsevier (2006)

van Teeffelen, S., Likos, C.N., L[ö]{}wen, H.: Colloidal [Crystal]{} [Growth]{} at [Externally]{} [Imposed]{} [Nucleation]{} [Clusters]{}. Phys. Rev. Lett. **100**(10), 108302 (2008)

Wilemski, G.: On the derivation of [Smoluchowski]{} equations with corrections in the classical theory of [Brownian]{} motion. J. Stat. Phys. **14**(2), 153–169 (1976)

Wittkowski, R., L[ö]{}wen, H.: Dynamical density functional theory for colloidal particles with arbitrary shape. Mol. Phys. **109**(23-24), 2935–2943 (2011)

Wolde, P.R.t., Frenkel, D.: Enhancement of [Protein]{} [Crystal]{} [Nucleation]{} by [Critical]{} [Density]{} [Fluctuations]{}. Science **277**(5334), 1975–1978 (1997)

Wolynes, P.G., Deutch, J.M.: Dynamical orientation correlations in solution. J. Chem. Phys. **67**(2), 733–741 (1977)

Wu, J., Li, Z.: Density-[Functional]{} [Theory]{} for [Complex]{} [Fluids]{}. Annu. Rev. Phys. Chem. **58**(1), 85–112 (2007)

Yatsyshin, P., Savva, N., Kalliadasis, S.: Spectral methods for the equations of classical density-functional theory: [Relaxation]{} dynamics of microscopic films. J. Chem. Phys. **136**(12), 124113 (2012)

Yatsyshin, P., Savva, N., Kalliadasis, S.: Geometry-induced phase transition in fluids: [Capillary]{} prewetting. Phys. Rev. E **87**(2), 020402(R) (2013)

Yatsyshin, P., Savva, N., Kalliadasis, S.: Density functional study of condensation in capped capillaries. J. Phys.: Condens. Matter **27**(27), 275104 (2015)

Yatsyshin, P., Savva, N., Kalliadasis, S.: Wetting of prototypical one- and two-dimensional systems: [Thermodynamics]{} and density functional theory. J. Chem. Phys. **142**(3), 034708 (2015)

Zhang, Z.X., Duijneveldt, J.S.v.: Isotropic-nematic phase transition of nonaqueous suspensions of natural clay rods. J. Chem. Phys. **124**(15), 154910 (2006)

[^1]: With $b$ and $B$ subscripts referring to the *bath* and *Brownian* (colloidal) particles respectively.
