---
abstract: 'We establish the global well-posedness of overdamped dynamical density functional theory (DDFT): a nonlinear, nonlocal integro-partial differential equation used in statistical mechanical models of colloidal flow and other applications including nonlinear reaction-diffusion systems and opinion dynamics. With no-flux boundary conditions, we determine the well-posedness of the full nonlocal equations including two-body hydrodynamic interactions (HI) through the theory of Fredolm operators. Principally, this is done by rewriting the dynamics for the density $\varrho$ as a nonlocal Smoluchowski equation with a non-constant diffusion tensor $\bm{D}$ dependent on the diagonal part ($\bm{Z}_1$) of the HI tensor, and an effective drift $\bm{A}[\vec{a}]$ dependent on the off-diagonal part ($\bm{Z}_2$). We derive a scheme to uniquely construct the mean colloid flux $\vec{a}(\vec{r},t)$ in terms of eigenvectors of $\bm{D}$, show that the stationary density $\varrho(\vec{r})$ is independent of the HI tensors, as well as proving exponentially fast convergence to equilibrium. The stability of the equilibria $\varrho(\vec{r})$ is studied by considering the bounded (nonlocal) perturbation of the differential (local) part of the linearised operator. We show that the spectral properties of the full nonlocal operator with no-flux boundary conditions can differ considerably from those with periodic boundary conditions. We showcase our results by using the numerical methods available in the pseudo-spectral collocation scheme 2DChebClass.'
author:
- '$^\dagger$,'
- '[^1] ,'
- '[^2]'
bibliography:
- 'bibYear3.bib'
title: '****'
---

dynamic density functional theory (DDFT), colloids, overdamped limit, hydrodynamic interactions, nonlocal-differential PDEs, interacting particle systems, McKean-Vlasov equation, phase transitions, bifurcation theory.

60F10; 60J75; 62P10; 92C37

[^3]

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

For suspended particles in a viscous fluid, the Navier-Stokes equations are not sufficient to model flows on a spatial scale comparable with the size of the individual particles. Instead, one requires a computationally tractable model that captures meso/macro-scale dynamics whilst also including physical effects driven by particle-level interactions. Dynamic density functional theories (DDFTs) are excellent candidates for modelling such systems [@marconi1999dynamic; @ArcherEvans04]. They are typically applied in condensed matter physics in the colloidal particle regime with particles of typical diameters 1nm$-$1$\mu$m. Recent advances have allowed the inclusion of inertia [@MarconiTarazonaCecconiMelchionna07; @archer2009dynamical], multiple species [@Archer05; @RothRauscherArcher09; @GNK13; @LichtnerArcherKlapp12], hydrodynamic interactions (HI) [@rex2009dynamical; @Rauscher10; @goddard2012general; @goddard2012unification], background flows [@RauscherDominguezKrugerPenna07], temperature gradients [@WittkowskiLowenBrand12; @anero2013functional], hard spheres [@Rosenfeld89; @RothEvansLangKahl02; @roth2010fundamental; @StopperMaroltRothHansen-Goos15], confined geometries [@goddard2016dynamical; @zimmermann2016flow], arbitrary shaped particles [@WittkowskiLowen11], and active microswimmers [@menzel2016dynamical; @hoell2017dynamical].

For equilibrium fluids, there is a rigorous mathematical framework proving the existence of nontrivial fluid densities, different from those found by classical fluid dynamical formalisms, by taking into account both many body effects and external force fields. This is commonly known as (classical) density functional theory (DFT) [@Mermin:1965lo]. It is able to predict effects driven by the microscale, e.g., the non-smooth droplet profiles which are formed at the gas-liquid-solid trijunction in contact line problems [@berim2009simple] and the coexistence of multiple fluid films at critical values of the chemical potential energy in droplet spreading [@pereira2012equilibrium]. It has been used to resolve the paradox of stress and pressure singularities normally found in classical moving contact line problems [@sibley2013contact]. What is more, DFT agrees well with molecular dynamics simulations; see, e.g., [@Lutsko10] and references therein. These advancements motivate more mathematical analysis, in particular, on the well-posedness of the underlying equations being used and on the number and structure of equilibrium states.

As a non-equilibrium extension to DFT for classical fluids, dynamic DFT (DDFT) has been applied to a wide range of problems: polymeric solutions [@PennaDzubiellaTarazona03], spinodal decomposition [@ArcherEvans04], phase separation [@Archer05], granular dynamics [@MarconiMelchionna07; @MarconiTarazonaCecconi07], nucleation [@vanTeeffelenLikosLowen08], liquid crystals [@WittkowskiLowenBrand10], and evaporating films [@ArcherRobbinsThiele10]. Recently, a stochastic version of DDFT has been derived [@Lutsko12], which allows the study of energy barrier crossings, such as in nucleation.

A crucial point is that the computational complexity of DDFT is (essentially) constant in the number of particles, which allows the treatment of macroscopically large systems, whilst retaining microscopic information. Furthermore, due to the universality of the underlying nonlinear, nonlocal partial differential equations, DDFT may be considered as a generalisation of a wider class of such models used in the continuum modelling of many natural phenomena consisting of complex, many body, multi-agent interparticle effects including: pattern formation [@camazine2003self], the flocking of birds, cell proliferation, the self organising of morphogenetic and bacterial species [@canizo2010collective; @carrillo2009double], nonlocal reaction-diffusion equations [@al2018dynamical] and even consensus modelling in opinion dynamics[@chazelle2017well]. Many of these applications are often described as systems of interacting (Brownian) particles and, in the case of hard particle viscous suspensions, bath-mediated HI effects may be included.

The HI are forces on the colloids mediated by the bath flow, generated by the motion of the colloidal particles. This in turn produces a nontrivial particle–fluid–particle hydrodynamic phenomenon, the inclusion of which has been shown to have substantial effects on the physics of many systems; for example, they have been found to be the underlying mechanism for the increased viscosity of suspensions compared to a pure bath [@Einstein06], the blurring of laning that arises in driven flow [@WysockiLowen11], the migration of molecules away from a wall [@HodaKumar07], and are particularly complex in confined systems [@happel2012low; @LaugaSquires05], and for active particles and microswimmers, which result in additional HI [@HuberKoehlerYang11].

Mathematically, the inter-particle forces and HI can be described through the hydrodynamic fields $\varrho$ and $\vec{v}$, the one-body density and one-body velocity fields, respectively. These fields, inherent to a continuum description of a collection of particles, are derived by considering successive moments (density, velocity, heat flux, …) of the underlying kinetic system [@gorban2014hilbert]. In particular, for systems of interacting Newtonian particles, when the momenta are non-negligible, the evolution of the phase space density $f(\vec{r}^N,\vec{p}^N, t)$ for a system of $N$ colloids determining the probability of finding the system in the state $(\vec{r}^N,\vec{p}^N)$ at time $t$ is described by the $N$-body Fokker-Planck equation and the dynamics of the hydrodynamic fields are defined by obtaining closed equations for $\{\varrho, \varrho\times \vec{v}\} := \int \mathrm{d}\vec{r}^{N-1}\,\mathrm{d}\vec{p}^N\, \{1, \vec{p}/m\}f(\vec{r}^N,\vec{p}^N, t)$, where $m$ is the particle mass. Here, $\vec{r}^N$ and $\vec{p}^N$ denote the $3N$-dimensional position and momentum vectors of all $N$ particles.

The inclusion of HI leads to a much richer hierachy of fluid equations compared to systems without HI; compare e.g. [@goddard2012unification] and [@archer2009dynamical]. In particular, see e.g. [@goddard2012unification], by integration over all but one particle position, the one-body Fokker-Planck equation may be obtained. If, in addition, two-body HI and interparticle interactions are assumed and the inertia of the colloids is considered small, a high friction limit $\gamma\to \infty$ may be taken [@goddard2012overdamped]. The result is that the velocity distribution converges to a Maxwellian, and one can eliminate the momentum variable through an adiabatic elimination process that is based on multiscale analysis [@pavliotis2008multiscale]. The final one-body Smoluchowski equation for $\varrho$ is a novel, nonlinear, nonlocal PDE shown to be independent of the unknown kinetic pressure term $\int \mathrm{d}\vec{r}\,\mathrm{d}\vec{p}\, m^{-2} \vec{p}\otimes\vec{p} f(\vec{r},\vec{p}, t)$, which normally persists at $\gamma = O(1)$ (see[@goddard2012overdamped], Theorem 4.1).

Existence, uniqueness and global asymptotic stability of the novel Smoluchowski equation in this overdamped limit has, until this work, remained unproven. It is the inclusion of HI that provides richness through additional nonlinearities in both the dissipation and convection terms. The inclusion of HI is interesting from both physical and mathematical standpoints. Physically, as above, the HI give rise to a much more complex evolution in the density. Mathematically, the convergence to equilibrium will depend inherently on the spectral properties of the effective diffusion tensor and effective drift vector arising from the HI. What is more, since the full $N$-body Fokker-Planck equation is a PDE in a very high dimensional phase space, well-posed nonlinear, nonlocal PDEs governing the evolution of the one-particle distribution function, valid in the mean field limit, describing the flow of nonhomogeneous fluids are desirable for computational reasons.

The equations studied in this paper are related to the McKean-Vlasov equation [@chayes2010mckean], a nonlinear nonlocal PDE of Fokker-Planck type that arises in the meanfield limit of weakly interacting diffusions. The novelty of the present problem lies in the space dependent diffusion tensor and nonlinear, nonlocal boundary conditions. Additionally, the problem that we study in this paper may in general not be written as a gradient flow, with the exception of the modelling assumption that the off-diagonal elements of the friction tensor $\bm{\Gamma}$ are zero. This choice is equivalent to setting $\bm{Z}_2$ to zero, and would be physically relevant for a diffuse system of particles with a strong hydrodynamic interaction with a wall but weak inter-particle hydrodynamic interactions [@goddard2016dynamical].

Description of the Model.
-------------------------

In this work we analyse the overdamped partial differential equation (PDE) associated to a system of interacting stochastic differential equations (SDEs) on $U$ an open, bounded subset of $\mathbb{R}^d$ of the following form, governing the positions $\vec{r}_i$ and momenta $\vec{p}_i$ of $i = 1,\dots, N$ colloidal particles immersed in a bath of many more, much smaller and much lighter particles:

$$\begin{aligned}
\frac{\mathrm{d}\vec{r}_i}{\mathrm{d}t} &= \frac{1}{m}\vec{p}_i,\label{eq:SDE_overdamped_pos}\\
\frac{\mathrm{d}\vec{p}_i}{\mathrm{d}t} &= -\nabla_{\vec{r}_i}  V(\vec{r}^N,t)-\sum_{j=1}^{N}\boldsymbol{\Gamma}_{ij}(\vec{r}^N)\vec{p}_j + \sum_{j=1}^N\boldsymbol{B}_{ij}(\vec{r}^N)\vec{f}_{j}(t) \label{eq:SDE_overdamped_mom}\end{aligned}$$

where $\vec{r}^N = (\vec{r}_1,\cdots,\vec{r}_N)$, $\boldsymbol{B} = \left( mk_BT\boldsymbol{\Gamma}\right)^{1/2}$, $\boldsymbol{\Gamma} = \gamma (\bm{1} + \tilde{\boldsymbol{\Gamma}})$ (where the tilde denotes the nondimensional tensor and $\bm{1}$ is the $3N\times 3N$ identity matrix), $V$ is a potential, $k_B$, $\mathrm{T}$, $\gamma$ are Boltzmann’s constant, temperature and friction, respectively, and $\vec{f}_i(t) = (\zeta^x_i(t),\zeta^y_i(t),\zeta^z_i(t))^\top$ is a Gaussian white noise term with mean and correlation given by $\langle \zeta_i^a(t)\rangle = 0$ and $\langle\zeta_i^a(t),\zeta_j^b(t) \rangle = 2\delta_{ij}\delta^{ab}\delta(t-t')$.

In $d = 3$ dimensions, the friction tensor $\boldsymbol{\Gamma}$ comprises $N^2$ positive definite $3\times 3$ mobility matrices $\boldsymbol{\Gamma}_{ij}$ for the colloidal particles. These couple the momenta of the colloidal particles to HI forces on the same particles, mediated by fluid flows in the bath. Typically, in the underdamped limit with dense suspensions, the HI may be short range lubrication forces, whereas in disperse systems in the overdamped limit, the HI are taken to be the long range forces given by the Rotne-Prager-Yamakawa tensor [@rotne1969variational]. However, we do not make any such assumptions on the form of the tensors here.

We have described a general set of coupled Langevin equations with spatially-dependent friction tensor $\boldsymbol{\Gamma}(\vec{r}^N)$. As we will see, the dynamics – tend towards an equilibrium given by the Gibbs probability measure, which we will show to be independent of the friction tensor. Instead of computing the trajectories of individual particles we consider the evolution of the density of particles $\varrho(\vec{r},t)$ given by the Smoluchowski equation in the high friction limit $\gamma\to \infty$, $$\begin{aligned}
\label{eq:mk-eq}
\qquad \partial_{t}\varrho(\vec{r},t) =-\tfrac{k_B\mathrm{T}}{m\gamma} \nabla_{\vec{r}}\cdot\vec{a}(\vec{r},[\varrho], t)  \qquad \text{ for } \vec{r} \in  U,\,t\in [0,T]\end{aligned}$$ where $\vec{a}(\vec{r},[\varrho], t)$ is the flux, $[\varrho]$ denotes functional dependence, $U\subseteq\mathbb{R}^d$ and $T<\infty$. Equation was derived rigorously as a solvability condition of the corresponding Vlasov-Fokker-Planck equation for the one-body density in position and momentum space $f(\vec{r},\vec{p},t)$ by writing $f$ as a Hilbert expansion in a small nondimensional parameter $\epsilon\propto\gamma^{-1}$ [@goddard2012overdamped]. Therein, $\epsilon$ has units length, and therefore a problem specific length scale must be introduced to make it truly nondimensional.

We are interested in global existence, uniqueness, positivity and regularity of the weak solution to when $\vec{a}(\vec{r},t)$ is given by the integral equation

$$\begin{aligned}
\vec{a}(\vec{r},t) + \boldsymbol{H}[\vec{a},\varrho](\vec{r},t)+\frac{\varrho(\vec{r},t)}{k_BT}\bm{D}(\vec{r},[\varrho],t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho](\vec{r},t)
=0,\label{eq:eqn_for_a}\end{aligned}$$

$$\begin{aligned}
\boldsymbol{H}[\vec{a},\varrho](\vec{r},t):=
\varrho(\vec{r},t)\bm{D}(\vec{r},[\varrho],t)\int_U\mathrm{d}\vec{r}'\, g(\vec{r},\vec{r}')\boldsymbol{Z}_2(\vec{r},\vec{r}')\vec{a}(\vec{r}',t),\label{eq:eqn_for_H}\end{aligned}$$

$$\begin{aligned}
&\frac{\varrho(\vec{r},t)}{k_BT}\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho](\vec{r},t) := [\nabla_{\vec{r}}+\tfrac{1}{k_B\mathrm{T}}\Big(\nabla_{\vec{r}}V_1(\vec{r},t)\nonumber\\
&\qquad\qquad\qquad\qquad\qquad\qquad+\int_U\mathrm{d}\vec{r}'\varrho(\vec{r}',t)g(\vec{r},\vec{r}')\nabla_{\vec{r}}V_{2}(\vec{r},\vec{r}')\Big)  ]\varrho(\vec{r},t),\label{eq:eqn_for_J}\end{aligned}$$

where to ease notation we have suppressed $[\varrho]$ in the argument of $\vec{a}$ and $\mathcal{F}$ is the free energy functional which will be defined in Section \[subsec:free\_energy\_framework\]. The functions $V_1$ and $V_2$ are the external and (two body) interparticle potentials respectively. Additionally, the non-constant diffusion tensor $$\begin{aligned}
\label{eq:def_diffusion_tensor}
\boldsymbol{D}(\vec{r},[\varrho],t):=\frac{k_{\text{B} }\mathrm{T}}{m \gamma}\Big[\boldsymbol{1}+\int\mathrm{d}\vec{r}'g(\vec{r},\vec{r}')\varrho(\vec{r}',t)\boldsymbol{Z}_1(\vec{r},\vec{r}')\Big]^{-1}\end{aligned}$$ will be considered; this is interesting from a physical point of view. It has been previously shown (see [@goddard2012overdamped]) that for $\boldsymbol{Z}_1$ being positive definite, $\boldsymbol{D}$ is also positive definite and therefore has positive, finite eigenvalues. The term $g(\vec{r},\vec{r}')$ (regarded as known) is the correlation function defined by the two-body density $\varrho^{(2)}(\vec{r},\vec{r}',t) =g(\vec{r},\vec{r}')\varrho(\vec{r},t)\varrho(\vec{r}',t)$ and the operator $\boldsymbol{H}[\cdot]$ describes terms corresponding to HI.

We note that if $\boldsymbol{D}$ were positive semidefinite, a zero eigenvalue of $\boldsymbol{D}$ is permitted, which physically-speaking would amount to the colloidal system possessing a zero diffusion rate in some subset of $U$ with nonzero measure. Such systems are interesting (for example, in many biological systems the physical domain $U$ could be a substrate including cuts, voids or interior walls) but are not considered in this paper. Throughout this work the largest and smallest eigenvalues of $\boldsymbol{D}$ will be denoted $\mu_{\max}$ and $\mu_{\min}$, respectively.

Furthermore, for two-body HI, $\boldsymbol{Z}_1$, $\boldsymbol{Z}_2$ are the diagonal and off-diagonal blocks respectively of the translational component of the grand resistance matrix originating in the classical theory of low Reynolds number hydrodynamics between suspended particles [@happel2012low], [@jeffrey1984calculation], related to the friction tensor by $$\begin{aligned}
\tilde{\boldsymbol{\Gamma}}_{ij}(\vec{r}^N) = \delta_{ij}\sum_{l\neq i}\boldsymbol{Z}_{1}(\vec{r}_i,\vec{r}_l)+(1-\delta_{ij})\boldsymbol{Z}_{2}(\vec{r}_i,\vec{r}_j).\end{aligned}$$

In $d = 3$ dimensions, and for the particular case $N=2$ (where $N$ is the number of particles in the system), $\bm{\Gamma}\in \mathbb{R}^{6\times 6}$ and $\bm{\Gamma}_{ij}$ may be seen as equivalent to the second-rank tensor of the translational part of the resistance matrix as found in [@jeffrey1984calculation] used to model lubrication forces. It should be noted however that the definition of those resistance matrices are formalism dependent, that is, the individual entries are scalar functions arising from the solution of Stokes equations for two-body lubrication interactions using multipole methods. Conversely, $\bm{\Gamma}_{ij}$ are general tensors, independent of the type of HI under consideration, and are therefore a more general representation of hydrodynamic phenomena of colloidal suspensions. Additionally, $\bm{\Gamma}_{ij}$ may be used to model not just lubrication forces between particles but also long range forces, wall effects and more. In the case of inter-particle HI, the diagonal blocks $\bm{\Gamma}_{ii}$ each represent the force exerted on the fluid due to the motion of particle $i$, which is simply the sum of all the pairwise HI from the perspective of particle $i$. The off-diagonal blocks $\bm{\Gamma}_{ij}$ represent the force on particle $i$ due to the motion of particle $j$.

The stationary equations for the equilibrium density $\varrho(\vec{r})$ and equilibrium flux $\vec{a}(\vec{r})$ are given by

$$\begin{aligned}
&\qquad\qquad\qquad\qquad\nabla_{\vec{r}}\cdot \vec{a}(\vec{r}) = 0,\label{eq:div_a_0}\\
&\vec{a}(\vec{r}) + \boldsymbol{H}[\vec{a},\varrho](\vec{r})+\frac{\varrho(\vec{r},t)}{k_BT}\bm{D}(\vec{r},[\varrho],t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho](\vec{r})
=0.\label{eq:eqn_for_a_in_equilibrium}
 \end{aligned}$$

Note that given a finite flux vector $\vec{a}$ solving -, it is not obvious that $\varrho$ is necessarily a minimiser of the free energy $\mathcal{F}-\int_U\mathrm{d}\vec{r}\,\mu_{c}\varrho$ (where $\mu_{c}$ is the chemical potential of the species). However, for the particular choice $\vec{a}\equiv \vec{0}$ (which is a natural and physically realistic solution), $\varrho$ is necessarily a minimiser of $\mathcal{F}-\int_U\mathrm{d}\vec{r}\,\mu_{c}\varrho$, and we will show that under reasonable assumptions these are indeed the only fixed points of the system.

Previous well-posedness studies of similar nonlinear, nonlocal PDEs focused on periodic boundary conditions; see, e.g., [@chazelle2017well; @greg_mckean_vlasov_torus]. In contrast, we are interested in the well-posedness of , - subject to no-flux boundary conditions. This choice admits the nontrivial effect of the two body forces generated by the potential $V_2$ interacting with density on the boundary of the physical domain. We also seek to understand the asymptotic stability of stationary states. The motivation for this choice of boundary condition is physical; it corresponds to a closed system of particles in which the particle number is conserved over time. It is clear that most applications of such equations will be in confined systems, rather than a periodic domain and, as such, no-flux boundary conditions are natural. We note that the choice of boundary condition is expected to have significant effects on the dynamics, including the form of the bifurcation diagram.

Free Energy Framework. {#subsec:free_energy_framework}
----------------------

Related to the system -, we define the free energy functional $\mathcal{F}:P^+_{\text{ac}}(U)\to \mathbb{R}$ where $P^+_{\text{ac}}$ is the set of strictly positive definite absolutely continuous probability measures on $U$. We define $$\begin{aligned}
\label{eq:def_of_F}
\mathcal{F}[\varrho]&:=\int_U\mathrm{d}\vec{r}\,\varrho(\vec{r},t)\log\varrho(\vec{r},t)+\int_U\mathrm{d}\vec{r}\,\varrho(\vec{r},t)\,\Big[V_1(\vec{r},t)+\tfrac{1}{2}(gV_2)\star\varrho \Big],\end{aligned}$$ where $\star$ denotes convolution in space. Here we assume the probability measure $\varrho$ has density with respect to the Lebesgue measure. Additionally we define the probability measure on $U$ $$\begin{aligned}
\label{eq:def_of_measure}
\mu(\mathrm{d}\vec{r}) = \mathrm{d}\vec{r}\,Z^{-1}e^{-\tfrac{(V_1+(gV_2)\star\varrho)}{k_BT}}\end{aligned}$$ where $Z = \int_U \mathrm{d}\vec{r}\,e^{-\tfrac{(V_1+(gV_2)\star\varrho)}{k_BT}}$ and $\varrho$ (when it exists) satisfies the nonlinear equation $$\varrho = Z^{-1}e^{-\tfrac{(V_1+(gV_2)\star\varrho)}{k_BT}}.$$ The existence of a probability density $\varrho$, and therefore a probability measure $\mu$ in , is obtained by Lemma \[thm:exis\_fix\_point\]. The functional $\mathcal{F}$ gives rise to the density minimising the free energy associated to the system - as $\gamma \to \infty$, which will be shown in Theorem \[thm:association \_of\_free\_energy\].

To make the connection between the free energy functional $\mathcal{F}$ in and the theory of non-uniform classical fluids, one may consider the Helmholtz free energy functional, which is the central energy functional of DFT [@evans1979nature] $$\begin{aligned}
\label{eq:ddft-helmholtz_func}
\mathcal{F}_{H}[\varrho] = \int_U\mathrm{d}\vec{r}\,\varrho(\vec{r},t)V_1(\vec{r},t)+ k_BT\int_U\mathrm{d}\vec{r}\,\varrho(\vec{r},t)[\log(\Lambda^3\varrho(\vec{r},t))-1]+\mathcal{F}_\text{ex}[\varrho]\end{aligned}$$ where $\mathcal{F}_{\text{ex}}$ is the excess over ideal gas term and $\Lambda$ the de Broglie wavelength, which turns out to be superfluous. The term $\mathcal{F}_{\text{ex}}$ is not in general known, the exception being for one dimensional hard rods [@percus1976equilibrium]. Using the free energy functional $\mathcal{F}_H$, the corresponding Euler-Lagrange equation is $$\begin{aligned}
\label{eq:chemical_potential_euler_lagrange_eqn}
\mu_{c}=V_1(\vec{r}) + k_BT [\log (\Lambda^3\varrho(\vec{r}))-1]+\tfrac{\delta\mathcal{F}_{\text{ex}}}{\delta\rho}[\varrho]\end{aligned}$$ where $\mu_{c}$ is the chemical potential which is constant at equilibrium. Note that $\mu_c$ should not be confused with the measure $\mu$ defined in . After taking the gradient of and multiplying by $\varrho$ we obtain $$\begin{aligned}
0=\varrho(\vec{r})\nabla_{\vec{r}}\tfrac{\delta\mathcal{F}}{\delta\rho}[\varrho]=k_BT\nabla_{\vec{r}}\varrho+\varrho(\vec{r})\nabla_{\vec{r}}\Big(V_1(\vec{r})+\tfrac{\delta\mathcal{F}_{\text{ex}}}{\delta\rho}[\varrho]\Big).\end{aligned}$$ At equilibrium, the sum rule holds (see, e.g. [@goddard2012unification]) $$\begin{aligned}
\label{eq:excess_free_energy_equilibrium_sum}
\varrho(\vec{r})\nabla_{\vec{r}}\tfrac{\delta\mathcal{F}_{\text{ex}}}{\delta\varrho}[\varrho]=\sum_{n=2}^N\int\mathrm{d}\vec{r}^{n}\nabla_{\vec{r}}V_n(\vec{r}^n)\varrho_n(\vec{r}^n).\end{aligned}$$ where $\varrho_n(\vec{r}^n)$ is the standard $n-$particle configuration distribution function in equilibrium. Limiting the particle interactions to two-body, for example with the approximation $\varrho_2(\vec{r},\vec{r}') = \varrho(\vec{r})\varrho(\vec{r}')g(\vec{r},\vec{r}',[\varrho])$, we take the first term in the above series to obtain the equality $\nabla_{\vec{r}}\mathcal{F}_H[\varrho]= \nabla_{\vec{r}}\mathcal{F}[\varrho]$. In this way wee see that the density minimising $\mathcal{F}_H$ will minimise $\mathcal{F}$.

When $\boldsymbol{Z}_2\equiv 0 $, and by using the adiabatic approximation that holds out of equilibrium, we note that PDE simplifies to (cf. [@rex2009dynamical]) $$\begin{aligned}
\label{eq:ddft-eq}
\partial_{t}\varrho = \nabla_{\vec{r}}\cdot \left[\boldsymbol{D}(\vec{r},t)\varrho(\vec{r},t)\,\nabla_{\vec{r}}\tfrac{\delta \mathcal{F}}{\delta \varrho}[\varrho]\right].\end{aligned}$$ From we conclude that the dynamics under the choice $\bm{Z}_2 \equiv 0$ has a gradient flow structure. When $\bm{Z}_2$ is not necessarily zero, one cannot in general write the full dynamics as a gradient flow and, hence, the inclusion of HI introduces a novel perturbation away from the classical theory of gradient flow structure. Additionally, one sees how the free energy functional gives rise to the concept of a local pressure variation by the term inside the divergence of . In particular, the term $\tfrac{k_\text{B}\mathrm{T}}{m}\varrho(\vec{r},t)\,\nabla_{\vec{r}}\tfrac{\delta \mathcal{F}}{\delta \varrho}[\varrho]$ represents the spatial variation of the energy available to change particle configurations per unit volume at fixed particle number, in other words, it is an analogue of a local pressure gradient for the particle density. We will show that $\mathcal{F}[\varrho]$ is associated to the PDE even when $\bm{Z}_2 \neq 0$, that is $\partial_t \varrho = 0$ implies $\varrho$ is a critical point of $\mathcal{F}$.

[0.48]{} ![(a). The bifurcation diagram for (a). $V_2(x,y) = x\cdot y$ and (b). $V_2(x,y) = -\cos\left(\frac{2\pi(x-y)}{\mathrm{L}}\right)$ in Section \[subsec:numericalexperiments\]: the solid [blue]{} line denotes the stable branch of solutions while the dotted [red]{} line denotes the unstable branch of solutions. In (a) the stationary density $e^{-x^2}/Z$ changes stability at the critical interaction energy $\kappa_2 = \kappa_{2\sharp}= -2.4$ and the new stable density is asymmetric adhering to one wall (Figure \[fig:bif\_fig\_right\]). In (b), in the absence of a confining potential, the uniform density becomes unstable at the critical interaction energy $\kappa_2 = \kappa_{2\sharp} = 0.4$ and the density may become multi-modal (Figure \[fig:bif\_stable\_equilibria\_V2\_cosine\]).[]{data-label="fig:bifurcation_diagram"}](bifurcation_diagram_for_Phi2.pdf "fig:"){width="\textwidth"}

 

[0.48]{} ![(a). The bifurcation diagram for (a). $V_2(x,y) = x\cdot y$ and (b). $V_2(x,y) = -\cos\left(\frac{2\pi(x-y)}{\mathrm{L}}\right)$ in Section \[subsec:numericalexperiments\]: the solid [blue]{} line denotes the stable branch of solutions while the dotted [red]{} line denotes the unstable branch of solutions. In (a) the stationary density $e^{-x^2}/Z$ changes stability at the critical interaction energy $\kappa_2 = \kappa_{2\sharp}= -2.4$ and the new stable density is asymmetric adhering to one wall (Figure \[fig:bif\_fig\_right\]). In (b), in the absence of a confining potential, the uniform density becomes unstable at the critical interaction energy $\kappa_2 = \kappa_{2\sharp} = 0.4$ and the density may become multi-modal (Figure \[fig:bif\_stable\_equilibria\_V2\_cosine\]).[]{data-label="fig:bifurcation_diagram"}](bifurcation_diagram_for_Phi2_cosine.pdf "fig:"){width="\textwidth"}

Description of Main Results and Organisation of the Paper.
----------------------------------------------------------

### Main Results {#main-results .unnumbered}

The main results of this work are threefold.

1.  We establish existence and uniqueness of weak solutions to DDFTs including two-body HI governed by equations , - with no-flux boundary conditions.

2.  We derive *a priori* convergence estimates of the density $\varrho(\vec{r},t)$ to equilibrium in $L^2$ and relative entropy.

3.  We study the stability of equilibrium states and construct bifurcation diagrams for two numerical applications.

These results are of particular interest for physical applications of colloidal systems where conservation of mass is either a desirable or necessary property of the system. Additionally, the stability theorem contrasts with simpler linear stability analyses of similar systems of gradient flow structure with periodic boundary conditions [@martzel2001mean], [@greg_mckean_vlasov_torus] which may be tackled by means of Fourier analysis.

### Organisation of the Paper {#organisation-of-the-paper .unnumbered}

The paper is organised as follows: in Section \[sec:preliminaries\] we present the boundary and initial conditions, introduce the main notation, nondimensionalise the main equations, state the stationary equation for the density, define the weak formulation of the Smoluchowski equation including full HI and provide a list of assumptions. In Section \[sec:statement\_of\_main\_results\] we state the main results of the present work in a precise manner. In Section \[sec:ex\_uni\_full\_HI\] we provide an existence and uniqueness theorem for the flux $\vec{a}$ when full HI are included. In Section \[sec:char\_stationary\_sol\] we characterise solutions of the stationary problem and convergence to equilibrium in $L^2$ as $t \to \infty$. In Section \[sec:global\_asymptotic\_stability\] we obtain results on the global asymptotic stability of the stationary densities by showing that the free energy is a continuous functional for all two-body interaction strengths. Additionally we prove an H- theorem for the equilibria, provide *a priori* convergence estimates in relative entropy, derive an asymptotic expansion of the equilibria for small interaction energy and perform a spectral analysis of the linearised nonlocal Smoluchowski operator. In Section \[sec:bifurcation\_theory\] we provide necessary and sufficient conditions for phase transitions in generalised DDFT-like systems with no-flux boundary conditions. In Section \[sec:manufactured\_bif\] we construct the bifurcation diagram for some example problems. In Section \[sec:existence\_uniqueness\_with\_partial\_HI\] we obtain an existence and uniqueness theorem for the Smoluchowski equation with non-constant diffusion tensor and effective drift vector dependent on the two-body HI tensors $\bm{Z}_1$ and $\bm{Z}_2$. In Section \[sec:discussion\] we present our concluding remarks and state some open problems. In Appendix \[sec:classical\_paraboliv\_pde\] we provide some technical results that are used in the proof of Theorem \[thm:exis\_uniq\_weak\_sol\_rho\]. Finally in Appendix \[app:nomenclature\] we provide a list of nomenclature.

Preliminaries {#sec:preliminaries}
=============

In this section we specify the nonlinear boundary conditions and initial data for the DDFT . We also nondimensionalise the governing equations and provide the assumptions on the regularity of the potentials, correlation function, diffusion tensor and initial data.

Boundary Conditions. {#subsec:boundary_conditions}
--------------------

When $U = \mathbb{R}^d$ we take $$\begin{aligned}
\label{bc:density_and_flux_decaying}
\left\{\begin{aligned}
\varrho(\vec{r},t)\to 0 \\
\vec{a}(\vec{r},t)\to \vec{0}
\end{aligned}\right. \quad \text{ as } \quad |\vec{r}|\to \infty,\end{aligned}$$ where we require $V_1$ to be growing at least quadratically as $\vec{r}\to \infty$. Physically-speaking this prevents the density from running out to infinity. When $U\subset \mathbb{R}^d$ is open and bounded we impose that the total mass of the system $M$ remains constant, in particular we have $$\begin{aligned}
\label{bc:mass_preserving}
\vec{a}(\vec{r},t)\cdot\vec{n}\bigg|_{\partial U\times [0,T]} = 0.\end{aligned}$$ The boundary condition may be viewed as a [*nonlinear*]{} Robin condition imposing the flux through the boundary $\partial U$ is zero for all time $t\in [0,T]$. If $\varrho$ is a number density then $\int \mathrm{d}\vec{r}\, \varrho = N$ for all time, however for the analysis in Section \[sec:ex\_uni\_full\_HI\] and onwards we will assume $\varrho$ is a probability density so that $\int \mathrm{d}\vec{r}\, \varrho = 1$. The rescaling between number and probability densities is discussed in the following section.

Initial Conditions.
-------------------

We will assume that the initial data has finite free energy and is consistent with the imposed boundary conditions. For example, one could prescribe initial data $(\varrho_0,\vec{a}_0)^\top$ such that $$\begin{aligned}
\label{eq:initial_data_for_a_and_rho}
\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho_0](\vec{r}) = \mu_{c},\qquad
\vec{a}_0 = \vec{0}.\end{aligned}$$ where $\mu_{c}$ is the chemical potential, constant at equilibrium. It is straightforward to check that $(\varrho_0,\vec{a}_0)^\top$ is an equilibrium point of the system . Commonly, one then drives the system out of equilibrium via a time-dependent external potential. In principle $\mu_{c}$ may be given and the equations , are well defined. In practice, for complicated particle configurations, $\mu_{c}$ is not known but can be computed by minimising the free energy along with the additional constraint $\int_U\mathrm{d}\vec{r}\, \varrho_0(\vec{r})=N$, where $N$ is the (expected) number of particles for a finite system and $\varrho_0$ is a number density. Note that $\mu_{c}$ is a potential, so by raising it one may force more particles into the system. We will assume that $\mu_c$ is constant to fix the number of particles. To ensure $\varrho$ (and $\varrho_0$) is a probability density one may rescale $\varrho/N = \tilde{\varrho}$, $N g = \tilde{g}$ and $\vec{a}/N^2 = \tilde{\vec{a}}$, where the tilde denotes the new variable, so that $\int_U\mathrm{d}\vec{r}\, \varrho_0(\vec{r})=1$ and equations - become independent of $N$.

This provides a method of converting back to the number density which is typically used in numerical modelling of finite colloidal systems [@goddard2016dynamical], [@goddard2012unification], [@goddard2012general]. Throughout however, since we will frequently use the integral of the density, we will assume $\varrho$ and $\varrho_0$ are probability densities to ease notation. With this, one has three equations for three unknowns $\mu_{c}$, $\varrho_0$, $\vec{a}_0$ and the initial density $\varrho_0$ can be computed. For the rest of paper it is convenient to work in dimensionless units. We now nondimensionalise the governing equations.

Evolution Equations.
--------------------

We now nondimensionalise our equations. Let $\mathrm{L}$, $\tau$, $\text{U}$ be characteristic length, time and velocity scales respectively, then by nondimensionalising $$\begin{aligned}
\vec{r}\sim \mathrm{L} \tilde{\vec{r}},\quad t\sim\tau\tilde{t}, \quad  \mathrm{U} = \tfrac{\mathrm{L}}{\tau}, \quad \varrho\sim\tfrac{1}{\mathrm{L}^d}\tilde{\varrho},\quad \mathcal{F}\sim k_BT \tilde{\mathcal{F}},\quad \vec{a}\sim \mathrm{A}\tilde{ \vec{a}}.\end{aligned}$$ where $d$ is the physical dimension and $\mathrm{A}$ is a characteristic flux scale. The system becomes (after dropping tildes) $$\begin{aligned}
\partial_t\varrho(\vec{r},t) = -\tfrac{1}{Fr}\times \tfrac{\tau^{-1}}{\gamma}\times \mathrm{A}\times \mathrm{L}^{d+1} \nabla_{\vec{r}}\cdot \vec{a}(\vec{r},t),\end{aligned}$$ where we have defined the Froude number $Fr = m\mathrm{U}^2/(k_BT)$. By choosing $Fr = 1$, $\tau = \gamma^{-1}$ and $\mathrm{A} = 1/\mathrm{L}^d$ we simplify the system of equations to the following boundary value problem.\

\[prop:non\_dimensional\_time\_evolving\_flux\_eqn\]

The non-dimensional one-body density $\varrho(\vec{r},t)$ and flux $\vec{a}(\vec{r},t)$ evolve according the the boundary value problem $$\begin{aligned}
\label{eq:evolution_eqn_for_a_dimensionless}
\begin{cases}
&\qquad\qquad\partial_t\varrho = -\nabla_{\vec{r}}\cdot \vec{a}(\vec{r},t), \\
&\vec{a}(\vec{r},t) + \bm{H}[\vec{a},\varrho] +\varrho(\vec{r},t)\bm{D}(\vec{r},[\varrho],t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]=0,
\\
&\qquad [\bm{H}[\vec{a},\varrho] +\varrho(\vec{r},t)\bm{D}(\vec{r},[\varrho],t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]]\cdot\vec{n}\big|_{\partial U}=0.
\end{cases}\end{aligned}$$

We note that when the off-diagonal HI tensor $\bm{Z}_{2}=0$, by using the definitions of $\mathcal{F}$ and $\bm{D}$ , the evolution equations in may be written as a nonlinear Smoluchowski equation (such as ) with non-constant diffusion coefficient. However we observe that even when $\bm{Z}_2 \neq 0$ the dynamics may be recast into a Smoluchowski equation for $\varrho$ under an effective drift vector dependent on $\bm{Z}_2$.\

\[prop:non\_dimensional\_time\_evolving\_rho\_eqn\]

The non-dimensional one-body density $\varrho(\vec{r},t)$ evolves according the the boundary value problem $$\begin{aligned}
\label{eq:evolution_eqn_for_rho_dimensionless}
\begin{cases}
&\qquad\partial_t\varrho=\,\nabla_{\vec{r}}\cdot\left[Pe^{-1}\bm{D}\nabla_{\vec{r}}\varrho+\varrho\,\bm{D}\left(\nabla_{\vec{r}} (\kappa _1V_1+\kappa _2\,(gV_2)\star\varrho) + \bm{A}[\bm{a}]\right)\right],
\\
&\qquad \qquad \qquad \qquad \qquad  \Pi [\varrho]\cdot\vec{n}\big|_{\partial U} = 0,\\
&\qquad\Pi[\varrho]:= \bm{D}\,\left(\nabla_{\vec{r}}\varrho + 
\varrho\, \nabla_{\vec{r}}(\kappa_1 V_1(\vec{r},t)+ \kappa_2 (gV_2)\star \varrho)+\bm{A}[\bm{a}]\right),
\end{cases}\end{aligned}$$ where $\bm{A}[\bm{a}]$ is an effective background flow induced by the hydrodynamic interactions defined by $$\begin{aligned}
\label{eq:def_of_V1_eff}
\bm{A}[\vec{a}]:=\int_U\mathrm{d}\vec{r}'\, g(\vec{r},\vec{r}') \bm{Z}_2(\vec{r},\vec{r}')\vec{a}(\vec{r'},t),\end{aligned}$$ $\kappa_1$, $\kappa_2$ are non-dimensional constants measuring the strength of confining and interaction potentials respectively, $Pe=\mathrm{L}\text{U}/\alpha$ is the P[é]{}clet number measuring the ratio of advection rates to diffusive rates and $\alpha = k_B\mathrm{T}/(m\gamma)$.

Corollary \[prop:non\_dimensional\_time\_evolving\_rho\_eqn\] is the general formulation of the nondimensional equations - when $\bm{Z}_2\neq0$, including a non-constant diffusion coefficient and an effective drift. Throughout this paper, to study the intermediate regime of equally strong advection and diffusion, we set $Pe = 1$. Additionally, we redefine the two-body potential to absorb the correlation function $g$ to ease notation, $V_2(\vec{r},\vec{r}'):= g(\vec{r},\vec{r}') V_2(\vec{r},\vec{r}')$. In practice, there are many choices for $g$, for example the hard sphere approximation takes $g(|\vec{r}- \vec{r}'|) = 0$ for $|\vec{r}-\vec{r}'| < 1$ and unity otherwise. Alternatively $g$ may be obtained numerically from microscopic dynamics. We consolidate the choices for equations , in Section \[subsec:assumptions\_definitions\].

The effective drift $\bm{A}[\vec{a}]$, dependent on $\bm{Z}_2$ and $\vec{a}$ may be determined once $\vec{a}(\vec{r},t)$ is solved from the second equation in . Note that the evolution equation in may be viewed as a generalised McKean-Vlasov equation with a non-constant diffusion tensor and confining potential. In particular the McKean-Vlasov equation may be recovered in the special case $\bm{Z}_1 = \bm{Z}_2 = V_1 = 0$, see for example [@greg_mckean_vlasov_torus], [@chayes2010mckean]. We will use Corollary \[prop:non\_dimensional\_time\_evolving\_rho\_eqn\], to write the full dynamics including full HI, to obtain our results on weak solutions for $\varrho(\vec{r},t)$ (see Theorem \[thm:eigenfn\_expansion\_of\_flux\], Section \[sec:ex\_uni\_full\_HI\] and Theorem \[thm:existence\_and\_uniqueness\], Section \[sec:existence\_uniqueness\_with\_partial\_HI\]). We continue to the next section by stating the stationary boundary value problem for equilibrium states $\varrho(\vec{r})$.

Stationary Equations.
---------------------

For general $\bm{Z}_2$ we will show in Theorem \[thm:association \_of\_free\_energy\] that the stationary density $\varrho(\vec{r})$ satisfies $$\begin{aligned}
\label{eq:stationary_eqn_for_rho_dimensionless}
\begin{cases}
&\qquad 0=\,\nabla_{\vec{r}}\cdot[\bm{D}\nabla_{\vec{r}}\varrho+\varrho\,\bm{D}\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2\,V_2\star\varrho)],
\\
&\qquad \qquad \qquad \qquad \qquad  \Pi [\varrho]\cdot\vec{n}\big|_{\partial U} = 0,\\
&\qquad\Pi[\varrho]:= \bm{D}\,(\nabla_{\vec{r}}\varrho + 
\varrho\, \nabla_{\vec{r}}(\kappa_1 V_1(\vec{r},t)+ \kappa_2 \,V_2\star \varrho)).
\end{cases}\end{aligned}$$ We now discuss regularity on the potentials and diffusion tensor.

Assumptions & Definitions. {#subsec:assumptions_definitions}
--------------------------

Typically for long range HI the $\boldsymbol{Z}_i$ exhibit singularities at the origin (particle centres) so the correlation function $g$ is a necessary inclusion and provides a way of smoothing $\boldsymbol{D}$ and we assume $g \in L^\infty(U)$. For $\varrho \geq 0$ the diffusion tensor $\boldsymbol{D}$ as a convolution with the density will then be a weakly differentiable function. For the existence and uniqueness theory in Appendix \[sec:classical\_paraboliv\_pde\] and Section \[sec:existence\_uniqueness\_with\_partial\_HI\] we require that first derivatives of $\bm{D}_{ij}$ to be bounded in $L^\infty(U)$ so that all coeeffiecients of the PDE are uniformly bounded.

Out of equilibirum, we will suppress the time dependence on $\bm{D}$, $V_1$ simply to ease notation. However at equilibrium $\bm{D}$, $V_1$ are assumed to be independent of time, indeed in order for equilibrium states of the density and flux to be well defined. We note that is $\bm{D}$ positive definite and symmetric, as it has been rigorously shown to be [@goddard2012overdamped]. In summary we have the following notational choices and assumptions for the evolution problem .

#### **[Notation]{}**

Throughout we ease notation on the two-body interaction potential.

-   The two-body interaction potential is redefined to absorb the correlation function $g$ $$\begin{aligned}
    \label{ass:V2_redef}
    V_2\stackrel{\text{redef}}{:=} g V_2. \tag{N1}\end{aligned}$$

For the dynamics we assume:

#### **[Assumptions D]{}**

-   The diffusion tensor $\bm{D}$ is symmetric, positive definite, and the first derivatives of $\bm{D}_{ij}$ are bounded in $L^\infty(U)$ $$\begin{aligned}
    \label{ass:D_pos_def_weak_diffable}
    \bm{D}_{ij}\in W^{1,\infty}(U). \tag{D1}\end{aligned}$$

-   The diagonal and off-diagonal blocks of the HI tensors are uniformly bounded in the sense $$\begin{aligned}
    \label{ass:Z2_uniformly_bd}
    \|gZ_2\|_{L^\infty(U)}<\infty, \quad  \|gZ_1\|_{L^\infty(U)}<\infty\tag{D2}\end{aligned}$$

-   The initial data $\varrho_0$ is a non-negative, square-integrable, absolutely continuous probability density $$\begin{aligned}
    \label{ass:rho_in_L2_P_ac}
    \varrho_0\in P_{ac}(U)\cap L^2(U). \tag{D3}\end{aligned}$$

-   The potentials each have two bounded derivatives $$\begin{aligned}
    \label{ass:V1_V2_in_W_1_inf}
     V_1,V_2\in W^{2,\infty}(U). \tag{D4}\end{aligned}$$ The functions $V_1$ and $V_2$ are the confining and two-body interaction potentials respectively, the former having explicit time dependence ($V_1 = V_1(\vec{r},t)$) only when we intend to drive - and out of equilibrium, and $V_1 = V_1(\vec{r})$ when we are concerned with the equilibrium properties of - and . This distinction will be important for the H Theorem and equilibrium theory in Section \[sec:global\_asymptotic\_stability\].

For the equilibrium problem we will assume:

#### **[Assumptions E]{}**

-   The potentials have first order weak derivatives in $L^2(U)$ $$\begin{aligned}
    \label{ass:V1_V2_in_H_1}
     V_1,V_2\in H^{1}(U). \tag{E1}\end{aligned}$$

In particular, Assumption will permit us to establish smooth stationary densities. Note that typical inter-particle potentials, such as Morse or Coulomb, are unbounded as the particle separation goes to zero. This is once again mitigated by the choice of $g$, which we recall has been absorbed into $V_2$ by assumption . In general we admit non-convex $V_1$ and $V_2$, for example multi-well potentials, except for in the convergence result of Theorem \[thm:rel\_entropy\_convergence\] where $V_1$ must be convex in order to invoke a log-Sobolev inequality on the measure $\mu$ given by .

The assumption that $\varrho_0\in P_{\text{ac}}(U)$ is included in order to cover a wider set of physically relevant scenarios. In particular we permit initial data such that $\varrho_0|_{A}=0$ for some $A\subset U$ where $A$ is non-empty. Physically speaking this system could correspond to, at time $t=0$, a box partitioned into closed regions with at least one region containing no particles. Then, instantaneously as soon as $t>0$, the partition is removed allowing the particles to move freely. At the end of Section \[sec:existence\_uniqueness\_with\_partial\_HI\] we will show by simple application of Harnack’s inequality, that we obtain strictly positive densities $\varrho(\vec{r},t_1)>0$ after an arbitrarily small time $t_1>0$. Principally this is provided by the property , since $\bm{D}$ is positive definite, the diffusion of density in the system is everywhere propagating in $U$.

Additionally, by the positive definite property in we may uniquely define the square root of $\bm{D}$ denoted $\bm{D}^{1/2}$ such that $$\bm{D}^{1/2}\bm{D}^{1/2}=\bm{D}$$ for every $\vec{r}\in U$, $t\in [0,T]$ and each $\varrho$.

We also define the eigenvalues $\mu_i\in \mathbb{R}^+$ of $\bm{D}$ and eigenvectors $\vec{e}_i(\vec{r},[\varrho],t)\in L^2(U)$ such that $$\begin{aligned}
\label{eq:mu_i_eigenvalue_eqn_D}
\bm{D}\vec{e}_i = \mu_i\vec{e}_i.\end{aligned}$$ Note that $\{\vec{e}_i\}_{i=1}^d$ forms an orthonormal basis of $\mathbb{R}^d$ (since $\bm{D}$ is a bounded, symmetric operator) for $i = 1,\cdots, d$ such that $$\begin{aligned}
\label{eq:orthogonality_of_ei}
\langle\vec{e}_i(\vec{r},[\varrho],t),\vec{e}_j(\vec{r},[\varrho],t)\rangle = \int_U\mathrm{d}\vec{r}\, \vec{e}_i(\vec{r},[\varrho],t)\cdot \vec{e}_j(\vec{r},[\varrho],t) = \delta_{ij}.\end{aligned}$$ We continue to the next section by defining a weak formulation of the dynamics .\

Weak Formulation.
-----------------

We provide the weak formulation of the full dynamics including HI for $\bm{Z}_1, \bm{Z}_2$ not necessarily zero.\

Let $\vec{a}(\vec{r},t)$ be a given flux. We say $\varrho\in L^2([0,T]; H^1(U))\cap L^\infty([0,T]; L^2(U))$ and $\partial_t\varrho\in L^2([0,T]; H^{-1}(U))$ is a weak solution to if for every $\eta\in L^2([0,T]; H^1(U))$ $$\begin{aligned}
\label{eq:weak_formulation_including_V1eff}
\int_0^T\mathrm{d}t\, \langle \partial_t \varrho(t), \, \eta(t)  \rangle +\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, \nabla_{\vec{r}}\eta\cdot \bm{D}\,[\nabla_{\vec{r}}\varrho
 +\varrho\,\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2\,V_2\star\varrho + \bm{A}[\vec{a}])]=0\end{aligned}$$ where $\varrho_0 = \varrho(0)$. Here, $\bm{A}[\vec{a}]$ is the effective drift induced by $\bm{Z}_2$ and is defined by equation .

It will be shown in the following sections (in particular Corollary \[cor:a\_is\_zero\_at\_equilibrium\]) that $\bm{A}[\vec{a}] \to \vec{0}$ as $t \to \infty$. We now state our main results in a precise manner.\

Statement of Main Results {#sec:statement_of_main_results}
=========================

Our main results concern existence, uniqueness and convergence to equilibrium of the density of colloids $\varrho$ and flux $\vec{a}$ on $U$ a compact subset of $\mathbb{R}^d$. The first result concerns existence of the flux $\vec{a}(\vec{r},t)$ with non-zero hydrodynamic interactions, the convergence of $\vec{a}(\vec{r},t)$ to zero at equilibrium and existence and uniqueness of fixed points of .\

\[them:exis\_uni\_of\_flux\]

Let $\bm{Z}_1,\bm{Z}_2$ be real, symmetric and $\mu_{\max}\|g\bm{Z}_2\|_{L^\infty(U)}<1$. Then\

1.  There exists a unique $\vec{a}(\vec{r},t)\in L^2(U)$ solving the evolution equation for each $\varrho(\vec{r},t)$. In particular $$\begin{aligned}
    \vec{a}(\vec{r},t) = \sum_{n=1}^d\delta_n\sum_{i=1}^d\frac{\psi_i}{\phi_n-\mu_i^{-1}} \vec{e}_i(\vec{r},[\varrho],t)\end{aligned}$$ where $\vec{e}_i(\vec{r},[\varrho],t)$ are eigenvectors of the diffusion tensor $\bm{D}(\vec{r},[\varrho],t)$ and $\delta_n$, $\phi_n$, $\psi_i$, $\mu_i^{-1}\in \mathbb{R}$.

2.  In addition, every stationary density $\varrho(\vec{r})$ and stationary flux $\vec{a}(\vec{r})$ are independent of the HI tensors and satisfy $$\begin{aligned}
    \varrho(\vec{r})\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho(\vec{r})] = \vec{0}, \quad \vec{a}(\vec{r}) = \vec{0}\end{aligned}$$ and consequently $\varrho(\vec{r})$ minimises the free energy $\mathcal{F}[\varrho](r)-\mu_{c} \int_U\mathrm{d}\vec{r}\,\varrho$, where $\mu_{c}$ is the chemical potential.

3.  If, in addition, $|\kappa_2|< \|V_2\|_{L^\infty(U)}^{-1}$ then $(\vec{a}_\star(\vec{r}),\varrho_\star) = (\vec{0}, \varrho_\infty)$ are the unique fixed points of and $\varrho_\infty(\vec{r})$ is given by the self-consistency equation $$\begin{aligned}
    \varrho_\infty(\vec{r}) = \frac{e^{-(\kappa_1 V_1(\vec{r})+\kappa_2 V_2\star \varrho_\infty)}}{Z(\varrho_\infty)}\end{aligned}$$

for $Z(\varrho_\infty) = \int_U\mathrm{d}\vec{r}\,e^{-(\kappa_1 V_1(\vec{r})+\kappa_2 V_2\star \varrho_\infty)}$.

For the evolution system we present the following second main result of the paper.\

\[thm:exis\_uniq\_weak\_sol\_rho\]

Let $\bm{Z}_1,\bm{Z}_2$ be real, symmetric and $\mu_{\max}\|g\bm{Z}_2\|_{L^{\infty}(U)}<1$, where $\mu_{\max}$ is the largest eigenvalue of $\bm{D}$, with $\varrho_0 \in C^\infty(U)$, $\varrho\geq 0$ and $\int_U \mathrm{d}\vec{r} \varrho_0(\vec{r})=1$. Then there exists a unique weak solution $\varrho\in L^\infty([0,T];L^2(U))\cap L^2([0,T]; H^1(U))$, with $\partial_t\varrho\in L^2([0,T]; H^{-1}(U))$ for , in the sense , and the following energy estimate holds $$\begin{aligned}
\|\varrho\|_{L^{\infty}([0,T];L^2(U))}+\|\varrho\|_{L^{2}([0,T];H^1(U))}+\|\partial_t\varrho\|_{L^2([0,T]; H^{-1}(U))}\leq C(T) \|\varrho_{0}\|_
{L^2(U)},\end{aligned}$$ where $C(T)$ is a constant dependent on $T$, $U$ and $\mu_{\max}$.

The existence and uniqueness is proved in Theorem \[thm:existence\_and\_uniqueness\], whilst the bound is shown in Lemma \[lem:weak\_conv\_results\].

Furthermore, we prove existence and uniqueness of the stationary density, and exponentially fast convergence in relative entropy.\

Let $\bm{Z}_1,\bm{Z}_2$ be real, symmetric and $\varrho$ be a solution to the DDFT with smooth initial data and smooth $V_1$, $V_2$. Then there exists stationary density $\varrho(\vec{r},t) = \varrho_0(\vec{r})$. If $|\kappa_2| \leq 1/4\times \|V_2\|_{L^\infty}^{-1}$ then the stationary solution is unique and is denoted by $\varrho_\infty$.

The proof of this result is standard, see [@dressler1987stationary].

The third main result of this paper concerns *a priori* estimates for exponential convergence of the density to stationarity.\

Let $\bm{Z}_1,\bm{Z}_2$ be real, symmetric and $\varrho$ be a solution to the DDFT with smooth initial data and smooth $V_1$, $V_2$. If $\kappa_2 \leq 1/4\times \|V_2\|_{L^\infty}^{-1}$ then\

1.  [**[Convergence in ]{}**]{}$\bm{L^2(U):}$ For $\kappa_1=0$ (in the absence of a confining potential) if $$\begin{aligned}
    \kappa_2^2<\frac{\mu_{\min}c^{-2}_{pw}\|\nabla_{\vec{r}}V_2\|^{-2}_{L^\infty(U)}}{2(1+e)\mu_{\max}},\end{aligned}$$ where $\mu_{\min}$ and $\mu_{\max}$ are the smallest and largest eigenvalues of the diffusion tensor $\bm{D}$, then $\varrho\to \varrho_\infty$ in $L^2(U)$ exponentially fast as $t\to \infty$. For $\kappa_1\neq 0$ the convergence criteria is modified to $$\begin{aligned}
    \mu_{\max}(\kappa_1^2\|\nabla_{\vec{r}}V_1\|_{L^\infty(U)}^2+2\kappa_2^2(1+e)\|\nabla_{\vec{r}}V_2\|_{L^\infty(U)}^2)<\frac{\mu_{\min}}{c_{pw}^2}.\end{aligned}$$

2.  [**[Convergence in Relative Entropy:]{}**]{} For any fixed confining potential $V_1$ such that the measure $\mu'(\mathrm{d}\vec{r}) = \mathrm{d}\vec{r}\,e^{-\kappa_1V_1}/Z$ satisfies a log-Sobolev inequality and provided $$\begin{aligned}
    \kappa_2^2<\frac{c^{-1}_{ls}}{2\|\nabla V_2\|^2_{L^{\infty(U)}}}\end{aligned}$$ then the measure $\mu$ in satisfies a log-Sobolev inequality and $\mathcal{H}(\varrho |\varrho_\infty) \to 0$ exponentially fast as $t\to \infty$ where $$\begin{aligned}
    \mathcal{H}(\varrho |\varrho_\infty) = \int_U\mathrm{d}\vec{r}\, \varrho \log \left(\tfrac{\varrho}{\varrho_\infty}\right)\end{aligned}$$ denotes the relative entropy.

For part 1, see Theorem \[lem:exp\_conv\_in\_L2\] and Proposition \[prop:general\_kappa1\_convergence\]. Theorem \[thm:rel\_entropy\_convergence\] gives the result for part 2.

The log-Sobolev inequality for $\mu$ is established by the Holley-Stroock perturbation lemma [@holley1986logarithmic]. The constants $c_{pw}$, $c_{ls}$ are the Poincar[é]{}-Wirtinger and log-Sobolev constants respectively. Nowhere do we assume parity on the two-body potential nor $V_2$ have zero mean. Additionally the optimal $c_{pw}$ is the inverse square root of the smallest eigenvalue of the Laplacian on the domain $U$ with no-flux boundary conditions.

We have the following conditions for the existence of bifurcating branches of steady states $\varrho(\vec{r})$.\

Fix $\kappa_1$ and let $\kappa_2\in(-\infty,\infty)$. Let $\mathcal{L}_{1} = \mathcal{A}_{\kappa_2}+\kappa_2\mathcal{B}$ denote the linearised operator to the stationary problem with eigenvalues $\lambda(\kappa_2)$ and eigenfunctions $w^{(\kappa_2)}(\vec{r})$. Denote by $\mathcal{A}_{\kappa_2}$, $\mathcal{B}$ the local and nonlocal parts of $\mathcal{L}_{1}$ respectively. Denote by $\gamma_k^{(\kappa_2)}$ the eigenvalues of $\mathcal{A}_{\kappa_2}$ with eigenvectors $v^{(\kappa_2)}$. If the solution $\kappa_2^\star(\lambda)$ of the equation $\lambda = \lambda_{k^\star}(\kappa_2^\star) $ exists, then it is unique and is given by the nonlinear equation $$\begin{aligned}
\kappa_2^\star(\lambda) = \left(\sum_{i=1}^\infty\tfrac{\theta_i^{(\kappa_2)}\gamma
_i^{(\kappa_2)}\beta_i^{(\kappa_2)}}{\lambda-\gamma_i^{(\kappa_2)}}\right)^{-1}. \end{aligned}$$

As a corollary we can determine the necessary condition on the interaction strength for a bifurcation of stable equilibrium densities solving .\

Provided that the spectral gap of $\mathcal{A}_{\kappa_2}$ is sufficiently large, that is, $$\begin{aligned}
|\kappa_2|<\frac{\min_{i,j\in\mathbb{N}}|\gamma_i^{(\kappa_2)}-\gamma_{j}^{(\kappa_2)}|}{2\|\mathcal{B}\|}\end{aligned}$$ then $\lambda(\kappa_2)\in\mathbb{R}$ and the point of critical stability $\kappa_{2_\sharp}$ occurs at the solution of the nonlinear equation $$\begin{aligned}
\kappa_{2_{\sharp}} = -\left(\sum_{i=1}^\infty\theta_i^{(\kappa_2)}\beta_i^{(\kappa_2)}\right)^{-1},\end{aligned}$$ where $\theta_i\beta_i$ are coefficients of the two-body potential expanded in the orthonormal basis of eigenvectors $\{v_k^{(\kappa_{2_\sharp})}\}_{k = 1}^\infty$.

The proofs of these results are given by Theorem \[thm:epsilon\_as \_a\_fn\_of\_lambda\] and the discussion immediately following it. We also obtain the following theorem for existence of bifurcations for the stationary equation .\

Let $\kappa_2\in (-\infty,\infty)$and let $\{\beta_{n}^{-1}\}_{n=1}^{\infty}$ be the eigenvalues of $\mathcal{R}$ with eigenfunctions $\{u_n\}_{n=1}^{\infty}$ where $$\begin{aligned}
\mathcal{R}[u_n] = -\varrho_{\kappa_2}(\vec{r})\int_U\mathrm{d}\vec{r}'\, V_2(\vec{r},\vec{r}')u_n(\vec{r}')\end{aligned}$$ and $\varrho_{\kappa_2}$ is a stationary solution to . If $|\kappa_2|\geq |\beta_1|$ then $\varrho_{\kappa_2}$ is unstable with respect to $\{u_n\}_{n=1}^{\infty}$ with $(\beta_1,u_1)$ a bifurcation point of where $\beta_1$ is the smallest eigenvalue of $\mathcal{R}^{-1}$ and $w_1$ is the eigenfunction of $\mathcal{R}$ associated to $\beta_1^{-1}$. There exists $\varrho_\ast>0$ such that $\mathcal{F}[\varrho_\ast]<\mathcal{F}[\varrho_{\kappa_2}]$.

We now give our arguments for Theorem \[them:exis\_uni\_of\_flux\].\

Existence & Uniqueness of Flux With Full HI {#sec:ex_uni_full_HI}
===========================================

We return to the full formulation of the overdamped DDFT with HI. The contraction condition $\mu_{\max}\|g\bm{Z}_2\|_{L^{\infty}(U)}<1$ can be seen as a necessary condition on the invertibility of an operator closely related to the positive definite grand friction tensor $\boldsymbol{\Gamma}$. Note that the flux equation in may be written more generally as $$\begin{aligned}
\label{eq:matrix_operator_friction_eqn}
(\bm{1}+ \mathcal{Z}_1^\varrho + \mathcal{Z}_2^\varrho)[\vec{a}(\vec{r},t)]=-\varrho(\vec{r},t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\end{aligned}$$ where the actions of the integral operators $\bm{1}+ \mathcal{Z}_1^\varrho$ and $\mathcal{Z}_2^\varrho$ are defined by

$$\begin{aligned}
(\bm{1}+ \mathcal{Z}_1^\varrho)[\vec{a}_1]&=\vec{a}_1(\vec{r},t)+\int_U\mathrm{d}\vec{r}' g(\vec{r},\vec{r}')\varrho(\vec{r}',t)\boldsymbol{Z}_1(\vec{r},\vec{r}')\times \vec{a}_1(\vec{r},t),\label{eq:friction_like_system_1plusZ1}\\
\mathcal{Z}_2^\varrho[\vec{a}_2]&=\varrho(\vec{r},t)\int_U\mathrm{d}\vec{r}'\, g(\vec{r},\vec{r}')\bm{Z}_2(\vec{r},\vec{r}')\vec{a}_2(\vec{r}',t).\label{eq:friction_like_system_Z2}\end{aligned}$$

Notice how the integral-matrix operators in - resemble the operators in the first row of the grand resistance matrix for a two particle system from classical hydrodynamics [@happel2012low]. The following lemma establishes a solvability condition for the flux equation in equation .

\[thm:cond\_converg\_fred\_det\]

Let $\bm{1}+ \mathcal{Z}_1^\varrho$ and $\mathcal{Z}_2^\varrho$ be bounded linear operators. Suppose $\mathcal{A}_\varrho:=(\bm{1}+ \mathcal{Z}_1^\varrho)^{-1}\mathcal{Z}_2^\varrho$ is compact in $L^2(U,\varrho^{-1}(\vec{r},t))$. If $\mu_{\max}\|g\bm{Z}_2\|_{L^\infty(U)}<1$ then the matrix integral operator $\bm{1}+ \mathcal{Z}_1^\varrho  + \mathcal{Z}_2^\varrho $ is invertible and the system is well-posed.

Since $\bm{1}+ \mathcal{Z}_1^\varrho$ is positive definite it is invertible, therefore may be rewritten $$\begin{aligned}
\label{eq:matrix_integral_op_eqn_rewrite}
(\bm{1}  + (\bm{1}+ \mathcal{Z}_1^\varrho )^{-1}\mathcal{Z}_2^\varrho) [\vec{a}(\vec{r},t)] = -(\bm{1}+ \mathcal{Z}_1^\varrho)^{-1} \varrho(\vec{r},t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho].\end{aligned}$$ We note that $\mathcal{A}^\varrho$ is a trace-class operator and that the left hand side of is an operator of the form $\bm{1}-\lambda\mathcal{A}_\varrho$. By classical theory [@fredholm1900nouvelle], [@lax2014functional] we have the identity $$\begin{aligned}
\label{eq:det(1+A)_in_terms_of_exp}
\det(\bm{1}-\lambda\mathcal{A}_\varrho) =  \exp\Big\{-\sum_{n=1}^\infty\tfrac{\text{Tr}(\mathcal{A}_\varrho^n)}{n}\lambda^n \Big\}.\end{aligned}$$

When $\lambda = -1$ we recover the determinant for the Fredholm operator on the left hand side of . Particularly, since for our consideration $|\lambda|=1$, the convergence of the infinite summation inside the argument of the exponential in will depend on the size of $\text{Tr}\mathcal{A}_\varrho$, and when $\lambda = -1$, the summand is an alternating sequence so we demand absolute convergence for the sum in to converge. We obtain results in $L^2(U, \varrho^{-1})$.

By definition of the trace we have $$\begin{aligned}
\text{Tr}\mathcal{A}^n_\varrho = \sum_{k=1}^d\langle \mathcal{A}^n_\varrho \vec{e}_k(\vec{r}), \vec{e}_k(\vec{r})\rangle_{L^2(U,\varrho^{-1})},\end{aligned}$$ where $\{\vec{e}_k\}_{k}$ are vectors such that their components form an orthonormal basis of $L^2(U,\varrho^{-1})$, in particular we choose the eigenvectors of the diffusion tensor $\bm{D}$. Since $\mathcal{A}$ is an integro-matrix operator the inner product is given by, for $n=1$ $$\begin{aligned}
\text{Tr}\mathcal{A}_\varrho &= \sum_{k=1}^d\langle \mathcal{A}_\varrho \vec{e}_k(\vec{r}), \vec{e}_k(\vec{r})\rangle_{L^2(U,\varrho^{-1})}\\
&=\sum_{k=1}^d\int_U\mathrm{d}\vec{r}\,\vec{e}_{k}(\vec{r})\cdot\bm{D}(\vec{r})\int_U\mathrm{d}\vec{r}'\,g(\vec{r},\vec{r}')\bm{Z}_2(\vec{r},\vec{r}')\vec{e}_{k}(\vec{r}')\\
&\leq \sum_{k=1}^d\mu_{\max}\|g \bm{Z}_2\|_{L^\infty(U)}\int_U\mathrm{d}\vec{r}\,\vec{e}_{k}(\vec{r})\cdot\int_U\mathrm{d}\vec{r}'\,\vec{e}_{k}(\vec{r}')\leq d\mu_{\max}\|g \bm{Z}_2\|_{L^\infty(U)}|U|\end{aligned}$$ where we have used $\int_U\mathrm{d}\vec{r}'\,\vec{e}_{k}(\vec{r}')\leq \|\vec{e}_k\|_{L^1(U)}\leq |U|^{1/2}\|\vec{e}_k\|_{L^2(U)} = |U|^{1/2}$ by orthonormality of the basis. Now for $n=2$ we have $$\begin{aligned}
\text{Tr}\mathcal{A}^2_\varrho &= \sum_{k=1}^d\langle \mathcal{A}^2_\varrho \vec{e}_k(\vec{r}), \vec{e}_k(\vec{r})\rangle_{L^2(U,\varrho^{-1})}\\
&=\sum_{k=1}^d\int_U\mathrm{d}\vec{r}\,\vec{e}_{k}(\vec{r})\cdot\varrho(\vec{r})^{-1}\int_U\mathrm{d}\vec{r}_1\,\varrho(\vec{r})\bm{D}(\vec{r})g(\vec{r},\vec{r}_1)\bm{Z}_2(\vec{r},\vec{r}_1)\\
&\quad \times\int_U\mathrm{d}\vec{r}_2\,\varrho(\vec{r}_1)\bm{D}(\vec{r}_1)g(\vec{r}_1,\vec{r}_2)\bm{Z}_2(\vec{r}_1,\vec{r}_2)\vec{e}_{k}(\vec{r}_2)\\
&\leq \sum_{k=1}^d\mu_{\max}^2\|g\bm{Z}_2\|_{L^\infty(U)}^2\int_U\mathrm{d}\vec{r}\,\vec{e}_{k}(\vec{r})\cdot\int_U\mathrm{d}\vec{r}_2\,\vec{e}_{k}(\vec{r}_2)\\
&\leq d\mu_{\max}^2\|g\bm{Z}_2\|_{L^\infty(U)}^2|U|,\end{aligned}$$ where we have used the fact that $\int_U\mathrm{d}\vec{r}_1\,\varrho(\vec{r}_1,t) = 1$ for $t\geq 0$ by the no-flux boundary condition (see Section \[subsec:boundary\_conditions\]). Iterating this argument one may obtain $$\begin{aligned}
\text{Tr}\mathcal{A}^n_\varrho &= \sum_{k=1}^d\langle \mathcal{A}^n_\varrho \vec{e}_k(\vec{r}), \vec{e}_k(\vec{r})\rangle_{L^2(U,\varrho^{-1})}\\
&\leq \mu_{\max}^n\|g\bm{Z}_2\|_{L^\infty(U)}^n\sum_{k=1}^d\int_U\mathrm{d}\vec{r}\,\vec{e}_{k}(\vec{r})\cdot\int_U\mathrm{d}\vec{r}_n\,\vec{e}_{k}(\vec{r}_n)\\
& = |U|d\times \mu_{\max}^n\|g\bm{Z}_2\|_{L^{\infty}(U)}^n.\end{aligned}$$

We observe that absolute convergence of requires $\mu_{\max}\|g\bm{Z}_2\|<1$. In particular, the sum of the absolute values of the terms is given by $$\begin{aligned}
\sum_{n=1}^\infty\Big|\tfrac{\text{Tr}(\mathcal{A}_\varrho^n)}{n}\Big| \leq d|U|\sum_{n=1}^\infty\tfrac{\mu_{\max}^n\|g\bm{Z}_2\|_{L^\infty(U)}^n}{n} = d|U|\log\left(\frac{1}{1-\mu_{\max}\|g\bm{Z}_2\|_{L^\infty(U)}}\right).\end{aligned}$$ Thus for $\mu_{\max}||g \bm{Z}_2||_{L^\infty(U)}<1$ the logarithm is finite and the determinant is positive, otherwise for the boundary case $\mu_{\max}||g \bm{Z}_2||_{L^\infty(U)} =1$ it may vanish, thus making $(\bm{I} + \mathcal{Z}_1^\varrho + \mathcal{Z}_2^\varrho)$ singular.

We now provide a scheme for computing solutions of equation for each time dependent $\varrho(\vec{r},t)$. The existence and uniqueness of $\varrho(\vec{r},t)$ is given in Section \[sec:existence\_uniqueness\_with\_partial\_HI\]. First we establish that $\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2$ is a compact self-adjoint operator in $L^2(U,\varrho^{-1})$.\

\[lem:flux\_operator\_is\_compact\_self\_adjoint\] Let $\lambda\in (-\infty,\infty)$ and assumption hold. Then $\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2$ is a compact and self-adjoint operator.

We let $\vec{a}\in L^1(U)$ and calculate $\|(\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2)[\vec{a}]\|_{L^1(U)}$. In particular we have $$\begin{aligned}
\|(\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2)[\vec{a}]\|_{L^1(U)} &= \int\mathrm{d}\vec{r}\,\Big|(\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2)[\vec{a}]\Big|\nonumber\\
&\leq \int\mathrm{d}\vec{r}\,\Big|(\bm{1}+\mathcal{Z}^\varrho_1)\vec{a}\Big|+|\lambda|\Big|\mathcal{Z}^\varrho_2[\vec{a}]\Big|\nonumber\\
&\leq (1+\|g \bm{Z}_1\|_{L^\infty(U)}+|\lambda|\|g\bm{Z}_2\|_{L^\infty(U)})\|\vec{a}\|_{L^1(U)}<\infty.\end{aligned}$$ Hence $\text{Im}(\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2)$ is bounded in $\mathbb{R}^3$. Now by Heine–Borel, the closure of $\text{Im}(\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2)$ is compact and hence $\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2$ is a compact operator.

We now show that $\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2$ is self-adjoint. The local part $\bm{1}+\mathcal{Z}^\varrho_1$ is a real, symmetric matrix and it is therefore self-adjoint, and in particular self-adjoint in $L^2(U,\varrho^{-1})$ . All that remains is to study the nonlocal part $\mathcal{Z}^\varrho_2$. By direct calculation we see that for $\vec{b}\in L^2(U)$ $$\begin{aligned}
\langle \vec{b},\mathcal{Z}_2^\varrho[\vec{a}]\rangle_{L^2(U,\varrho^{-1})} &= \int_U\mathrm{d}\vec{r}\,\vec{b}(\vec{r})\cdot \int_U\mathrm{d}\vec{r}'\,g(\vec{r},\vec{r}')\bm{Z}_2(\vec{r},\vec{r}')\vec{a}(\vec{r}')\\
&=\int_U\mathrm{d}\vec{r'}\,\vec{b}(\vec{r})^\top \int_U\mathrm{d}\vec{r}\,g(\vec{r},\vec{r}')\bm{Z}_2(\vec{r},\vec{r}')^\top\vec{a}(\vec{r}')\\
&=\int_U\mathrm{d}\vec{r'}\,\int_U\mathrm{d}\vec{r}\,\left(g(\vec{r},\vec{r}')\bm{Z}_2(\vec{r},\vec{r}') \vec{b}(\vec{r})\right)^\top \vec{a}(\vec{r}')\\
&=\int_U\mathrm{d}\vec{r'} \vec{a}(\vec{r}')\cdot\int_U\mathrm{d}\vec{r}\,g(\vec{r},\vec{r}')\bm{Z}_2(\vec{r},\vec{r}') \vec{b}(\vec{r})\\
&=\langle \mathcal{Z}_2^\varrho[\vec{b}],\vec{a}\rangle_{L^2(U,\varrho^{-1})}\end{aligned}$$ where we have used the symmetry of $\bm{Z}_2$, and on the last line used Fubini’s theorem to interchange the order of the integration between the $\vec{r}'$ and $\vec{r}$ variables. Hence the lemma is proved.

Since we have now established that $\bm{1}+\mathcal{Z}^\varrho_1-\lambda\mathcal{Z}^\varrho_2$ is a compact and self-adjoint operator we may use its eigenvectors as a complete basis of $\mathbb{R}^3$ to expand the flux $\vec{a}(\vec{r},t)$.\

\[thm:eigenfn\_expansion\_of\_flux\]

Let $\bm{Z}_2$ be symmetric and real and $\mu_{\max}\|g\bm{Z}_2\|_{L^\infty(U)}<1$ and let $\vec{e}_i(\vec{r},[\varrho],t)$ and $\mu_i^{-1}$ be the eigenvectors and eigenvalues of $\bm{D}^{-1}(\vec{r},[\varrho],t)$ where $[\cdot]$ denotes functional dependence and $i = 1,\cdots, d$. Then there is a unique $\vec{a}(\vec{r},t)\in L^2(U)$ solving given by the eigenfunction expansion $$\begin{aligned}
\vec{a}(\vec{r},t) = \sum_{n=1}^d\delta_n\vec{w}_{n}(\vec{r},t)\label{eq:a_expanded_in_wn}.\end{aligned}$$ Here, $\vec{w}_n$ are eigenfunctions of $(\bm{1}+\mathcal{Z}_1^\varrho -\lambda\mathcal{Z}_2^\varrho)$ obtained by a second expansion in $\vec{e}_i(\vec{r},[\varrho],t)$ of the form $$\begin{aligned}
\label{eq:wn_expanded_in_es}
\vec{w}_n(\vec{r},t) = \sum_{i=1}^d\frac{\psi_i}{\mu_i^{-1}-\phi_n} \vec{e}_i(\vec{r},[\varrho],t).\end{aligned}$$

Additionally, the expansion coefficients $\delta_n$ are given by the formula $$\begin{aligned}
\label{eq:eqn_for_deltans}
\delta_n = \frac{1}{\phi_n}\sum_{i=1}^d \frac{\psi_i}{\phi_n-\mu_i^{-1}}\int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t)\vec{e}_i(\vec{r},[\varrho],t)\cdot\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\end{aligned}$$ where $\{\phi_n\}_{n=1}^d$ are the discrete set of eigenvalues of $(\bm{1}+\mathcal{Z}_1^\varrho -\lambda\mathcal{Z}_2^\varrho)$ given by roots of the equation $\lambda(\phi_n) = -1$, where the function $\lambda(\cdot)$ is defined by $$\begin{aligned}
\lambda(\phi_n) := \left[\sum_{l=1}^d\frac{\eta_l \,\psi_l}{\mu_l^{-1}-\phi_n}\right]^{-1}.\end{aligned}$$

Finally, $\psi_k$ and $\eta_l$ are the expansion coefficients defined by $$\begin{aligned}
\varrho(\vec{r},t)g(\vec{r},\vec{r}')\bm{Z}_2(\vec{r},\vec{r}') = \sum_{k=1}^d\sum_{l=1}^d
\psi_k\eta_l \vec{e}_k(\vec{r},[\varrho],t)\otimes\vec{e}_l(\vec{r}',[\varrho],t)\label{eq:rho_g_Z2_expanded_in_e}\end{aligned}$$ and each $\varrho$ is obtained from the continuity equation and no-flux condition $$\begin{aligned}
\partial_t \varrho &= -\nabla_{\vec{r}}\cdot \vec{a},\\
0 & =\Pi[\varrho]\cdot \vec{n}|_{\partial U}.\end{aligned}$$

The scalars $\mu_i^{-1}$, $\psi_k$, $\eta_l$ (and by proxy $\delta_n$) each have functional dependence on $\varrho$ since they are obtained by integrals involving $\vec{e}_j(\vec{r},[\varrho],t)$, for $i,j,k,l = 1,\cdots, d$. The eigenvalues $\phi_n$ are so called ‘moving eigenvalues’ of $(\bm{1}+\mathcal{Z}_1^\varrho -\lambda\mathcal{Z}_2^\varrho)$ (cf. [@davidson2006spectral]). If $\bm{Z}_2 = 0$ then $\phi_i = \mu_i^{-1}$ for each $i=1,\cdots, d$. In general, for $\bm{Z}_2 \neq 0$, an eigenvalue of $\bm{D}^{-1}$ may also be an eigenvalue of $(\bm{1}+\mathcal{Z}_1^\varrho -\lambda\mathcal{Z}_2^\varrho)$ and this occurs on the line $\lambda = 0$. Since $\bm{Z}_2$ is symmetric it can be diagonalised, and therefore the kernel of the operator $\mathcal{Z}_2$ can be decomposed into a finite (of length $d$) sum of products of continuous functions and has at most $d$ eigenvalues. The equation $\lambda(\phi_n) = -1$ may be rearranged into a characteristic polynomial equation in $\phi_n$ with coefficients dependent on $\eta_l$, $\psi_l$ and $\mu_l$ and since $(\bm{1}+\mathcal{Z}_1^\varrho -\lambda\mathcal{Z}_2^\varrho)$ is assumed to be real and symmetric, each $\phi_n\in\mathbb{R}$. Finally, the condition $\mu_{max}\|g\bm{Z}_2\|_{L^\infty(U)}<1$ ensures $\phi_n \neq 0$ for any $n\in \mathbb{N}$.

We consider the more general operator $(\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho)$ where $\lambda \in \mathbb{R}$. One may think of this operator as a nonlocal matrix operator where $(\bm{1}+\mathcal{Z}_1^\varrho)$ is the local part and $\mathcal{Z}_2^\varrho$ is the nonlocal part. Here $\lambda$ is a perturbation parameter measuring the distance of the full operator $(\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho)$ from locality. Since $\bm{Z}_1$ and $\bm{Z}_2$ are real and symmetric and $\lambda\in \mathbb{R}$ , $\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho$ coincides with its adjoint in $L^2(U,\varrho^{-1})$. For the homogeneous adjoint equation $$\begin{aligned}
\label{eq:adjoint_eqn}
(\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho)^\dagger\bm{z} = 0\end{aligned}$$ we know from Lemma \[thm:cond\_converg\_fred\_det\] that when $\mu_{\max}\|g\bm{Z}_2\|_{L^\infty(U)}<1$ there is no $\lambda\in \mathbb{R}$ satisfying $\det((\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho)^\dagger) = 0$ and therefore the only solution to the homogeneous adjoint equation is $\bm{z} = \vec{0}$. Therefore by the Fredholm alternative there is a unique solution to .

Now consider the eigenvalue problem $$\begin{aligned}
\label{eq:eigenvalue_problem_(1+lambdaA)}
(1+\mathcal{Z}_1^\varrho-\lambda \mathcal{Z}_2^\varrho)[\vec{w}_n(\vec{r},t)] = \phi_n\vec{w}_n(\vec{r},t)\end{aligned}$$ for eigenvalues $\phi_n\in\mathbb{R}$ and eigenvectors $\vec{w}_n\in \mathbb{R}^d$. We write $$\begin{aligned}
\label{eq:wn_expanded_in_ejs}
\vec{w}_n = \sum_{j=1}^d\alpha_{j,n}\vec{e}_j(\vec{r},[\varrho],t).\end{aligned}$$ By inserting into we obtain $$\begin{aligned}
\label{eq:eigenvalue_problem_expanded}
(\bm{1}+\mathcal{Z}_1^\varrho)\sum_{j=1}^d\alpha_{j,n}\vec{e}_j(\vec{r},[\varrho],t) - \lambda \mathcal{Z}_2^\varrho\left[ \sum_{j=1}^d\alpha_{j,n}\vec{e}_j(\vec{r},[\varrho],t)\right]  = \phi_n\sum_{j=1}^d\alpha_{j,n}\vec{e}_j(\vec{r},[\varrho],t).\end{aligned}$$ Now by inserting the expansion into we obtain $$\begin{aligned}
\sum_{j=1}^d\alpha_{j,n}(\mu_j^{-1}-\phi_n)\vec{e}_{j}(\vec{r},[\varrho],t) - \lambda\sum_{k,l=1}^d\psi_k\,\eta_l\int_U\mathrm{d}\vec{r}'\vec{e}_k(\vec{r},[\varrho],t)\otimes \vec{e}_l(\vec{r}',[\varrho],t)\vec{w}_n(\vec{r}',t) = 0.\end{aligned}$$ Taking the inner product of this equation with $\vec{e}_i(\vec{r},[\varrho],t)$ and integrating we obtain $$\begin{aligned}
\alpha_{i,n}(\mu^{-1}_i-\phi_n) -\lambda\,\psi_i\sum_{l=1}^d\eta_l\int_U\mathrm{d}\vec{r}'\vec{e}_l(\vec{r},[\varrho],t)\cdot \vec{w}_n(\vec{r}',t)=0,\end{aligned}$$ which may be rearranged to obtain $$\begin{aligned}
\label{eq:separtion_eqn_for_lambda}
\lambda  = \frac{\alpha_{i,n}(\mu_i^{-1}-\phi_n)}{\psi_i\,\sum_{l=1}^d\eta_l\int_U\mathrm{d}\vec{r}'\vec{e}_l(\vec{r},[\varrho],t)\cdot \vec{w}_n(\vec{r}',t)}.\end{aligned}$$ Since both the left hand side of and $\sum_{l=1}^d\eta_l\int_U\mathrm{d}\vec{r}'\vec{e}_l(\vec{r},[\varrho],t)\cdot \vec{w}_n(\vec{r}',t)$ are independent of the index $i$ it must be that $$\begin{aligned}
\frac{\alpha_{i,n}(\mu_i^{-1}-\phi_n)}{\psi_i} = K\end{aligned}$$ for some constant $K$ for which, without loss of generality, we choose $K=1$. With this we obtain an expression for the coefficients $\alpha_{i,n}$ $$\begin{aligned}
\label{eq:eqn_for_alpha_in}
\alpha_{i,n} = \frac{\psi_i}{\mu_i^{-1}-\phi_n}.\end{aligned}$$

We may also obtain a scheme to determine the $\phi_n$. In particular by and we have $$\begin{aligned}
\lambda &= \left(\sum_{l=1}^d\eta_l\int_U\mathrm{d}\vec{r}'\vec{e}_l(\vec{r},[\varrho],t)\cdot \vec{w}_n(\vec{r}',t)\right)^{-1}\nonumber\\
& = \left(\sum_{l=1}^d\eta_l\int_U\mathrm{d}\vec{r}'\vec{e}_l(\vec{r},[\varrho],t)\cdot \sum_{j=1}^d\frac{\psi_j}{\mu_j^{-1}-\phi_n}\vec{e}_j(\vec{r}',[\varrho],t)\right)^{-1} = \left(\sum_{l=1}^d\frac{\eta_l\,\psi_l}{\mu_l^{-1}-\phi_n}\right)^{-1}\end{aligned}$$ hence we have that the eigenvalues of $(\bm{1}+\mathcal{Z}_1^\varrho+\mathcal{Z}_2^\varrho)$ are given by the roots of the equation $\lambda(\phi_n) = -1$.

We now return to the inhomogeneous problem and expand $\vec{a}(\vec{r},t)$ in eigenfunctions $\vec{w}_n(\vec{r},t)$. We propose an expansion of the form and insert into to obtain $$\begin{aligned}
\sum_{n=1}^d \delta_n\phi_n\vec{w}_n(\vec{r},t) = -\varrho(\vec{r},t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho].\end{aligned}$$ Now by taking the inner product with some $\vec{w}_k(\vec{r},t)$ and integrating we obtain $$\begin{aligned}
\delta_k\phi_k = -\int_U\mathrm{d}\vec{r}\,\varrho(\vec{r},t)\vec{w}_k(\vec{r},t)\cdot \nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho].\end{aligned}$$ By inserting the definition of $\vec{w}_k$ from we deduce $$\begin{aligned}
\label{eq:eqn_for_delta_k_coeffs}
\delta_k\phi_k = \sum_{i=1}^d \frac{\psi_i}{\phi_k-\mu_i^{-1}}\int_U\mathrm{d}\vec{r}\, \vec{e}_i(\vec{r},[\varrho],t)\cdot\varrho(\vec{r},t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho].\end{aligned}$$

Now we would like to divide through by $\phi_k$ but must check that no $\phi_k$ is zero for each $k = 1,\cdots, d$. This is a consequence of the condition $\mu_{\max}\|g \bm{Z}_2\|<1$. In particular, using properties of the determinant, we have that $$\begin{aligned}
\det(\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho)  = \det(\bm{1}+\mathcal{Z}_1^\varrho) \times\det(\bm{1}-\lambda(\bm{1}+\mathcal{Z}_1^\varrho)^{-1}\mathcal{Z}_2^\varrho).\end{aligned}$$ Now since $\bm{D}$ is positive definite, so is $\bm{1}+\mathcal{Z}_1^\varrho$ and therefore $\det(\bm{1}+\mathcal{Z}_1^\varrho) >0$ because the determinant is simply the product of its (strictly positive) eigenvalues. Additionally, since $\mu_{\max}\|g \bm{Z}_2\|<1$, we have by Lemma \[thm:cond\_converg\_fred\_det\] $\det(\bm{1}-\lambda(\bm{1}+\mathcal{Z}_1^\varrho)^{-1}\mathcal{Z}_2^\varrho)>0$ therefore $\det(\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho)>0$ and $\phi_k \neq 0$ for all $k\in \mathbb{N}$. We may now divide by $\phi_k$ to obtain . Finally $\vec{a}(\vec{r},t)\in L^2(U)$ may be seen by squaring , integrating over $\mathrm{d}\vec{r}$ and using .

Theorem \[thm:eigenfn\_expansion\_of\_flux\] provides a scheme for computing the unique flux $\vec{a}(\vec{r},t)$, given $\varrho$ satisfying $\partial_t\varrho = -\nabla_{\vec{r}}\cdot \vec{a}$ over time. We now use this result to show that the free energy functional $\mathcal{F}[\varrho]$ may be associated to the full system even when $\bm{Z}_2 \neq 0$. In particular, that $\varrho(\vec{r},t)$ solving implies $\varrho$ is a critical point of the free energy $\mathcal{F}[\varrho]$.\

\[thm:association \_of\_free\_energy\]

Let $\mu_{\max}\|g\bm{Z}_2\|_{L^\infty(U)}<1$ $V_1 = V_1(\vec{r})$ be a time independent confining potential so that $\varrho(\vec{r})$ is a stationary density to the system then $\varrho(\vec{r})$ is a critical point of $\mathcal{F}[\varrho]$.

Let $\varrho(\vec{r},t) = \varrho(\vec{r})$ be a stationary density. Then by equation one has $$\begin{aligned}
(\bm{1}+(\bm{1}+\mathcal{Z}_1^\varrho)^{-1}\mathcal{Z}_2^\varrho)[\vec{a}(\vec{r})]&=-\varrho(\vec{r})\bm{D}(\vec{r},[\varrho],t)\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho],\label{eq:stationary_eqn_for_a_dimensionless}\\
\nabla_{\vec{r}}\cdot \vec{a}(\vec{r})&=0.\label{eq:incomp_a_at_equilibrium}\end{aligned}$$ We have that for each $\lambda$, the operator $\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho$ is compact self-adjoint in $L^2(U,\varrho^{-1}(\vec{r},t))$ (by Lemma \[lem:flux\_operator\_is\_compact\_self\_adjoint\]). We also have that $\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho$ is positive definite for $\mu_{\max}\|g Z_2\|_{L^\infty(U)}<1$. In particular, $\phi_n\neq 0$ for every $n = 1,\cdots d $ and $\phi_n(\lambda)$ is continuous function of $\lambda$ such that $\phi_n(0) = \mu_{n}^{-1}>0$ for each $n$. Hence we may invert $\bm{1}+\mathcal{Z}_1^\varrho+\mathcal{Z}_2^\varrho$ given $\mu_{\max}\|g Z_2\|_{L^\infty(U)}<1$. With this, by using equations , we have $$\begin{aligned}
\label{eq:div_stationary_eqn_in_terms_of_F}
0 = \nabla_{\vec{r}}\cdot \vec{a} = \nabla_{\vec{r}}\cdot \left(\varrho(\vec{r})(\bm{1}+\mathcal{Z}_1^\varrho+\mathcal{Z}_2^\varrho)^{-1}\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\right). \end{aligned}$$

Now, assuming $\varrho$ is stationary we see that $$\begin{aligned}
0 = \Big\langle \frac{\delta\mathcal{F}}{\delta\varrho}[\varrho],\partial_t\varrho\Big\rangle = -\int_U\mathrm{d}\vec{r}\, \frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\nabla_{\vec{r}}\cdot \vec{a} = \int_U\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\cdot \vec{a}\end{aligned}$$ where we have used the no-flux boundary condition. Now since $(\bm{1}+\mathcal{Z}_1^\varrho+\mathcal{Z}_2^\varrho)^{-1}$ is strictly positive definite and self-adjoint in $L^2(U,\varrho^{-1}(\vec{r},t))$ it possesses a unique strictly positive definite self-adjoint square root in $L^2(U,\varrho^{-1}(\vec{r},t))$ (see [@wouk1966note]). We define $\mathcal{X}_\varrho = (\bm{1}+\mathcal{Z}_1^\varrho+\mathcal{Z}_2^\varrho)^{-1}$ and $\mathcal{X}_\varrho^{1/2}\mathcal{X}_\varrho^{1/2} = \mathcal{X}_\varrho$. Then we find $$\begin{aligned}
0 &= \int_U\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\cdot \vec{a} = \int_U\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\cdot \mathcal{X}_\varrho \left[\varrho\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\right]\nonumber\\ 
&= \Big\langle \varrho \nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho],\mathcal{X}_\varrho \left[\varrho\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\right]\Big\rangle_{L^2(U,\varrho^{-1})}\nonumber\\
&=\Big\langle \mathcal{X}_\varrho^{1/2} \left[\varrho\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\right],\mathcal{X}_\varrho^{1/2} \left[\varrho\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\right]\Big\rangle_{L^2(U,\varrho^{-1})}\nonumber\\
&=\Big\|\mathcal{X}_\varrho^{1/2} \left[\varrho\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\right]\Big\|_{L^2(U,\varrho^{-1})}^2\label{eq:X_sqrt_energy_argument}\end{aligned}$$ where we have used the self-adjoint property of $\mathcal{X}_\varrho^{1/2} $. From the above we deduce that, since the integrand in the last line of is positive, that the stationary density $\varrho(\vec{r})$ satisfies $$\begin{aligned}
\label{eq:grad_F_rho_is_zero}
\varrho(\vec{r})\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho(\vec{r})] = \vec{0}.\end{aligned}$$ Therefore we obtain that $\varrho$ is a critical point the free energy $\mathcal{F}[\varrho]$.

Let $\mu_{\max}\|g\bm{Z}_2\|_{L^\infty(U)}<1$ $V_1 = V_1(\vec{r})$ be a time independent confining potential. If $\varrho$ is a stationary density then it is a critical point of the free energy $\mathcal{F}[\varrho]-\int\mathrm{d}\vec{r}\mu_c\varrho$.

From Theorem \[thm:association \_of\_free\_energy\] we obtain the following two corollaries. In particular, the proof shows rigorously how the diffusion tenor decouples from the stationary density.\

\[cor:stationary\_density\_independent\_of\_D\]

Let $\bm{Z}_1$, $\bm{Z}_2$ be real and symmetric. Then the stationary density $\varrho$ is independent of $\bm{Z}_1$, $\bm{Z}_2$ and, as a consequence, of $\bm{D}$.

Additionally since holds in equilibrium even when $\bm{Z}_2 \neq 0$ and the condition $\mu_{\max}\|g\bm{Z}_2\|_{L^
\infty(U)}<1$ implies that the operator $(\bm{1}+(\bm{1}+\mathcal{Z}_1^\varrho)^{-1}\mathcal{Z}_2^\varrho)$ has no zero eigenvalue, the homogeneous problem $(\bm{1}+(\bm{1}+\mathcal{Z}_1^\varrho)^{-1}\mathcal{Z}_2^\varrho) = \vec{0}$ (i.e. at equilibrium) must have only the trivial solution $\vec{a}(\vec{r}) = \vec{0}$. In addition by equation , at equilibrium one has $$\begin{aligned}
\label{eq:stationary_eqn_rigorously_derived}
\nabla_{\vec{r}}\cdot \left(\varrho(\vec{r})\bm{D}(\vec{r},[\varrho])\nabla_{\vec{r}}\frac{\delta\mathcal{F}}{\delta\varrho}[\varrho]\right)=0\end{aligned}$$ where $\bm{D}(\vec{r},[\varrho])$ is the time limiting diffusion tensor.\

\[cor:a\_is\_zero\_at\_equilibrium\]

Let $\bm{Z}_1$, $\bm{Z}_2$ be real and symmetric. Then $\vec{a}(\vec{r}) = \vec{0}$ is the unique stationary flux. In particular, there do not exist stationary densities which are advected by the existence of some non-zero flux, hence the only stationary states are equilibrium states.

We remark that Corollary \[cor:stationary\_density\_independent\_of\_D\] is related to the well-known result that for finite dimensional reversible diffusions, i.e. Langevin dynamics of the form $dX_t = - (D((X_t)) \nabla V((X_t))) \, dt + \nabla \cdot D(X_t) \, dt + \sqrt{2 D(X_t)} \, dW_t$ for an arbitrary strictly positive definite mobility matrix $D$, $V$ a confining potential and Wiener process $W_t$, the invariant measure $\mu(dx) = \frac{1}{Z} e^{-V(x)} \, dx$ is independent of $D$. We refer to [@pavliotis2014stochastic Sec 4.6]. To our knowledge, this is the first instance where such a result is proved in the context of DDFT.

In the following Sections \[sec:char\_stationary\_sol\], \[sec:global\_asymptotic\_stability\], \[sec:bifurcation\_theory\], we consider the global asymptotic stability of the stationary equation (equivalently ) for which, we have shown by Corollary \[cor:stationary\_density\_independent\_of\_D\], that is the equation determining the equilibrium density the dynamics driven to equilibrium by the HI tensors $\bm{Z}_1,\bm{Z}_2$.\

\[rem:final\_assumptions\]

Out of equilibrium, the effective drift is augmented by $\bm{A}[\vec{a}]$ (as defined in ), the flow induced by the HI. In order to simplify the presentation of the calculations needed for the proofs of several results presented later on, (Theorem \[lem:exp\_conv\_in\_L2\], Proposition \[prop:general\_kappa1\_convergence\], all results in Sections \[sec:existence\_uniqueness\_with\_partial\_HI\] and \[sec:classical\_paraboliv\_pde\]) we suppress $\bm{A}[\vec{a}]$ because it may trivially included as a linear contribution which is bounded in $L^1(U)$: $$\begin{aligned}
\|\bm{A}[\vec{a}]\|_{L^1(U)} = \int_U\mathrm{d}\vec{r}\,\Big| \int_U\mathrm{d}\vec{r}'\, \bm{Z}_2(\vec{r},\vec{r}')\vec{a}(\vec{r}',t)\Big|\leq \|\bm{Z}_2\|_{L^\infty(U)}\|\vec{a}\|_{L^1(U)}<\infty\end{aligned}$$ where we have used and the fact that, by Theorem \[thm:eigenfn\_expansion\_of\_flux\], $\|\vec{a}\|_{L^1(U)}^2\leq |U|\|\vec{a}\|^2_{L^2(U)}<\infty$. Hence all the coefficients of remain uniformly bounded and the existence and uniqueness results of Section \[sec:existence\_uniqueness\_with\_partial\_HI\] may be easily obtained with $\bm{A}[\vec{a}]$ included.

Additionally, since we have shown that at equilibrium $\bm{A}[\vec{a}]= \vec{0}$ uniquely, the results of Sections \[sec:global\_asymptotic\_stability\], \[sec:bifurcation\_theory\] hold for the dynamics tending to equilibrium including the effects of the HI.

Given this remark, we now discuss the existence of stationary solutions to .\

Characterisation of Stationary Solutions {#sec:char_stationary_sol}
========================================

We now define the stationary problem. \[sec:stationary\_problem\] We seek classical solutions $\varrho\in C^2(\bar{U})$ of

$$\begin{aligned}
\nabla_{\vec{r}}\cdot \left[\bm{D}\left( \nabla_{\vec{r}}\varrho+\varrho\nabla_{\vec{r}}[\kappa_1 V_1+\kappa_2V_2\star\varrho]\right) \right]   &= 0 \qquad \vec{r}\in U, \label{eq:stationary_problem}\\
\Pi[\varrho]\cdot \vec{n} &= 0 \qquad \vec{r} \text{ on } \partial_U. \label{eq:stationary_problem_bc}\end{aligned}$$

where $$\begin{aligned}
\Pi[\varrho]:=\bm{D}\left( \nabla_{\vec{r}}\varrho+\varrho\nabla_{\vec{r}}[\kappa_1 V_1+\kappa_2V_2\star\varrho]\right).\end{aligned}$$ The existence and uniqueness for the stationary problem is based on a fixed point argument for the nonlinear map, defined by integrating equation . In particular we find the stationary distribution satisfies the nonlinear map (the self-consistency equation) $$\begin{aligned}
\label{eq:self_cons_eq}
\varrho(\vec{r}) = \frac{e^{-(\kappa_1V_1(\vec{r})+ \kappa_2V_2\star\varrho(\vec{r}))}}{Z},\end{aligned}$$ where $Z = \int_U\mathrm{d}\vec{r}\, \exp\{-(\kappa_1V_1(\vec{r})+ \kappa_2V_2\star\varrho(\vec{r}))\}$. Note that the stationary distribution is independent of the diffusion tensor (see Corollary \[cor:stationary\_density\_independent\_of\_D\]). We now present our first result concerning the existence and uniqueness of the solutions to the self-consistency equation.\

\[thm:exis\_fix\_point\]

The stationary equation with boundary condition has a smooth, non-negative solution with $\|\varrho\|_{L^1(U)}=1$. When the interaction energy is sufficiently small, $| \kappa_2 | \leq 1/4\times \|V_2\|_{L^\infty}^{-1}$, the solution is unique.

The proof follows Dressler et al. [@dressler1987stationary]. The main idea is to show that the right hand side of equation is a contraction map on $C^2(U)$, and for sufficiently small interaction energy $\kappa_2$, $\varrho_\infty\in L^1(U)$ is the unique invariant measure which is a non-negative function with unit mean.

Consider the stationary problem such that Assumption holds. Then we have that

1.  There exists a weak solution $\varrho\in H^1(U)\cap P_{ac}(U)$ to as a fixed point of the equation .

2.  Any weak solution $\varrho\in H^1(U)\cap P_{ac}(U)$ is smooth and strictly positive, that is $\varrho\in C^{\infty}(\bar{U})\cap P^+_{ac}(U)$.

The proof is similar to [@greg_mckean_vlasov_torus Theorem 2.3] but one must check the conclusions of the theorem hold with no flux boundary conditions and a confining potential $V_1$. This result is similar to arguments in [@tamura1984asymptotic] but here we consider a compact domain $U$. The weak formulation of is $$\begin{aligned}
\label{eq:weak_form_of_stationary_eqn}
-\int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}\eta\cdot \bm{D}\nabla_{\vec{r}}\varrho -\kappa_1 \int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}\eta\cdot \varrho\bm{D}\nabla_{\vec{r}}V_1 - \kappa_2 \int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}\eta\cdot \varrho\bm{D}\nabla_{\vec{r}}V_2\star \varrho = 0,\end{aligned}$$ for $\eta \in H^1(U)$ where we have used the no-flux boundary condition in on $\varrho$ and we seek solutions $\varrho\in H^1(U)\cap P_{ac}(U)$. Now define $F:P_{ac}(U)\to P_{ac}(U)$ by $$\begin{aligned}
\label{eq:def_of_T_rho}
F \varrho = \frac{1}{Z(\varrho,\kappa_2)}e^{-(\kappa_1V_1+ \kappa_2V_2\star \varrho)}, \quad Z(\varrho,\kappa_2) = \int_U\mathrm{d}\vec{r}\,e^{-(\kappa_1V_1+ \kappa_2V_2\star \varrho)}.\end{aligned}$$ By we see that $$\begin{aligned}
\label{eq:T_on_E}
\|F\varrho\|_{L^2(U)}^2 \leq \frac{1}{|U|} e^{4(|\kappa_1|\|V_1\|_{L^\infty(U)}+|\kappa_2|\|V_2\|_{L^\infty(U)})} =: E_0,\end{aligned}$$ and therefore we seek solutions to in the set $E := \{\varrho\in L^2(U) \,:\, \|\varrho\|_{L^2(U)}^2 \leq E_0\}$. Note that $E$ is a closed, convex subset of $L^2(U)$ and therefore we may redefine $T$ to act on $E$. Additionally we see that for $\varrho\in E$ $$\begin{aligned}
\|F\varrho\|_{H^1(U)}^2 &= \|F\varrho\|_{L^2(U)}^2 + \|\nabla_{\vec{r}}T\varrho\|_{L^2(U)}^2\nonumber\\
&\leq E_0\left(1+2|\kappa_1|^2\|\nabla V_1\|_{L^2(U)}^2+|\kappa_2|^2\|\nabla V_2\|_{L^2(U)}^2E_0\right),\label{eq:H1_estimate_of_rho}\end{aligned}$$ where we have used that $\varrho\in L^1(U)$ by Lemma \[thm:exis\_fix\_point\] and $V_1,V_2\in H^1(U)$. Similarly to [@greg_mckean_vlasov_torus Theorem 2.3] we have by that $F(E)\subset E$ and by $F(E)$ is uniformly bounded in $H^1(U)$. Therefore by Rellich’s compactness theorem, $F(E)$ is relatively compact in $L^2(U)$, and therefore in $E$, since $E$ is closed.

We may show using similar calculations to [@dressler1987stationary Theorem 1] that the non-linear map in is Lipschitz continuous in $E$, and by Schauder fixed point theorem there exists $\varrho\in E$ solving which by is in $H^1(U)$. By inserting the expression for $F\varrho$ into we obtain (1). Also note that solutions $\varrho\in E$ to are bounded below by $E_0^{-1}/|U|^2$ giving positivity of solutions.

We now show that every weak solution in $\varrho\in H^1(U)\cap P_{ac}(U)$ is a fixed point of the nonlinear map in . Consider the frozen weak formulation $$\begin{aligned}
\label{eq:weak_form_of_stationary_eqn_frozen}
-\int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}\eta\cdot \bm{D}\nabla_{\vec{r}}\theta -\kappa_1 \int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}\eta\cdot \bm{D}\nabla_{\vec{r}}V_1\theta - \kappa_2 \int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}\eta\cdot \bm{D}\nabla_{\vec{r}}V_2\star \varrho\, \theta = 0.\end{aligned}$$ This is the weak formulation of the PDE (for the unknown function $\theta$) $$\begin{aligned}
\nabla_{\vec{r}}\cdot\left(\bm{D}\nabla_{\vec{r}}\theta + \theta\bm{D}(\nabla_{\vec{r}}V_1+\nabla_{\vec{r}}V_2\star\varrho)\right) = 0, \quad \text{ s.t. }  \nabla_{\vec{r}}\left((F\varrho)^{-1}\theta\right)\cdot\vec{n}|_{\partial U} = 0.\end{aligned}$$ We note that we may rewrite the weak formulation as $$\begin{aligned}
\label{eq:weak_form_of_stationary_eqn_frozen_rewrite}
-\int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}\eta\cdot \bm{D}\nabla_{\vec{r}}h\,F\varrho = 0\end{aligned}$$ for every $\eta\in H^1(U)$ and where $h = \theta/(F\varrho)$. This holds true for any $\eta$, in particular $\eta = h$ hence we find $$\begin{aligned}
-\int_U\mathrm{d}\vec{r}\,\Big| (F\varrho)^{1/2}\bm{D}^{1/2}\nabla_{\vec{r}}h\Big|^2 = 0\end{aligned}$$ where we have used that $\bm{D}$ is positive definite by and $F\varrho$ is strictly positive. All in all we obtain $\nabla_{\vec{r}}h = 0$ a.e. and hence $\theta = F\varrho$ up to normalisation. But if $F\varrho$ is a probability density we must have $\theta \equiv F\varrho$ and we conclude that since $\varrho = F\varrho$, any weak solution $\varrho\in H^1(U)\cap P^+_{ac}(U)$ of must be such that $\varrho = F\varrho$. The regularity of $\varrho$ follows from the same bootstrapping argument of [@greg_mckean_vlasov_torus Theorem 2.3].

We can also obtain an estimate on the rate of convergence to the equilibrium density in $L^2(U)$ as $t\to \infty$ with the following theorem. In order to forgo additional assumptions on the initial data $\varrho_0$ we restrict ourselves to the case where the equilibrium density is unique and given by $\varrho_\infty$.\

\[lem:exp\_conv\_in\_L2\]

Let $\varrho\in C^1([0,\infty];C^2(U)) $ be a solution of with initial data $\varrho_0\in L^2(U)$ a probability density. For $\kappa_1=0$, if $$\begin{aligned}
\kappa_2^2 < 
\min \Big\{  \frac{\mu_{\min}c_{pw}^{-2}\|\nabla_{\vec{r}}V_2\|^{-2}_{L^\infty}}{2(1+e)\mu_{\max}},
\frac{1}{4\|V_2\|_{L^\infty}} \Big\},\end{aligned}$$ where $c_{pw}$ is a Poincar[' e]{}$-$Wirtinger constant on the domain $U$ and $\mu_{\max}$ and $\mu_{\min}$ are the largest and smallest eigenvalues of $\bm{D}$, then $\varrho\to \varrho_\infty$ in $L^2(U)$ exponentially as $t\to \infty$. In particular the convergence in $L^2(U)$ is given by $$\begin{aligned}
\|\varrho(\cdot, t)-\varrho_\infty(\cdot)\|^2_{L^2(U)}\leq \|\varrho_0(\cdot)-\varrho_\infty(\cdot)\|_{L^2(U)}^2 e^{-r_{\kappa_2}t}\end{aligned}$$ as $t \to \infty$ where $r_{\kappa_2} = \mu_{\min}c^{-2}_{pw}- 2\mu_{\max}|\kappa_2|^2(e+1)\|\nabla_{\vec{r}}V_2\|_{L^{\infty}(U)}^2$ is the rate of convergence.

Let $\psi  = \varrho-\varrho_\infty$, then the evolution equation for $\psi$ may be written $$\begin{aligned}
\label{eq:exp_conv_in_l2_bound_1}
\partial_t\psi - \nabla_{\vec{r}}\cdot[\bm{D}\,\nabla_{\vec{r}}\psi] = \kappa_2 \nabla_{\vec{r}}\cdot[\bm{D}\,(\varrho_\infty\nabla_\vec{r}V_2\star \psi +\psi\nabla_{\vec{r}}V_2\star \varrho)].\end{aligned}$$ Multiplying by $\psi$, integrating and using the boundary condition $\Pi[\psi]\cdot\vec{n} = 0$ on $\partial U\times [0,T]$ we obtain $$\begin{aligned}
&\tfrac{1}{2}\der[]{t}\|\psi (t)\|^2_{L^2(U)} +\|\bm{D}^{1/2}\nabla_{\vec{r}}\psi \|_{L^2(U)}^2 \nonumber\\
&\leq \int_U\mathrm{d}\vec{r}\,|\bm{D}^{1/2}\nabla_{\vec{r}}\psi |\times |\kappa_2 \bm{D}^{1/2}(\varrho_\infty\nabla_\vec{r}V_2\star \psi +\psi\nabla_{\vec{r}}V_2\star \varrho)|. \end{aligned}$$ Using H[ö]{}lder’s inequality on the right hand side this becomes $$\begin{aligned}
&\tfrac{1}{2}\der[]{t}\|\psi (t)\|^2_{L^2(U)} +\|\bm{D}^{1/2}\nabla_{\vec{r}}\psi \|_{L^2(U)}^2 \nonumber\\
&\leq \|\bm{D}^{1/2}\nabla_{\vec{r}}\psi \|_{L^2(U)}\times \|\kappa_2 \bm{D}^{1/2}(\varrho_\infty\nabla_\vec{r}V_2\star \psi +\psi\nabla_{\vec{r}}V_2\star \varrho)\|_{L^2(U)}. \end{aligned}$$ Now using Young’s inequality twice on the right hand side we obtain $$\begin{aligned}
&\tfrac{1}{2}\der[]{t}\|\psi (t)\|^2_{L^2(U)} +\|\bm{D}^{1/2}\nabla_{\vec{r}}\psi \|_{L^2(U)}^2\nonumber\\
& \leq \tfrac{1}{2}\|\bm{D}^{1/2}\nabla_{\vec{r}}\psi \|_{L^2(U)}^2 
+  \tfrac{1}{2} \|\bm{D}^{1/2}(\varrho_\infty\nabla_\vec{r}V_2\star \psi +\psi\nabla_{\vec{r}}V_2\star \varrho)\|_{L^2(U)}^2 \nonumber \\
& \leq \tfrac{1}{2}\|\bm{D}^{1/2}\nabla_{\vec{r}}\psi \|_{L^2(U)}^2 +  |\kappa_2|^2 \|\varrho_\infty\bm{D}^{1/2}\nabla_{\vec{r}}V_2\star \psi \|^2_{L^2(U)} + |\kappa_2|^2 \|\psi\bm{D}^{1/2}\nabla_{\vec{r}}V_2\star \varrho \|^2_{L^2(U)}.\label{eq:l2_trend_inq_1}\end{aligned}$$ From the positive definiteness and boundedness of the diffusion tensor, we have $\mu_{\min}\leq \|\bm{D}\|_{L^\infty(U)}\leq \mu_{\max}$.

We also have the following bounds in terms of $\|\psi\|_{L^2(U)}^2$ $$\begin{aligned}
&\|\psi\bm{D}^{1/2}\nabla_{\vec{r}}V_2\star \varrho \|^2_{L^2(U)}\leq\mu_{\max}\|\nabla_{\vec{r}}V_2\|_{L^\infty(U)}^2\|\psi\|_{L^2(U)}^2\label{eq:l2_trend_formula_1a}\\
&\|\varrho_\infty\bm{D}^{1/2}\nabla_{\vec{r}}V_2\star \psi \|^2_{L^2(U)}\leq |U|\mu_{\max}\|\varrho_\infty\|_{L^2(U)}^2\|\nabla_{\vec{r}}V_2\|^2_{L^\infty(U)}\|\psi \|_{L^2(U)}^2\label{eq:l2_trend_formula_1b}\end{aligned}$$ where $|U|$ denotes the size of $U$ and in we have used that $\nabla V_2 \star \varrho \leq  \| \nabla V_2\|_{L^\infty(U)} \| \varrho \|_{L^1(U)}$ and the fact that $\varrho$ is a probability density with $\| \varrho \|_{L^1} =1$ (see Corollary \[cor:L\_1\_varrho\_is\_1\]). To obtain we use that $$\|\varrho_\infty\bm{D}^{1/2}\nabla_{\vec{r}}V_2\star \psi \|^2_{L^2(U)}\leq
    \mu_{\max}\|\varrho_\infty\|_{L^2(U)}^2\|\nabla_{\vec{r}}V_2\|^2_{L^\infty(U)}
    \int_U \mathrm{d}\vec{r}\, \Big| \rho_\infty(\vec{r}) \int \mathrm{d}\vec{r} ' \, \psi(\vec{r}') \Big|^2.$$ We then note that, by Hölder’s inequality, $ \int \mathrm{d}\vec{r} ' \, \psi(\vec{r}') \leq \|\psi\|_{L^2}\|1\|_{L^2(U)} = 
|U|^{1/2} \|\psi\|_{L^2}$, which gives the result. For it remains to bound the non explicit stationary distribution $\varrho_\infty$ in $L^2(U)$, to do this we observe that by the self-consistency equation $$\begin{aligned}
\|\varrho_\infty\|_{L^2(U)}^2\leq \frac{|U|\times e^{2|\kappa_2\||V_2\|_{L^\infty}}}{|U|^2\times e^{-2|\kappa_2\||V_2\|_{L^\infty}}}.\label{eq:l2_trend_formula_2}\end{aligned}$$

Using , and the bounds on $\bm{D}$, inequality becomes $$\begin{aligned}
\tfrac{1}{2}\der[]{t}\|\psi (t)\|^2_{L^2(U)}&\leq -\tfrac{\mu_{\min}}{2}\|\nabla_{\vec{r}}\psi \|_{L^2(U)}^2\nonumber\\
&\quad+ \mu_{\max}|\kappa_2|^2(e^{4|\kappa_2\||V_2\|_{L^\infty}}+1)\|\nabla_{\vec{r}}V_2\|_{L^{\infty}(U)}^2\|\psi \|^2_{L^2(U)}. \end{aligned}$$ Now since $\psi$ has mean zero we may use the Poincar[' e]{}–Wirtinger inequality to write $$\begin{aligned}
\der[]{t}\|\psi (t)\|^2_{L^2(U)}&\leq -\mu_{\min}c^{-2}_{pw}\|\psi \|_{L^2(U)}^2\nonumber\\
&\quad+ 2\mu_{\max}|\kappa_2|^2(e^{4|\kappa_2\||V_2\|_{L^\infty}}+1)\|\nabla_{\vec{r}}V_2\|_{L^{\infty}(U)}^2\|\psi \|^2_{L^2(U)}.\end{aligned}$$ Finally, by Gr[ö]{}nwall’s lemma [@evans2002partial], we obtain $$\begin{aligned}
&\|\psi (t)\|^2_{L^2(U)}\nonumber\\
&\leq \|\psi(0)\|_{L^2(U)}^2\exp\left\lbrace -(\mu_{\min}c^{-2}_{pw}- 2\mu_{\max}|\kappa_2|^2(e^{4|\kappa_2\||V_2\|_{L^\infty}}+1)\|\nabla_{\vec{r}}V_2\|_{L^{\infty}(U)}^2)t\right\rbrace.  \label{eq:l2_trend_gronwall_1}\end{aligned}$$

Therefore for any $\varrho_\ast$ a stationary density the necessary condition for exponential convergence $\varrho\to\varrho_\ast$ in $L^2(U)$ as $t\to\infty$ is $$\label{eq:gen_kappa_2_ineq_for_exp_convergence}
\mu_{\min}c^{-2}_{pw}- 2\mu_{\max}|\kappa_2|^2(e^{4|\kappa_2\||V_2\|_{L^\infty}}+1)\|\nabla_{\vec{r}}V_2\|_{L^{\infty}(U)}^2>0.$$

It will now be seen that, under the assumption that $\varrho_\infty$ is the unique stationary density with $\kappa_2 \leq \|V_2\|_{L^\infty}^{-1}/4$, we may obtain an explicit condition for $|\kappa_2|$. In particular becomes $$\begin{aligned}
\|\psi (t)\|^2_{L^2(U)}\leq \|\psi(0)\|_{L^2(U)}^2 \exp\left\lbrace -(\mu_{\min}c^{-2}_{pw}- 2\mu_{\max}|\kappa_2|^2(e+1)\|\nabla_{\vec{r}}V_2\|_{L^{\infty}(U)}^2)t\right\rbrace.  \end{aligned}$$\[eq:l2\_trend\_gronwall\_2\] Then to ensure the argument in the exponential remains negative, we require $$\begin{aligned}
|\kappa_2|^2<\frac{\mu_{\min}c_{pw}^{-2}\|\nabla_{\vec{r}}V_2\|^{-2}_{L^\infty}}{2(1+e)\mu_{\max}}.\end{aligned}$$ This completes the proof of the theorem.

We remark that $\psi\in \left\lbrace u\in H^1(U)\, | \,\int_U\mathrm{d}\vec{r}\,u = 0 \right\rbrace $, therefore, we may determine that the sharpest value of $c_{pw}$ conincides with the Poincar[é]{} constant as found by Steklov [@kuznetsov2015sharp], equal to $\nu_1^{-1/2}$ where $\nu_1$ is the smallest eigenvalue of the problem $$\begin{aligned}
\Delta u &= -\nu u \quad \text{ in } U,\\
\partial_{\vec{n}}u & = 0 \qquad \,\,\,\,  \text{ on } \partial U.\end{aligned}$$ Here $\partial_{\vec{n}}$ is the directional derivative along the unit vector $\vec{n}$ pointing out of the domain $U$. Additionally Payne and Weinberger [@payne1960optimal] proved that for convex domains in $\mathbb{R}^n$ one has $c_{pw}\leq \tfrac{\mathrm{diam}(U)}{\pi}$.

One may obtain a similar convergence result including a confining potential as given by the following corollary.\

\[prop:general\_kappa1\_convergence\]

Let $\kappa_1\neq 0$ and let $\varrho\in C^1([0,\infty];C^2(U)) $ be a solution of with initial data $\varrho_0\in L^2(U)$ a probability density. Then the exponential convergence $\varrho\to\varrho_\infty$ in $L^2$ criteria is modified to $$\begin{aligned}
\mu_{\max}\kappa_1^2\|\nabla_{\vec{r}}V_1\|_{L^\infty(U)}^2<r_{\kappa_2}\end{aligned}$$ along with $| \kappa_2 | \leq 1/4\times \|V_2\|_{L^\infty(U)}^{-1}$. In particular the convergence in $L^2$ is given by $$\begin{aligned}
\|\varrho(\cdot, t)-\varrho_\infty(\cdot)\|^2_{L^2(U)}\leq \|\varrho_0(\cdot)-\varrho_\infty(\cdot)\|_{L^2(U)}^2 e^{-(r_{\kappa_2}-\mu_{\max}\kappa_1^2\|\nabla_{\vec{r}}V_1\|_{L^\infty(U)}^2)t}.\end{aligned}$$

Since the inclusion of an external field is linear in the PDEs , the proof is similar to Lemma \[lem:exp\_conv\_in\_L2\], the only term to resolve for the evolution equation for $\psi$ first occurring at being $$\begin{aligned}
\kappa_1^2\|\psi\boldsymbol{D}^{1/2}\nabla_{\vec{r}}V_1\|_{L^2(U)}^2 \leq \kappa_1^2\mu_{\max}\|\nabla V_1\|_{L^\infty(U)}^2\|\psi\|_{L^2(U)}^2.\end{aligned}$$ The remainder of the calculations to derive a Gr[ö]{}nwall type inequality including this term are similar.

Global Asymptotic Stability {#sec:global_asymptotic_stability}
===========================

In this section we study the stability properties of stationary states. We start by showing the free energy is a strictly convex functional, provided $\kappa_2$ is sufficiently small, and that $\mathcal{F}$ is bounded below. Recall the free energy functional $\mathcal{F}:P^+_{\text{ac}}(U)\to \mathbb{R}$ is given by $$\begin{aligned}
\mathcal{F}[\varrho]&:=\int_U\mathrm{d}\vec{r}\,\varrho(\vec{r})\log\varrho(\vec{r})+\kappa_1\int_U\mathrm{d}\vec{r}\,V_1(\vec{r})\varrho(\vec{r})+\frac{\kappa_2}{2}\int_U\mathrm{d}\vec{r}\,\varrho(\vec{r})V_2\star\varrho(\vec{r}). \nonumber\end{aligned}$$\

\[prop:F\_is\_strictly\_convex\]

For $|\kappa_2|\in [0,\|V_2\|_{L^\infty(U)}^{-1})$ the free energy functional $\mathcal{F}$ is strictly convex. Additionally there exists a positive constant $B_0<\infty$ for every $\varrho\in P_{\text{ac}}^+$ such that $|\mathcal{F}[\varrho]| \geq B_0$.

Suppose $\varrho_1$ and $\varrho_2$ satisfy with $\Pi[\varrho_1]\cdot\vec{n} = \Pi[\varrho_2]\cdot\vec{n} = 0$ on $\partial_U$ for all $t\in [0,\infty)$. Letting $\zeta = \varrho_2-\varrho_1$ and $\varrho_s = (1-s)\varrho_1+s\varrho_2$ we compute $\tfrac{\mathrm{d}^2}{\mathrm{d}s^2}\mathcal{F}_H[\varrho_s]$ by direct calculation $$\begin{aligned}
\frac{\mathrm{d}^2}{\mathrm{d}s^2}\mathcal{F}_H[\varrho_s] &= \frac{\mathrm{d}}{\mathrm{d}s}\frac{\mathrm{d}}{\mathrm{d}s}\left[\int_U\mathrm{d}\vec{r}\,\varrho_s\log\varrho_s +\kappa_1 \int_U\mathrm{d}\vec{r}\,\varrho_s V_1 +\frac{\kappa_2}{2} \int_U\mathrm{d}\vec{r}\,\varrho_s V_2\star\varrho_s\right]\nonumber\\
&= \frac{\mathrm{d}}{\mathrm{d}s}\left[\int_U\mathrm{d}\vec{r}\,\zeta \log\varrho_s +\zeta \right.\nonumber\\
& \qquad \left.  + \kappa_1\int_U\mathrm{d}\vec{r}\,\zeta V_1 + \frac{\kappa_2}{2}\int_U\mathrm{d}\vec{r}\,\zeta V_2\star\varrho_s + \frac{\kappa_2}{2}\int_U\mathrm{d}\vec{r}\,\varrho_s V_2\star\zeta\right]\nonumber\\
&=\int_U\mathrm{d}\vec{r}\,\frac{\zeta^2}{\varrho_s} +  \kappa_2\int_U\mathrm{d}\vec{r}\,\zeta V_2\star \zeta.\end{aligned}$$

Now using the measure $\mathrm{d}\mu = \varrho_s\mathrm{d}\vec{r}$ we have, by Jensen’s inequality, $$\begin{aligned}
\int_U\mathrm{d}\vec{r}\,\frac{\zeta^2}{\varrho_s}  = \int_U\mathrm{d}\mu\,\frac{\zeta^2}{\varrho_s^2} \geq \left(\int_U\mathrm{d}\vec{r}\,|\zeta|\right)^2.\end{aligned}$$ We also have that $V_2$ is bounded below by the negative of its its essential supremum from . Combining these facts we find $$\begin{aligned}
\label{eq:convexity_condition_of_F}
\frac{\mathrm{d}^2}{\mathrm{d}s^2}\mathcal{F}_H[\varrho_s] \geq (1-|\kappa_2\|V_2\|_{L^\infty(U)}\})\left(\int_U\mathrm{d}\vec{r}\,|\zeta|\right)^2\end{aligned}$$ and we therefore find that, for $\kappa_2$ such that $|k_2|\leq \tfrac{1}{4}\|V_2\|_{L^\infty(U)}^{-1}$, the free energy functional $\mathcal{F}$ is strictly convex.

Now let $\varrho\in P^+_{\text{ac}}$ and observe that $$\begin{aligned}
\mathcal{F}[\varrho]\geq -\Big|\int_U\mathrm{d}\vec{r}\,\varrho\log\varrho\Big| -|\kappa_1|\|V_1\|_{L^\infty(U)}-\tfrac{|\kappa_2|}{2}\|V_2\|_{L^\infty(U)}.\end{aligned}$$ The entropy $\varrho\log\varrho$ is continuous and bounded below on $U$ and therefore we have that $$\begin{aligned}
\mathcal{F}[\varrho]\geq \int\mathrm{d}\vec{r}\,|\varrho\log\varrho|-|\kappa_1|\|V_1\|_{L^\infty(U)}-\tfrac{|\kappa_2|}{2}\|V_2\|_{L^\infty(U)}>-\infty.\end{aligned}$$ where we have used the assumptions on the potentials in . Hence $\mathcal{F}[\cdot]$ is bounded below.

Note that the convexity condition in Proposition \[prop:F\_is\_strictly\_convex\] holds independently of the confining potential $V_1$. We therefore have the following Corollary for the total free energy $\mathcal{F}-\int_U\mathrm{d}\vec{r}\,\mu_{c}\varrho$.\

\[cor:the\_free\_energy\_is\_convex\_and\_bounded\_below\] The total free energy $\mathcal{F}-\int\mathrm{d}\vec{r}\,\mu_{c}\varrho$ is strictly convex for $|\kappa_2|\in [0,\|V_2\|_{L^\infty(U)}^{-1})$ and bounded below.

We now provide a useful Lemma which will be used eventually to show that $\mathcal{F}$ always has a minimiser, for any $\kappa_2$ (see Lemma \[lem:minimisers\_always\_exist\]).\

\[lem:can\_pick\_a\_smaller\_free\_energy\_density\] Let $V_1$, $V_2$ satisfy the assumptions then there exists a positive constant $B_0$ such that for every $\varrho\in P_{ac}(U)$ with $\|\varrho\|_{L^\infty(U)}>B_0$ there exists some $\varrho^\dagger \in P_{ac}(U)$ with $\|\varrho^{\dagger}\|_{L^\infty(U)}\leq B_0$ such that $$\begin{aligned}
\mathcal{F}(\varrho^\dagger)<\mathcal{F}(\varrho).\end{aligned}$$

For a proof see [@greg_mckean_vlasov_torus Lemma 2.5] or [@chayes2010mckean Lemma 2.1], the only modification required is to include $V_1$ which by assumption is bounded below and the proof follows a similar argument.

We now show that minimisers of $\mathcal{F}$ exist for all $\kappa_2$. First we define the integral operator $\mathcal{R}$ which will be useful for the following calculations.\

\[def:def\_of\_R\_op\] Let $\mathcal{R}:L^1(U)\to L^1(U)$ be given by $$\begin{aligned}
\label{eq:def_of_R_operator}
 \mathcal{R}u = -\varrho V_2\star u.
 \end{aligned}$$ We note that $\mathcal{R}$ is a compact (since $V_2$ is uniformly bounded in $L^\infty(U)$) self-adjoint operator in $L^2(U,\varrho^{-1})$. We label its eigenvalues $\{\beta_n^{-1}\}_{n=1}^\infty$ and eigenfunctions $\{u_n\}_{n=1}^\infty$ satisfying $$\begin{aligned}
\label{eq:eigenvalue_problem_for_R}
 \mathcal{R}u_n = \beta_n^{-1}u_n.\end{aligned}$$

\[lem:minimisers\_always\_exist\] Let $\kappa_2\in (-\infty,\infty)$ and let $V_1$, $V_2$ satisfy the assumptions . Then there exists a $\varrho\in P_{ac}(U)$ that minimizes $\mathcal{F}$.

Since $\mathcal{F}$ is bounded below there exists a minimising sequence $\{\varrho_j\}_{j=1}^\infty \in P_{ac}(U)$ so that $\mathcal{F}(\varrho_{j})<\mathcal{F}(\varrho_{j+1})$. Therefore, by Lemma \[lem:can\_pick\_a\_smaller\_free\_energy\_density\] $\{\varrho_j\}_{j = 1}^\infty$ may be chosen such that $\|\varrho_j\|_{L^2}(U)\leq \|\varrho_j\|^2_{L^\infty(U)}|U|$. Now by the Eberlein-Smuljan theorem, since $\{\varrho_j\}_{j = 1}^\infty$ is bounded, there exists a subsequence (which we will denote again by $\{\varrho_j\}_{j = 1}^\infty$) such that $\varrho_j  \rightharpoonup \varrho_\ast$ weakly in $L^2$ to some $\varrho_\ast$. Therefore $\lim_{j\to \infty}\int_U\mathrm{d}\vec{r}\,\eta(\varrho_j-\varrho_\ast) = 0$ for every $\eta\in L^2(U)$, so in particular for $\eta = 1$ we obtain $\lim_{j\to \infty}\int_U\mathrm{d}\vec{r}\, \varrho_j = 1 = \int_U\mathrm{d}\vec{r}\, \varrho_\ast$. Additionally we note that $|\varrho_j| \rightharpoonup |\varrho_\ast|$ in $L^2(U)$, and therefore $\|\varrho_\ast\|_{L^1(U)} = 1$, which is enough to show that $\varrho_\ast\geq 0 $ a.e. by standard arguments (see, for example, the proof of Corollary \[cor:L\_1\_varrho\_is\_1\]).

We define $\Lambda:P_{ac}\to \mathbb R$ such that $$\begin{aligned}
\Lambda(z) := \int_U\mathrm{d}\vec{r}\,zV_2\star z.\end{aligned}$$ Now let $\varrho_{\beta_{n}}\in L^1(U)$ be a solution to , which is known to exist by Lemma \[thm:exis\_fix\_point\]. Note that $\varrho_{\beta_{n}}$ need not be a minimiser of $\mathcal{F}$ and may be an inflection point or local maximum. Additionally since $\varrho_{\beta_{n}}\in L^1(U)$ solves , we have that $\varrho_{\beta_{n}}>e^{-(|\kappa_1|\|V_1\|_{L^\infty(U)}+|\beta_{n}|\|V_2\|_{L^\infty(U)})}/Z>0$ (where $Z$ is a normalisation constant) and therefore there exists $\delta\in \mathbb{R}^+$ such that $\varrho_{\beta_{n}}>\delta$ for every $\vec{r}\in U$.

Now we estimate the interaction energy difference by $$\begin{aligned}
|\Lambda(\varrho_j)-\Lambda(\varrho_\ast)| &\leq \sum_{n=1}^N|\beta_n^{-1}|\Big|\langle \varrho_j,w_n\rangle_{L^2(U,\varrho^{-1}_{\beta_{n}})}-\langle \varrho_\ast,w_n\rangle_{L^2(U,\varrho^{-1}_{\beta_{n}})}\Big| + 2 |\beta_N^{-1}|\delta^{-1}B_0\nonumber\\
& \leq 2\delta^{-1} B_0  \sum_{n=1}^N \langle\varrho_j -\varrho_\ast, w_n\rangle_{L^2(U)} + 2 |\beta_N^{-1}|\delta^{-1} B_0\end{aligned}$$ where we have used the fact that the integrand of $\Lambda(z)$ is equal to $\mathcal{R}$ acting on $z \in P_{ac}$. Additionally we have used that $\mathcal{R}$ is self-adjoint in $L^2(U,\varrho^{-1}_{\beta_{n}})$, to write $\mathcal{R}$ as a projection onto its eigenvectors $\{w_n\}_{n=1}^{\infty}$ and bounded the tail of the infinite sum using Bessel’s inequality.

Now since $\mathcal{R}$ is self-adjoint in $L^2(U, \varrho_{\beta_n}^{-1})$ we have that (after reordering) $|\beta_n^{-1}|\to 0$ as $n \to \infty$ so the second term may be made arbitrarily small. The first term may be made arbitrarily small by taking the limit $j \to \infty$ inside the finite sum and using that $\varrho_j\rightharpoonup \varrho_{\ast}$ weakly in $L^2(U)$. This shows that $\Lambda(\cdot)$ is continuous in $\varrho$.

Additionally, for the external energy, we have $$\begin{aligned}
\Big|\int_U\mathrm{d}\vec{r}\, V_1(\vec{r})\varrho_j(\vec{r})-\int_U\mathrm{d}\vec{r}\, V_1(\vec{r})\varrho_\ast(\vec{r})\Big| & = \Big|\int_U\mathrm{d}\vec{r}\, V_1(\vec{r})(\varrho_j(\vec{r})-\varrho_\ast(\vec{r}))\Big|  \nonumber\\
& \leq  \Big|\int_U\mathrm{d}\vec{r}\, V_1(\vec{r})(\varrho_j(\vec{r})-\varrho_\ast(\vec{r}))\Big| \to 0\end{aligned}$$ as $j\to \infty$. The lower semicontinuity of the entropy term in follows from standard results [@jost1998calculus Lemma 4.3.1]. Therefore the free energy $\mathcal{F}[\varrho]$ has a minimiser $\varrho$ over $P_{ac}(U)$.

We may refine this result to show that minimisers are attained in $P^+_{ac}(U)$ with the following lemma.\

\[lem:minimisers\_are\_positive\] Let $\varrho\in P_{ac}(U)\backslash P^+_{ac}(U)$. Then there exists $\varrho^\dagger\in P^+_{ac}(U)$ such that $\mathcal{F}[\varrho^\dagger]<\mathcal{F}[\varrho]$.

The proof is similar to [@greg_mckean_vlasov_torus Lemma 2.6]. One must show that the potential energy for a $P^+_{ac}(U)$ density may be bounded by the potential energy of a $ P_{ac}(U)$ density. We let $\epsilon>0$ and define the competition state $\varrho_\epsilon$ such that $$\begin{aligned}
\varrho_\epsilon(\vec{r}) = \frac{(\varrho(\vec{r})+\epsilon\mathbb{I}_{\mathbb{B}_0}(\vec{r}))}{1+\epsilon |\mathbb{B}_0|}\end{aligned}$$ where $\mathbb{B}_0 = \{\vec{r}\in U \, :\, \varrho(\vec{r}) = 0\}$ and since by assumption $\varrho \notin  P^+_{ac}(U)$ one has $|\mathbb{B}_0|>0$ and $\varrho_\epsilon\in P^+_{ac}(U)$. Then we obtain that $$\begin{aligned}
\int_U\mathrm{d}\vec{r}\, V_1 \varrho_\epsilon \leq \int_U\mathrm{d}\vec{r}\, V_1 \varrho + \epsilon |\mathbb{B}_0|.\end{aligned}$$ Using this bound, together with the result [@greg_mckean_vlasov_torus Lemma 2.6] we obtain the required result.

Exponential Convergence to Equilibrium in Relative Entropy.
-----------------------------------------------------------

In this section we derive an H-theorem which guarantees that the time evolution of the dynamics converges to the equilibrium distribution given by the self-consistency equation. First consider the time derivative of the integral of the free energy $$\begin{aligned}
&\der[]{t} \mathcal{F}[\varrho] = \int_U \mathrm{d}\vec{r}\, \partial_t \varrho\, \tfrac{\delta\mathcal{F}[\varrho] }{\delta \varrho} = \int_U \mathrm{d}\vec{r}\,\nabla\cdot\left[\bm{D}(\vec{r},t)\varrho(\vec{r},t)\nabla_{\vec{r}}\tfrac{\delta\mathcal{F}[\varrho] }{\delta \varrho} \right]\tfrac{\delta\mathcal{F}[\varrho] }{\delta \varrho} \nonumber\\
&=-\int_U \mathrm{d}\vec{r}\,\Big|\bm{D}(\vec{r},t)^{1/2}\varrho(\vec{r},t)^{1/2}\nabla_{\vec{r}}\tfrac{\delta\mathcal{F}[\varrho] }{\delta \varrho} \Big|^2,\end{aligned}$$ where we have integrated by parts and used the boundary condition $\Pi\varrho \cdot\vec{n}\big|_{\partial_U}= 0$ or $\varrho\to 0$ as $|\vec{r}|\to\infty$ for bounded and unbounded domains respectively. Here we see that as long as both $\bm{D}(\vec{r},t)$ and $\varrho(\vec{r},t)$ remain positive definite then the free energy is monotonically decreasing in time. Indeed the diffusion tensor $\bm{D}$ is positive definite as proven in [@goddard2012overdamped] and we will show strict positivity of $\varrho(\vec{r},t)$ in Section \[subsec:strict\_pos\_rho\]. We now introduce the relative entropy functional $$\begin{aligned}
\label{eq:def_of_rel_entropy_func}
\mathcal{H}[\varrho |\varrho_\infty]:=\int_{U}\mathrm{d}\vec{r}\, \varrho \log\left(\frac{\varrho}{\varrho_\infty}\right),\end{aligned}$$ and obtain the following theorem for convergence to equilibrium in relative entropy.

\[thm:rel\_entropy\_convergence\]

Let $V_1$ be convex, $|\kappa_2|<\tfrac{1}{4}\|V_2\|_{L^\infty(U)}^{-1}$ and $\varrho\in C^1([0,\infty];C^2(U))$ be a classical solution to equation . If $\kappa_2^2 < \frac{c_{ls}^{-1}}{2 \| \nabla V_2 \|_{L^\infty(U)}^2}$ then $\varrho$ is exponentially stable in relative entropy and it holds that $$\begin{aligned}
\mathcal{H}[\varrho |\varrho_\infty]\leq \mathcal{H}[\varrho_0|\varrho_\infty] e^{-\tfrac{1}{2}(c_{ls}^{-1}-2|\kappa_2|^2 \|\nabla V_2\|_{L^\infty(U)}^2)t},\end{aligned}$$ where $c_{ls}>0$ is the log-Sobolev constant for the measure $\mu$.

By direct calculation we find $$\begin{aligned}
\label{eq:H_time_der}
\frac{\mathrm{d} \mathcal{H}[\varrho |\varrho_\infty]}{\mathrm{d}t} &= \int_U\mathrm{d}\vec{r}\, \partial_t\left( \varrho\log\left(\frac{\varrho}{\varrho_\infty}\right)\right) =   \int_U\mathrm{d}\vec{r}\, \partial_t\varrho \log\left(\frac{\varrho}{\varrho_\infty}\right) + \int_U\mathrm{d}\vec{r}\,\partial_t \varrho\nonumber\\
& = \int_U\mathrm{d}\vec{r}\, \partial_t\varrho \log\left(\frac{\varrho}{\varrho_\infty}\right) + \frac{\mathrm{d}M}{\mathrm{d}t} =  - \int_U\mathrm{d}\vec{r}\, \varrho\nabla \frac{\delta \mathcal{F}[\varrho]}{\delta \varrho}\cdot \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right) + 0 \nonumber\\
&= - \int_U\mathrm{d}\vec{r}\, \varrho \left(\nabla \log\varrho + \kappa_1 \nabla V_1 +\kappa_2\nabla V_2\star \varrho\right)\cdot   \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right)\nonumber\\
& = - \int_U\mathrm{d}\vec{r}\, \varrho \left(\nabla \log\varrho +\kappa_2\nabla V_2\star \varrho-\left( \nabla \log \varrho_\infty +\kappa_2\nabla V_2\star \varrho_\infty\right)\right)\cdot   \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right)\nonumber\\
& =- \int_U\mathrm{d}\vec{r}\, \varrho \left(\nabla \log\left(\frac{\varrho}{\varrho_\infty}\right) +\kappa_2\nabla V_2\star (\varrho-\varrho_\infty)\right)\cdot   \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right)\end{aligned}$$ where we have used the no-flux boundary condition and the self-consistency equation $\nabla \log\varrho_\infty + \kappa_1 \nabla V_1 +\kappa_2\nabla V_2\star \varrho_\infty = 0$. Note that the contribution from the $V_1$ term is constant, independent of $\rho$, and so cancels after using the self-consistency equation.

Continuing by expanding out the integrand and using H[ö]{}lder’s inequality we obtain $$\begin{aligned}
\frac{\mathrm{d} \mathcal{H}[\varrho |\varrho_\infty]}{\mathrm{d}t} 
&= - \int_U\mathrm{d}\vec{r}\, \varrho \Big| \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right) \Big|^2 + \kappa_2 \int_U\mathrm{d}\vec{r}\, \varrho \nabla\log \left(\frac{\varrho}{\varrho_\infty}\right)  \cdot \nabla V_2\star (\varrho- \varrho_\infty)\nonumber\\
&\leq -\int_U\mathrm{d}\vec{r}\, \varrho \Big| \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right)\Big|^2\nonumber\\
&\qquad + \left[\int_U\mathrm{d}\vec{r}\, \varrho \Big|\nabla\log \left(\frac{\varrho}{\varrho_\infty}\right)\Big|^2\right]^{1/2} 
 \times  \left(\kappa_2^2 \int_U\mathrm{d}\vec{r}\, \varrho |\nabla V_2\star (\varrho-\varrho_\infty)|^2\right)^{1/2}.\end{aligned}$$ Now, by Young’s inequality, $$\begin{aligned}
\frac{\mathrm{d} \mathcal{H}[\varrho |\varrho_\infty]}{\mathrm{d}t} \leq -\tfrac{1}{2}\int_U\mathrm{d}\vec{r}\, \varrho \Big| \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right)\Big|^2 +\tfrac{\kappa_2^2}{2} \int_U\mathrm{d}\vec{r}\, \varrho |\nabla V_2\star (\varrho-\varrho_\infty)|^2\end{aligned}$$ and we may estimate the second term on the right hand side (in particular using that $\int_U \rho = 1$ from Corollary \[cor:L\_1\_varrho\_is\_1\]), giving $$\begin{aligned}
\label{eq:rel_entropy_proof_first_bound}
\frac{\mathrm{d} \mathcal{H}[\varrho |\varrho_\infty]}{\mathrm{d}t} \leq -\tfrac{1}{2}\int_U\mathrm{d}\vec{r}\, \varrho \Big| \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right)\Big|^2 + \tfrac{\kappa_2^2}{2}\|\nabla V_2\|_{L^\infty(U)}^2\|\varrho-\varrho_\infty \|_{L^1(U)}^2\end{aligned}$$

We bound the first term as follows. Since $V_1$ is convex, we have $$\begin{aligned}
\nabla_{\vec{r}}^2V_1\geq \theta_1>0\end{aligned}$$ for some $\theta_1\in \mathbb{R}^+$. Now by the Bakry–[É]{}mery criterion (see [@menz2014poincare Sec 3, Theorem 3.1], and [@malrieu2001logarithmic]) the measure $\mu'(\mathrm{d}\vec{r}) = \mathrm{d}\vec{r}\,e^{-\kappa_1 V_1}/Z$ where $Z$ is a normalisation constant satisfies a log-Sobolev inequality (LSI) with constant $c_{ls}'$ such that $$\begin{aligned}
\frac{1}{c_{ls}'}\geq \theta_1\kappa_1.\end{aligned}$$ However since $V_2$ is not general a convex function, we cannot use the Bakry–[É]{}mery criterion for $\mu$ as defined in . However we may deduce a LSI using the Holley–Stroock perturbation lemma [@menz2014poincare Sec 3, Theorem 3.2] since $V_1+V_2\star\varrho_\infty$ is a bounded perturbation of $V_1$, in particular $$\begin{aligned}
\Big|V_1+V_2\star\varrho_\infty\Big| \leq \Big|V_1\Big|+ \|V_2\|_{L^\infty}\|\varrho_\infty\|_{L^1(U)} <\infty.\end{aligned}$$ Therefore $\mu$ as defined in with $\varrho = \varrho_\infty$ (after appropriate nondimensionalisation) is unique and satisfies a LSI with constant $$\begin{aligned}
c_{ls}^{-1} \geq \exp\left(-\kappa_1\kappa_2\,\text{Osc}\left[V_2\star \varrho_\infty\right]\right)\frac{1}{c_{ls}'}\end{aligned}$$ where $$\begin{aligned}
\text{Osc}\left[V_2\star \varrho_\infty\right]   = \sup V_2\star \varrho_\infty - \inf V_2\star \varrho_\infty.\end{aligned}$$

The constant $c_{ls}$ is such that such that for each $f:U\to \mathbb{R}^+$ one has $$\begin{aligned}
\label{eq:LSI}
    \int_U f^2\log f^2\mathrm{d}\mu -\int_U f^2 \log \left(\int_U f^2 \mathrm{d}\mu\right)\mathrm{d}\mu  \leq c_{ls}\int_U |\nabla f|^2\mathrm{d}\mu = c_{ls} \int_U f^2 | \nabla \log f^2 |^2 \mathrm{d} \mu.\end{aligned}$$ We let $f = \sqrt{\varrho/\varrho_\infty}$ and $\mathrm{d}\mu = \varrho_\infty \mathrm{d}\vec{r}$ and the second term on the left hand side of is zero (since, again $\int_U \rho = 1$). Hence this shows that $$\mathcal{H}[\varrho |\varrho_\infty] = 
    \int_U f^2\log f^2\mathrm{d}\mu 
    \leq c_{ls} \int_U\mathrm{d}\vec{r}\, \varrho \Big| \nabla \log\left(\frac{\varrho}{\varrho_\infty}\right)\Big|^2.$$

We combine the LSI with Pinsker’s inequality [@bolley2005weighted] to deduce $$\begin{aligned}
\frac{\mathrm{d} \mathcal{H}[\varrho |\varrho_\infty]}{\mathrm{d}t}\leq -\tfrac{1}{2}(c_{ls}^{-1}-2 \kappa_2^2\|\nabla V_2\|_{L^\infty(U)}^2) 
\mathcal{H}[\varrho|\varrho_\infty].
$$ Thus we obtain, by Gr[ö]{}nwall’s inequality, $$\begin{aligned}
\mathcal{H}[\varrho |\varrho_\infty]\leq \mathcal{H}[\varrho_0|\varrho_\infty] \exp[ -\tfrac{1}{2}(c_{ls}^{-1}-2\kappa_2^2 \|\nabla V_2\|_{L^\infty(U)}^2) t ]\end{aligned}$$ and the theorem is proved.

The constant $c_{ls}$ is not known explicitly but may be estimated in terms of the convexity of $V_1$, $V_2$ and the curvature of $U$ [@chen1997estimates]. We now consider asymptotic expansions of the steady states for small interaction energy $\kappa_2$.\

Asymptotic Expansion of the Steady States For Weak Interactions.
----------------------------------------------------------------

We begin this section by recalling that steady states satisfy the self-consistency equation $$\begin{aligned}
\label{eq:self_consis}
\varrho = \frac{e^{-(\kappa_1 V_1+\kappa_2 V_2\star \varrho)}}{Z},\end{aligned}$$ where $Z = \int_U \mathrm{d}\vec{r}\,e^{-(\kappa_1V_1+\kappa_2V_2\star\varrho)}$. We know from Lemma \[thm:exis\_fix\_point\] that for sufficiently weak interactions, i.e. $|\kappa_2|<1/4\|V_2\|_{L^\infty(U)}^{-1}$, the stationary distribution is unique; equivalently, the nonlinear equation has a unique fixed point. Let $\kappa_2 \ll 1$, then the stationary solution $\varrho(\vec{r}) = \varrho_\infty(\vec{r})$ has the form $$\begin{aligned}
\varrho(\vec{r})  = \frac{e^{-\kappa_1 V_1(\vec{r})}}{Z(\varrho)}(1 + O(\kappa_2)),\end{aligned}$$ where the first order correction may be obtained explicitly as follows.

Recall the stationary equation for $\varrho$:

$$\begin{aligned}
\nabla_{\vec{r}} \cdot [\bm{D} ( \nabla_{\vec{r}}\varrho+\kappa_1\varrho\,\nabla_{\vec{r}} V_1(\vec{r}) +\kappa_2\varrho\,\nabla_{\vec{r}} V_2\star \varrho)]=0 & \text{ on } U, \label{eq:stat_pde_for_rho}\\
\bm{D}(\nabla_{\vec{r}}\varrho+\kappa_1\varrho\,\nabla_{\vec{r}}V_1 +\kappa_2\varrho\,\nabla_{\vec{r}} V_2\star \varrho)\cdot\vec{n}=0 & \text{ on } \partial_U. \label{eq:stat_pde_for_rho_bc}\end{aligned}$$

Fix $\kappa_1=1$ and insert the perturbation expansion $$\begin{aligned}
\varrho(\vec{r}) = \sum_{k=0}^\infty \kappa_2^k \varrho_k(\vec{r}).\end{aligned}$$ We find at the first order of $\kappa_2$

$$\begin{aligned}
\mathcal{L}_0\varrho_0:=\nabla_{\vec{r}} \cdot (\bm{D}\nabla_{\vec{r}}\varrho_0+\bm{D}\left(\varrho_0\nabla_{\vec{r}}V_1)\right)   = 0 & \text{ on } U,\label{eq:fredholm_1d_eqn}\\
\bm{D}(\nabla_{\vec{r}}\varrho_0+\varrho_0\nabla_{\vec{r}}V_1) \cdot \vec{n} =0 & \text{ on } \partial_U,\label{eq:fredholm_1d_eqn_bc}\end{aligned}$$

from which we deduce $$\begin{aligned}
\varrho_0(\vec{r}) = \frac{e^{-V_1(\vec{r})}}{Z_0}\end{aligned}$$ for $Z_0 = \int_{U}\mathrm{d}\vec{r}\,e^{-V_1(\vec{r})}$.

Note that $\mathcal{L}_0$ is self-adjoint in the space $L^2(U,\varrho_0^{-1})$. We may also show that the resolvent of $\mathcal{L}_0$ is compact in $L^2(U,\varrho^{-1}_0)$.\

\[lem:L0\_has\_compact\_resolvent\] The operator $\mathcal{L}_0$ has a compact resolvent in $L^2(U,\varrho_0^{-1})$.

We let $\phi\in C^2(U)$, by direct calculation we have that $$\begin{aligned}
\mathcal{L}_0\phi = [\nabla_{\vec{r}}\cdot \bm{D}]\cdot \nabla_{\vec{r}}\phi + \text{Tr}\left[\bm{D}\nabla_{\vec{r}}^2\phi \right] + [\bm{D}\nabla_{\vec{r}}V_1]\cdot \nabla_{\vec{r}}\phi + \phi\left[[\nabla_{\vec{r}}\cdot \bm{D}]\cdot \nabla_{\vec{r}}V_1+\text{Tr}[\bm{D}\nabla_{\vec{r}}^2V_1]\right]  \end{aligned}$$ then we have that $$\begin{aligned}
\|\mathcal{L}_0\phi\|_{L^2(U,\varrho_0^{-1})}^2   &= \int_U\mathrm{d}\vec{r}\, \Big|\nabla_{\vec{r}} \cdot (\bm{D}\nabla_{\vec{r}}\phi+\bm{D}\left(\phi\nabla_{\vec{r}}V_1)\right)\Big|^2\varrho_0^{-1}\nonumber\\
&\leq C(U;\bm{D};V_1)\sum_{n = 0}^2\sup_{\vec{r}\in U}\Big|\phi^{(n)}(\vec{r})\Big|<\infty\end{aligned}$$ where the constant $C(U;\bm{D};V_1)$ is dependent on $U$, the diffusion tensor $\bm{D}$ and the first weak derivatives of its entries (bounded in $L^{\infty}(U)$ by ), and the confining potential $V_1$ and its first two weak derivatives (bounded in $L^{\infty}(U)$ by ).

Therefore there exists $C\in \mathbb{R}^+$ such that $\|\mathcal{L}_0\|_{L^2(U,\varrho^{-1}_0)}<C$ and the spectrum of $\mathcal{L}_0$ is bounded. Now let $z\in \rho(\mathcal{L}_0)$ with $|z|>C$, where $\rho(\cdot)$ denotes the resolvent set, then we may write the resolvent $R(z;\mathcal{L}_0)$ of the operator $\mathcal{L}_0$ as $$\begin{aligned}
R(z;\mathcal{L}_0) = -z^{-1}\sum_{k = 0}^\infty z^{-k} \mathcal{L}_0^k.\end{aligned}$$ We now show that $R$ is compact. First consider the sequence $(R^N)_{N\geq 1}$ defined by $$\begin{aligned}
R^N(z;\mathcal{L}_0) := -z^{-1}\sum_{k = 0}^N z^{-k} \mathcal{L}_0^k,\end{aligned}$$ then let $(\phi_j)_{j\geq 1}$ be a sequence in $C^2(U)$. We have that $(\phi_j)_{j\geq 1}$ is a bounded sequence in $C^2(U)$ and $$\begin{aligned}
\|R^N(z;\mathcal{L}_0)[\phi_j]\|_{L^2(U,\varrho_0^{-1})} &\leq |z|^{-1}\sum_{k = 0}^N |z|^{-k} \|\mathcal{L}_0^k[\phi_j]\|_{L^2(U,\varrho_0^{-1})}\nonumber\\
&\leq |z|^{-1}\sum_{k = 0}^N |z|^{-k} C^{K}.\end{aligned}$$ Hence, as long as $|z|>C$ then, $\|R^N(z;\mathcal{L}_0)[\phi_j]\|_{L^2(U,\varrho_0^{-1})}$ converges for all $N$ and $\text{Im}\left(R^N\right)$ is relatively compact in $L^2(U,\varrho^{-1}_0)$. It is then a standard result that the limit of a sequence of compact operators is compact, hence $R$ is compact.

Thus we have a complete set of orthonormal basis functions $\{v^{(0)}_k \}_{k=0}^{\infty}$ and corresponding eigenvalues $\{\gamma^{(0)}_n\}_{n\geq 1}$. Note that $v^{(0)}_0 =  \varrho_0$ and $\gamma_0^{(0)} = 0$. At the next order of $\kappa_2$ we obtain $$\begin{aligned}
\label{eq:lin_stab_pert_equilib_order_eps}
\mathcal{L}_0\varrho_1+f(\varrho_0) = 0,\end{aligned}$$ where $$\begin{aligned}
f(\varrho_0):=-\nabla_{\vec{r}} \cdot (\bm{D}\varrho_0\,\nabla_{\vec{r}}V_2\star \varrho_0),\end{aligned}$$ subject to $$\begin{aligned}
\bm{D}(\nabla_{\vec{r}}\varrho_1+\varrho_1\,\nabla_{\vec{r}}V_1 +\varrho_0\,\nabla_{\vec{r}} V_2\star \varrho)\cdot\vec{n}=0 & \text{ on } \partial_U.\end{aligned}$$ The solvability condition for then becomes $$\begin{aligned}
\label{eq:solvability_condition_V2_star_rho0_dot_n}
0=\langle f(\varrho_0),\,v_0^{(0)}\rangle_{L^2(U,\varrho_0^{-1})} & =\int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}} \cdot (\varrho_0\,\bm{D}\nabla_{\vec{r}}V_2\star \varrho_0) = \int_{\partial U}\mathrm{d}S\,\vec{n}\cdot\varrho_0\,\bm{D}\nabla_{\vec{r}}V_2\star \varrho_0.\end{aligned}$$ If the solvability condition is satisfied then, by the Fredholm alternative, there exists a solution to .

We may then write $\varrho_1$ in an eigenfunction expansion $$\begin{aligned}
\varrho_1(\vec{r}) = \sum_{j=0}^{\infty}\alpha_j v_{j}^{(0)} 
\quad \text{ where } \quad 
\alpha_j = -\frac{1}{\gamma_j\|v_j^{(0)}\|^2_{L^2_{\varrho_0^{-1}}}}\langle f(\varrho_0),\,v_j^{(0)} \rangle_{L^2_{\varrho_0^{-1}}}.\end{aligned}$$ This yields that $$\begin{aligned}
\varrho(\vec{r}) = \frac{e^{-V_1(\vec{r})}}{Z_0} + \kappa_2 \sum_{j=0}^{\infty} \frac{\langle  \nabla_{\vec{r}} \cdot (\bm{D} \varrho_0 \,\nabla_{\vec{r}} V_2 \star \varrho_0),\,v_j^{(0)} \rangle_{L^2_{\varrho_0^{-1}}} v_{j}^{(0)} (\vec{r})}{\gamma_j\|v_j^{(0)}\|^2_{L^2_{\varrho_0^{-1}}}} + O(\kappa_2^2).\end{aligned}$$\
We now consider a linear stability analysis of the equilibrium density solving .\

Linear Stability Analysis.
--------------------------

We first investigate the spectrum of the linearised operator $\mathcal{L}_1$ in terms of the eigenspace of its local part. We determine a scheme for computing the eigenvalues of $\mathcal{L}_1$ explicitly. Writing $\varrho = \varrho +\epsilon\,\omega + O(\epsilon^2)$ where $\epsilon \ll 1$ is an arbitrary parameter and not equal to $\kappa_2$, we obtain\

#### $O(\epsilon^{0})$:

$$\begin{aligned}
\mathcal{L}\varrho = 0
 \end{aligned}$$

where we have set $\varrho = \varrho_\infty$ (the unique stationary state) to ease notation and $$\begin{aligned}
 \mathcal{L}\varrho = \nabla\cdot (\bm{D}\,\nabla \varrho) +\kappa_1\nabla\cdot(\bm{D}\,\varrho\nabla V_1) +\kappa_2\nabla\cdot (\varrho\bm{D}\,\nabla V_2\star \varrho).
 \end{aligned}$$

#### $O(\epsilon^{1})$:

$$\begin{aligned}
\label{eq:dynamical_eqn_for_omega}
\dot{\omega} = \mathcal{L}_1w\end{aligned}$$

where $$\begin{aligned}
\label{eq:def_of_linearised_L}
\mathcal{L}_1\omega:=\nabla\cdot (\bm{D}\,\nabla \omega) +\kappa_1\nabla\cdot(\bm{D}\,\omega\nabla V_1) +\kappa_2\nabla\cdot (\varrho\bm{D}\,\nabla V_2\star \omega) + \kappa_2\nabla\cdot (\omega\bm{D}\,\nabla V_2\star \varrho).\end{aligned}$$ We remark that the operator $\mathcal{L}_1$ is different to the one found in the linear stability analysis of [@greg_mckean_vlasov_torus Sec 3.3] due to the difference in boundary conditions.

Perturbations must be mean zero, that is $\int_U\mathrm{d}\vec{r}\,\omega=0$, which may be determined by observing that $$\begin{aligned}
1 = \int\mathrm{d}\vec{r}\,\varrho + \epsilon \int\mathrm{d}\vec{r}\,w + O(\epsilon^2).\end{aligned}$$ Equally, all higher order perturbations must have mean zero. Physically speaking this is a compatibility condition with the no-flux boundary condition in to ensure that perturbations do not change the mass of the system.

Additionally by linearising the self-consistency equation we find that mean zero perturbations $w$ satisfy the integral equation $$\begin{aligned}
\label{eq:int_eqn_for_omega}
w = -\varrho_\infty \kappa_2 V_2\star w.\end{aligned}$$ We linearise the nonlinear boundary condition to find that $$\begin{aligned}
\Pi_1[\omega]\cdot \vec{n}|_{\partial_U}:= 0 \end{aligned}$$ where $$\begin{aligned}
\label{bc:linearised_bc_for_omega}
\Pi_1[\omega] =  \bm{D}\,(\nabla_{\vec{r}}\omega + 
\omega\nabla_{\vec{r}}(\kappa_1 V_1(\vec{r},t)+ \kappa_2 \,V_2\star \varrho) + \kappa_2 \varrho\,\nabla_{\vec{r}}V_2\star \omega).\end{aligned}$$ We note that if any such $\omega$ exist for , then trivially holds, and equation is underdetermined. In order to properly determine $\omega$ we let $\omega \in L^2_{c}(U,\varrho^{-1})$ where $$\begin{aligned}
 L^2_{c}(U,\varrho^{-1}) := \Big\{ u \in  L^2(U,\varrho^{-1}) \, : \, \nabla_{\vec{r}}(\varrho^{-1}u)\cdot \vec{n}|_{\partial U} = 0\Big\}.\end{aligned}$$ The choice $\omega\in  L^2_{c}(U,\varrho^{-1})$ preserves the boundary condition $\Pi_1[\varrho]\cdot \vec{n} |_{\partial U} =0$ and we will show in Lemma \[lem:self\_adjoint\_A\] that it is the most general restriction to ensure that the local part of $\mathcal{L}_1$ is self-adjoint in $L^2(U,\varrho^{-1})$. With this we write

$$\begin{aligned}
\mathcal{L}_{1} &=  \mathcal{A}_{   \kappa_2} + {   \kappa_2}\mathcal{B},\label{eq:def_of_L1_op}\\
\mathcal{A}_{\kappa_2}w &: =\nabla_{\vec{r}}\cdot[ \bm{D} (\nabla_{\vec{r}}w+ w \nabla_{\vec{r}}\varphi_{\kappa_2} )] ,\label{eq:def_of_A_op}\\
\mathcal{B}w&: = \nabla_{\vec{r}}\cdot\left(\varrho\bm{D}\,\nabla_{\vec{r}}V_2\star w\right),\label{eq:def_of_B_op}\\
\varphi_{\kappa_2} &:= \kappa_1 V_1+{\kappa_2} V_2\star \rho.\label{eq:def_of_varphi_k2}\end{aligned}$$

Here, $\mathcal{A}_{\kappa_2}$ and $\mathcal{B}$ are the local and nonlocal parts of $\mathcal{L}_{1}$, respectively. Note however that $\mathcal{A}_{\kappa_2}\neq \mathcal{L}_0$ by definition since $\kappa_2$ is no longer small. All operators $\mathcal{A}_{\kappa_2}$, $\mathcal{B}$, $\mathcal{L}_{1}$ are maps $H^2(U, \varrho_\infty^{-1})\to L^2(U)$. We now show that $A_{\kappa_2}$ is a self-adjoint operator in the space $L^2_{c}(U,\varrho^{-1})$.

[0.485]{} ![Plots of a) The eigenfunctions of $A_{\kappa_2}$ in $L^2([-1,1],\varrho_\infty^{-1})$ as computed with pseudospectral methods for $\kappa_2 = 1$ and $N=100$ spectral points, b) the inner product between pairs of eigenfunctions showing orthogonality of the $\{v_k^{(\kappa_2)}\}$, c) the eigenfunction expansion of $V_2(r) =  1/2(-\tanh((r-1/2)/.05)+\tanh((r+1/2)/.05))$ and d) the absolute error between the expansion $V_{2e}$ and $V_{e}$. The $L^2$ error between $V_2$ and its expansion in eigenfunctions $V_{2e}$ was found to be 5.761e-9.[]{data-label="fig:eigenfunction_figs"}](eigenfuns_k2_1.pdf "fig:"){width="\textwidth"}

 

[0.485]{} ![Plots of a) The eigenfunctions of $A_{\kappa_2}$ in $L^2([-1,1],\varrho_\infty^{-1})$ as computed with pseudospectral methods for $\kappa_2 = 1$ and $N=100$ spectral points, b) the inner product between pairs of eigenfunctions showing orthogonality of the $\{v_k^{(\kappa_2)}\}$, c) the eigenfunction expansion of $V_2(r) =  1/2(-\tanh((r-1/2)/.05)+\tanh((r+1/2)/.05))$ and d) the absolute error between the expansion $V_{2e}$ and $V_{e}$. The $L^2$ error between $V_2$ and its expansion in eigenfunctions $V_{2e}$ was found to be 5.761e-9.[]{data-label="fig:eigenfunction_figs"}](inner_prod_eigenfunctions_on_interval.pdf "fig:"){width="\textwidth"}

 

[0.485]{} ![Plots of a) The eigenfunctions of $A_{\kappa_2}$ in $L^2([-1,1],\varrho_\infty^{-1})$ as computed with pseudospectral methods for $\kappa_2 = 1$ and $N=100$ spectral points, b) the inner product between pairs of eigenfunctions showing orthogonality of the $\{v_k^{(\kappa_2)}\}$, c) the eigenfunction expansion of $V_2(r) =  1/2(-\tanh((r-1/2)/.05)+\tanh((r+1/2)/.05))$ and d) the absolute error between the expansion $V_{2e}$ and $V_{e}$. The $L^2$ error between $V_2$ and its expansion in eigenfunctions $V_{2e}$ was found to be 5.761e-9.[]{data-label="fig:eigenfunction_figs"}](V2_expanded.pdf "fig:"){width="\textwidth"}

 

[0.485]{} ![Plots of a) The eigenfunctions of $A_{\kappa_2}$ in $L^2([-1,1],\varrho_\infty^{-1})$ as computed with pseudospectral methods for $\kappa_2 = 1$ and $N=100$ spectral points, b) the inner product between pairs of eigenfunctions showing orthogonality of the $\{v_k^{(\kappa_2)}\}$, c) the eigenfunction expansion of $V_2(r) =  1/2(-\tanh((r-1/2)/.05)+\tanh((r+1/2)/.05))$ and d) the absolute error between the expansion $V_{2e}$ and $V_{e}$. The $L^2$ error between $V_2$ and its expansion in eigenfunctions $V_{2e}$ was found to be 5.761e-9.[]{data-label="fig:eigenfunction_figs"}](error_V2_and_expanded.pdf "fig:"){width="\textwidth"}

\

\[lem:self\_adjoint\_A\] $\mathcal{A}_{\kappa_2}$ is self-adjoint in $L^2_{c}(U,\varrho^{-1})$.

First note that from and we have that $\nabla_{\vec{r}} \varphi_{\kappa_2} = \rho \nabla_{\vec{r}} \rho^{-1}$ and so $
    \mathcal{A}_{\kappa_2}w = \nabla_{\vec{r}} \cdot [ \rho \bm{D} \nabla_{\vec{r}} ( \rho^{-1} w )].
$ Let $u\in L^2_c(U,\varrho^{-1})$ then $$\begin{aligned}
\langle u, \, \mathcal{A}_{\kappa_2} w\rangle_{L^2(U, \varrho^{-1})} 
&= \int_U\mathrm{d}\vec{r}\, \varrho^{-1} u \mathcal{A}_{\kappa_2} w \nonumber\\
& = \int_U\mathrm{d}\vec{r}\, \varrho^{-1} u \nabla_{\vec{r}} \cdot [ \rho \bm{D} \nabla_{\vec{r}} ( \rho^{-1} w )] \nonumber \\
&=\int_{\partial U}\mathrm{d}S\vec{n}\cdot u\bm{D}\nabla_{\vec{r}}(\varrho^{-1}\omega)-\int_U \mathrm{d}\vec{r}\,  \nabla_{\vec{r}}\left[ \varrho^{-1}u \right]\cdot\left[\varrho \bm{D}\nabla_{\vec{r}}\left(\varrho^{-1}w\right)  \right] \nonumber\\
&= -\int_{\partial U}\mathrm{d}S\vec{n}\cdot w\bm{D}\nabla_{\vec{r}}(\varrho^{-1}u) + \int_U \mathrm{d}\vec{r}\, \nabla_{\vec{r}} \cdot [ \rho \bm{D} \nabla_{\vec{r}} ( \rho^{-1} u )]   \varrho^{-1}w\nonumber\\
&= \langle A_{\kappa_2} u, \,  w\rangle_{L^2(\mathbb{R}, \varrho^{-1})} \nonumber\end{aligned}$$ where we have integrated by parts twice and used that $\bm{D}$ is symmetric and the fact that $u, w\in L^2_c(U,\varrho^{-1})$ to eliminate the boundary terms.

We have established that $\mathcal{A}_{\kappa_2}$ is self-adjoint in $L^2_c(U,\varrho^{-1})$. Additionally we observe that $\mathcal{A}_{\kappa_2}$ has a compact resolvent in $L^2_c(U,\varrho^{-1})$ by a similar result to Lemma \[lem:L0\_has\_compact\_resolvent\]. The spectral theorem therefore provides a complete basis of orthonormal eigenfunctions $v_k^{({\kappa_2})}$ spanning $L^2_c(U,\varrho^{-1})$ with corresponding eigenvalues $\gamma_k^{({\kappa_2})}$ such that $$\begin{aligned}
\label{eq:A_k2_eigenvalue_problem}
A_{\kappa_2} v_{k}^{({\kappa_2})} = \gamma_k^{({\kappa_2})}v_{k}^{(\kappa_2)}.\end{aligned}$$

We note that the operator $\mathcal{B}$ as defined in is, in general, not self-adjoint. From now on we assume that the set of eigenfunctions $\{v_{k}^{({\kappa_2})}\}_{k=1}^\infty$ are normalised to form an orthonormal basis.

[0.485]{} ![Moving eigenvalues for various differential nonlocal operators.[]{data-label="fig:lambda_paths"}](davidson_dodds_example_1.pdf "fig:"){width="\textwidth"}

 

[0.485]{} ![Moving eigenvalues for various differential nonlocal operators.[]{data-label="fig:lambda_paths"}](davidson_dodds_example_2.pdf "fig:"){width="\textwidth"}

 

[0.485]{} ![Moving eigenvalues for various differential nonlocal operators.[]{data-label="fig:lambda_paths"}](tanh_pot_epsilon_lambda.pdf "fig:"){width="\textwidth"}

 

[0.485]{} ![Moving eigenvalues for various differential nonlocal operators.[]{data-label="fig:lambda_paths"}](gaussian_pot_epsilon_lambda.pdf "fig:"){width="\textwidth"}

The stability of the equilibrium density will depend on the spectrum of the operator $\mathcal{L}_1$ so that perturbations evolving according to either grow or decay. We now study the spectrum of $\mathcal{L}_1$. We fix $\kappa_1$ and consider $\kappa_2$, not necessarily small, as a perturbation parameter from the differential part of $\mathcal{L}_1$. The following theorem establishes the parametrisation of the eigenvalues $\lambda$ by $\kappa_2$.\

\[thm:epsilon\_as \_a\_fn\_of\_lambda\]

Suppose that $\lambda \neq \gamma_k^{(\kappa_2)}$ for all $k \in \mathbb{N}$. If the solution $\kappa_2^\star(\lambda)$ of the equation $\lambda = \lambda_{k^\star}(\kappa_2^\star) $ exists, then it is given by $$\begin{aligned}
\label{eq:kappa2_parametrisation}
\kappa_2^\star(\lambda) = \left(\sum_{i=0}^\infty\tfrac{\theta_i^{(\kappa_2)}\gamma
_i^{(\kappa_2)}\beta_i^{(\kappa_2)}}{\lambda-\gamma_i^{(\kappa_2)}}\right)^{-1},\end{aligned}$$ where $\theta_j^{(\kappa_2)}$ and $\beta_j^{(\kappa_2)}$ are given by $$\begin{aligned}
\label{eq:expansion_coeffs_of_V2}
\theta_k^{(\kappa_2)}\beta_l^{(\kappa_2)} = \int_U\mathrm{d}\vec{r}\,v_l^{(\kappa_2)}(\vec{r})V_2\star v_{k}^{(\kappa_2)}(\vec{r}).\end{aligned}$$

Let $$\begin{aligned}
V_2(\vec{r}-\vec{r}') &= \varrho^{-1}(\vec{r})  \varrho^{-1}(\vec{r'}) \sum_{j,k=0}^\infty \beta_j^{(\kappa_2)} v_j^{(\kappa_2)}(\vec{r})\theta_k^{(\kappa_2)} v_k^{(\kappa_2)}(\vec{r}'),\nonumber\\
w(\vec{r}) &= \sum_{i=0}^\infty \alpha_iv_i^{(\kappa_2)}(\vec{r}).\end{aligned}$$ Inserting these expressions into the eigenvalue problem $\mathcal{L}_1 w = \lambda w$ we find $$\begin{aligned}
\sum_{i=1}^\infty \alpha_i(\gamma^{(\kappa_2)}_i-\lambda)v_i^{(\kappa_2)}(\vec{r}) + \kappa_2\nabla_{\vec{r}}\cdot\left(\varrho\bm{D} \nabla_{\vec{r}} V_2\star w\right)  =0.\end{aligned}$$ Multiplying this equation by $v^{(\kappa_2)}_n$ and integrating against the weight function $\varrho^{-1}$ we obtain $$\begin{aligned}
0 &= \alpha_n(\gamma^{(\kappa_2)}_n-\lambda)+\kappa_2\int_U\mathrm{d}\vec{r}\, \varrho^{-1}(\vec{r}) v^{(\kappa_2)}_n(\vec{r})\nabla_{\vec{r}}\cdot\left(\varrho(\vec{r}) \bm{D} \nabla_{\vec{r}} V_2\star w\right) \nonumber\\
&  =\alpha_n(\gamma^{(\kappa_2)}_n-\lambda)
- \kappa_2 \int_U\mathrm{d}\vec{r}\, \left[\nabla_{\vec{r}}v^{(\kappa_2)}_n(\vec{r}) + v^{(\kappa_2)}_n(\vec{r})\nabla_{\vec{r}}\varphi_{\kappa_2} \right] \cdot \bm{D}\, \nabla_{\vec{r}} V_2\star w,\end{aligned}$$ where there is no boundary term since $\nabla_{\vec{r}} V_2\star w = -\kappa_2^{-1}\nabla_{\vec{r}}(\varrho^{-1}w)$ is zero on the boundary of $U$ because $w\in L^2_c(U,\varrho^{-1})$.

Continuing by integrating by parts we find $$\begin{aligned}
0 &= \alpha_n(\gamma^{(\kappa_2)}_n-\lambda)+\kappa_2\int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}} \cdot [\bm{D} ( \nabla_{\vec{r}}v^{(\kappa_2)}_n(\vec{r}) +  v^{(\kappa_2)}_n(\vec{r})\nabla_{\vec{r}}\varphi_{\kappa_2})] V_2\star w\nonumber\\
&= \alpha_n(\gamma^{(\kappa_2)}_n-\lambda)-\kappa_2\int_U\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\cdot(\varrho\bm{D}\nabla_{\vec{r}}(\varrho^{-1}v^{(\kappa_2)}_n(\vec{r}))) V_2\star w\nonumber\\
&= \alpha_n(\gamma^{(\kappa_2)}_n-\lambda)+\kappa_2\gamma^{(\kappa_2)}_n\int_U\mathrm{d}\vec{r}\, v^{(\kappa_2)}_n(\vec{r}) V_2\star w\end{aligned}$$ where we have used $\nabla_{\vec{r}}(\varrho^{-1}v_n^{\kappa_2}) = 0$ on $\partial_U$ to eliminate the boundary term, and in the last line used the fact that $v_n^{\kappa_2}$ is an eigenfunction of $\mathcal{A}_{\kappa_2}$. Inserting the expansion for $V_2$ and using the orthonormality of the $v_i^{(\kappa_2)}$ gives $$\begin{aligned}
\kappa_2 = \frac{(\lambda-\gamma^{(\kappa_2)}_n)\alpha_n}{\gamma_n^{(\kappa_2)}\theta_n^{(\kappa_2)}\sum_{j=0}^\infty \int_U\mathrm{d}\vec{r}'\, \varrho_\infty^{-1}(\vec{r}') \beta_j^{\kappa_2} v_j^{(\kappa_2)}(\vec{r}') w(\vec{r}')}.\end{aligned}$$ This holds for all $V_2$ and, in particular, for all $\theta_n^{(\kappa_2)}\neq 0$ so it be must be the case that $$\begin{aligned}
\frac{(\lambda-\gamma^{(\kappa_2)}_n)\alpha_n}{\gamma_n^{(\kappa_2)}\theta_n^{(\kappa_2)}}=K,\end{aligned}$$ for some constant $K$, independent of $n$. Without loss of generality we can take $K=1$. Hence we have $$\begin{aligned}
w(x) = \sum_{i=0}^\infty \tfrac{\gamma_n^{(\kappa_2)}\theta_n^{(\kappa_2)}}{\lambda-\gamma^{(\kappa_2)}_n}v_n^{(\kappa_2)}(x)\end{aligned}$$ and it follows that $$\begin{aligned}
\label{eq:lem_lambda_eqn_fin_line}
\kappa_2 = \left( \sum_{j=0}^\infty \int_U\mathrm{d}\vec{r}'\, \varrho^{-1}(\vec{r}') \beta_j^{(\kappa_2)} v_j^{(\kappa_2)}(\vec{r}') w(\vec{r}')\right)^{-1} =\left(\sum_{i=0}^\infty\tfrac{\theta_i^{(\kappa_2)}\gamma^{(\kappa_2)}
_i\beta_i^{(\kappa_2)}}{\lambda-\gamma^{(\kappa_2)}_i}\right)^{-1}.\end{aligned}$$ Hence the theorem is proved.

   $k$   $-\gamma_k^{(\epsilon)}\cdot 10^3$   $\nabla_{\vec{r}}(\varrho^{-1}v_k^{(\kappa_2)})\cdot\vec{n}\big|_{x = -1}$   $\nabla_{\vec{r}}(\varrho^{-1}v_k^{(\kappa_2)})\cdot\vec{n}\big|_{x = 1}$
  ----- ------------------------------------ ---------------------------------------------------------------------------- ---------------------------------------------------------------------------
    1            0.042470931917315                                      -0.326629834290770e-10                                                      -0.049430114879555e-10
    2            0.161343622578368                                      0.177791115163473e-11                                                       -0.638172005587933e-11
    3            0.359053918066979                                      -0.532729416136135e-11                                                      -0.705155659306588e-11
    4            0.635777464488092                                      0.040225600628219e-10                                                        -0.10593000196636e-10
    5            0.991543834922971                                      -0.143929312912405e-11                                                      -0.982251087217608e-11
    6            1.426361079097481                                      -0.421040979858844e-11                                                       0.332564751050043e-11
    7            1.940232081076136                                      0.130651045537888e-11                                                        0.966037681231493e-11
    8            2.533158075256826                                      -0.417399448338074e-11                                                      -0.401360089043934e-11
    9            3.205139660109592                                      -0.019828583219805e-11                                                      -0.919929559744713e-11
   10            3.956177153938488                                      -0.687472301308389e-11                                                       0.423840601914807e-11

  : The first 10 eigenvalues $-\gamma_k^{(\kappa_2)}\cdot 10^3$ and boundary condition values of the corresponding eigenvectors $v_{k^{(\kappa_2)}}$ for $\kappa_2 = .05$.[]{data-label="tab:eigenstat_table"}

The expression for $\kappa_2$ allows the paths of the eigenvalues $\lambda_k(\kappa_2)$ to be computed. For practical purposes, it may be sufficient to use a truncation of the series or, if $w(\vec{r})$ can be computed explicitly, the first expression in can be used. As shown in [@kato2013perturbation Section II-5.1], the eigenvalues of $\mathcal{L}_1$ will remain real as long as $$\begin{aligned}
\label{eq:kappa2_bound_spectral_gap}
|\kappa_2|<\frac{\min_{i,j\in\mathbb{N}}|\gamma_i^{(\kappa_2)}-\gamma_{j}^{(\kappa_2)}|}{2\|\mathcal{B}\|}.\end{aligned}$$ We see from that the point of critical stability (if it exists) $\kappa_{2\sharp}$ occurs at $$\begin{aligned}
\label{eq:critical_kappa_2}
\kappa_{2\sharp} = -\left(\sum_{i=0}^\infty\theta_i^{(\kappa_2)}\beta_i^{(\kappa_2)}\right)^{-1}\end{aligned}$$ and is independent of $\gamma_n^{(\kappa_2)}$ (the eigenvalues of the local operator $A_{\kappa_2}$). The critical point of stability will have implicit dependence on $\bm{D}$, $V_1$ and $V_2$ through . As long as $\kappa_2$ remains sufficiently small, Lemma \[thm:epsilon\_as \_a\_fn\_of\_lambda\] provides a nonlinear map to compute $\kappa_2$ parametrised by $\lambda$ therefore permitting the paths of the moving eigenvalues to be calculated. In particular by fixing $\lambda\in\mathbb{R}$ we have the iterative problem $$\begin{aligned}
\label{eq:iteration_for_the_mving_lambdas}
\begin{cases}
\frac{1}{\kappa_2^{l+1}}  =\sum_{i=1}^\infty\tfrac{\theta_i^{(\kappa_2^l)}\gamma^{(\kappa_2^l)}
_i\beta_i^{(\kappa_2^l)}}{\lambda-\gamma^{(\kappa_2^l)}_i},\\
\gamma_{n}^{(\kappa_2^0)} = \gamma_{n}^{(0)}.
\end{cases}\end{aligned}$$ We note that the eigenvalues $\lambda_{k}^{(\kappa_2)}$ are implicitly dependent on the diffusion tensor $\bm{D}$ and confining potential $V_1$.

Figure \[fig:eigenfuns\_k2\_1\] shows typical eigenfunctions $v_k^{(\kappa_2)}$ of the local part of the linearised operator $\mathcal{L}$. Figure \[fig:inner\_prod\_eigenfunctions\_on\_interval\] shows the pairwise $L^2_c(U,\varrho^{-1})$ inner product of the $v_k^{(\kappa_2)}$ demonstrating orthogonality of the basis functions. Figure \[fig:V2\_expanded\] shows the expansion of the two-body function $V_2$ (here a Morse like potential) in terms of the eigenfunctions $v_k$ meanwhile Figure \[fig:error\_V2\_and\_expanded\] shows the error between the expansion and $V_2$. We also demonstrate the accuracy of the collocation scheme in computing eigenvalues and eigenfunctions of $\mathcal{A}_{\kappa_2}$ in Table \[tab:eigenstat\_table\]. In particular, $\mathcal{A}_{\kappa_2}$ is composed of dense first and second order differentiation matrices and the value $\nabla_{\vec{r}}(\varrho^{-1}v_k^{(\kappa_2)})\cdot\vec{n}$ is very small on the boundary using only 100 collocation points.

In Figures \[fig:tanh\_pot\_moving\_eigenvals\], \[fig:gaussian\_pot\_epsilon\_lambda\], we plot various paths $\kappa_2^{\star}(\lambda)$ as solutions to the equation $\lambda = \lambda_{k^\star}(\kappa_2^*)$ for $k$ the wave number by numerically solving for different two-body potentials. We also reproduce figures from Davidson & Dodds [@davidson2006spectral] in Figures \[fig:davidson\_dodds\_example\_1\], \[fig:davidson\_dodds\_example\_2\], verifying our numerical procedure for computing the spectra of similar nonlocal differential operators. Note however that operators in [@davidson2006spectral] do not contain convolution type integral operators, and, with Dirichlet boundary conditions, their spectra differ substantially from those considered here (for example Figures \[fig:tanh\_pot\_moving\_eigenvals\], \[fig:gaussian\_pot\_epsilon\_lambda\]). The intersection through the $\lambda$ axis in each Figure \[fig:davidson\_dodds\_example\_1\]–\[fig:gaussian\_pot\_epsilon\_lambda\] gives the local eigenvalues $\gamma^{(0)}_k$ for the corresponding nonlocal differential operator. Note that it is not necessary for $\gamma^{(0)}_k$ to lie on the moving path for every $k$.

The numerical solution of involves both a truncation of the infinite series and a numerical tolerance for the zeros of the nonlinear function $f(\kappa_2) = \kappa_2-\kappa_2(\lambda_k)$. Note that $\mathcal{L}$ is self-adjoint in $L^2_c(U, \varrho^{-1})$ (with real eigenvalues) only for $\kappa_2 = 0$. The $\lambda$’s are otherwise complex and the curves plotted show when the paths drop to the real plane. When $|\kappa_2|$ is sufficiently large, that is when is violated, the $\lambda$’s have non-zero imaginary part.\
\
We now investigate the spectrum of the linearised operator $\mathcal{L}$ in terms of the eigenspace of its nonlocal part. We determine necessary conditions for bifurcations.

Bifurcation Theory {#sec:bifurcation_theory}
==================

We now provide our first result of the section which relates the stability of equilibrium density to the two-body interaction potential.\

\[thm:stability\_of\_rho\_inf\]

Let $\kappa_2\in(-\infty,\infty)$ and suppose $\varrho$ is a solution to the self-consistency equation . Let $\mathcal{R}$ be given by $$\begin{aligned}
\label{eq:R_int_operator_defn}
 \mathcal{R}w = -\varrho V_2\star w,
 \end{aligned}$$ where $w\in L^2(U,\varrho^{-1})$ is mean zero. If $\mathcal{R}$ is positive definite and $\kappa_2<\beta_1$ where $\beta_1$ is the smallest eigenvalue of $\mathcal{R}^{-1}$, then equilibrium densities formed from repulsive two–body kernels $V_2$ are stable. Conversely if $\mathcal{R}$ is negative definite and $\kappa_2>\beta_1$ where $\beta_1$ is the largest eigenvalue of $\mathcal{R}^{-1}$, then equilibrium densities formed from attractive two–body kernels $V_2$ are stable.

We observe that $\mathcal{L}_1$ is self-adjoint in $L^2(U,\varrho^{-1})$ only when there is no interaction ($\kappa_2 = 0$). We may however expand the eigenfunctions of $\mathcal{L}_1$ in the eigenfunctions of $\mathcal{R}$, $\{u_n\}_{n = 1}^\infty$ which form an orthonormal basis of $L^2(U,\varrho^{-1})$ (see Definition \[def:def\_of\_R\_op\]). We write $w_n = \sum_{i = 1}\alpha_{n_i}u_i$. By the definition of the eigenvalue problem for $\mathcal{L}_1$ $$\begin{aligned}
\mathcal{L}_1w_n = \lambda_nw_n.\end{aligned}$$ Now inserting the expansion in $u_i$’s we obtain $$\begin{aligned}
\lambda_n\sum_{i = 1}\alpha_{n_i}u_i &=\mathcal{L}_1\sum_{i = 1}\alpha_{n_i}u_i \\
&=\left[\mathcal{A}_{\kappa_2}+\kappa_2\mathcal{B}\right]  \sum_{i = 1}\alpha_{n_i}u_i\\
&=\sum_{i = 1}\alpha_{n_i}\left\lbrace\mathcal{A}_{\kappa_2}u_i+\kappa_2\mathcal{B}u_i \right\rbrace \\
&=\sum_{i = 1}\alpha_{n_i}\left\lbrace \nabla_{\vec{r}}\cdot \left(\bm{D}\varrho\left(\nabla_{\vec{r}}(\varrho^{-1}u_i)\right)  \right)-\kappa_2\nabla_{\vec{r}}\cdot \left(\bm{D}\varrho\left(\nabla_{\vec{r}}(\varrho^{-1}(\varrho\mathcal{R}u_i))\right)  \right)\right\rbrace\\
&=\sum_{i = 1}\alpha_{n_i}\left\lbrace\nabla_{\vec{r}}\cdot \left(\bm{D}\varrho\left(\nabla_{\vec{r}}(\varrho^{-1}u_i)\right)  \right)-\kappa_2\nabla_{\vec{r}}\cdot \left(\bm{D}\varrho\left(\nabla_{\vec{r}}(\varrho^{-1}(\varrho\mathcal{R}u_i))\right)  \right)\right\rbrace \\
&=\sum_{i = 1}\alpha_{n_i}\left\lbrace 1-\frac{\kappa_2}{\beta_i}\right\rbrace \nabla_{\vec{r}}\cdot \left(\bm{D}\varrho\left(\nabla_{\vec{r}}(\varrho^{-1}u_i)\right)\right),\end{aligned}$$ where we have used the definitions , and and that each $u_i$ is an eigenfunction of $\mathcal{R}$. Now by multiplying my $\varrho^{-1}u_j$ and integrating we obtain $$\begin{aligned}
\lambda_n\alpha_{n_j}\|u_j\|_{L^2(U,\varrho^{-1})}^2 = \sum_{i = 1}\alpha_{n_i} \left\lbrace 1-\frac{\kappa_2}{\beta_i}\right\rbrace\int_{U}\mathrm{d}\vec{r}\,u_j \nabla_{\vec{r}}\cdot \left(\bm{D}\varrho\left(\nabla_{\vec{r}}(\varrho^{-1}u_i)\right)\right)\end{aligned}$$

Now by integrating by parts, using Gauss’s theorem and the condition that $\nabla_{\vec{r}}(\varrho^{-1}u_i)$ is zero on the boundary of $U$, we obtain $$\begin{aligned}
\lambda_n =\alpha_{n_j}^{-1} \sum_{i = 1}\alpha_{n_i}\left(\frac{\kappa_2}{\beta_i}-1\right)\int_U\mathrm{d}\vec{r}\, \Big|\varrho^{1/2}\bm{D}^{1/2}\nabla_{\vec{r}}\left(\varrho^{-1}u_i\right)\Big|^2  \end{aligned}$$ for every $j = 1,\cdots $. Hence, a bifurcation from the equilibrium density $\varrho$ may occur when $\kappa_2$ coincides with $\beta_j$, for some $j = 1,\cdots $ and perturbations $w_n$ are linear combinations of $u_j$. To ensure $\varrho$ is stable one must have, for every $j\in \mathbb{N}$ $$\begin{aligned}
\begin{cases}
\kappa_2<\beta_j \quad \text{ if } \mathcal{R} \text{ is positive definite,}\\
\beta_j<\kappa_2 \quad \text{ if } \mathcal{R} \text{ is negative definite.}
\end{cases}\end{aligned}$$ Now by the spectral theorem, the $\{\beta_n^{-1}\}_{n\geq 1}$ are discrete, countable and may be ordered such that $|\beta_n^{-1}|\to 0$. Therefore to ensure the stability of $\varrho$ we require $\kappa_2<\beta_1$ if $\mathcal{R}$ is positive definite and $\beta_1<\kappa_2$ if $\mathcal{R}$ is negative definite. This completes the proof of the theorem.

We now relate theorem \[thm:stability\_of\_rho\_inf\] to the H-stability result in [@greg_mckean_vlasov_torus].\
\

\[rem:estimate\_betans\] We remark on the consistency with the H-stability condition of [@greg_mckean_vlasov_torus] with periodic boundary conditions, the equilibrium density may bifurcate if the interaction kernel has a negative Fourier mode. In the present work, the distribution of the eigenvalues of the operator $\mathcal{R}$ determines whether the equilibrium density is stable with respect to $\{u_j\}_{j = 1}^{\infty}$. In particular, if $\mathcal{R}$ has a negative eigenvalue then equilibrium densities formed from repulsive $V_2$ may become unstable.

We may obtain an estimate for the eigenvalues $\beta_n^{-1}$ in terms of $V_2$ and $\varrho$ in the following way, by the eigenvalue problem we have $$\begin{aligned}
|\beta_n^{-1}| &= |\beta_n^{-1}|\langle u_n,u_n\rangle_{L^{2}(U,\varrho^{-1})} = |-\langle u_n,V_2\star u_n\rangle_{L^{2}(U)}| \nonumber\\
&\leq \|V_2\|_{L^\infty(U)} \|u_n\|^2_{L^1(U)}= \|V_2\|_{L^\infty(U)}\|\varrho^{1/2}(\varrho^{-1/2})u_n\|_{L^1}^2\nonumber\\
&\leq \|V_2\|_{L^\infty(U)}\|\varrho\|_{L^1(U)}^2\|(\varrho^{-1/2})u_n\|_{L^2(U)}^2= \|V_2\|_{L^\infty(U)},\end{aligned}$$ where we have used the Cauchy-Schwarz inequality and the fact that the $\{u_n\}_{n = 1}^\infty$ are orthonormal in $L^2(U,\varrho^{-1})$. From this we obtain the lower bound $\|V_2\|_{L^\infty(U)}^{-1}\leq |\beta_n|$, this lower bound shows that the bifurcation point coincides with the boundary of the interval in which free energy $\mathcal{F}$ is convex (c.f. Proposition \[prop:F\_is\_strictly\_convex\]).

\[thm:necessary\_bifurcation\_conditions\] Let $\{\beta_{n}^{-1}\}_{n=1}^{\infty}$ be the ordered eigenvalues of $\mathcal{R}$. If $|\kappa_2|\geq |\beta_1|$ then $(\beta_1,w_1)$ is a bifurcation point of where $w_1$ is the eigenfunction of $\mathcal{R}$ associated to $\beta_1^{-1}$ and there exists $0<\varrho_\ast\neq \varrho_{\infty}$ solving .

Let $\varrho_{\kappa_2}$ denote the solution to for a given $\kappa_2$ which is known to exist by Theorem \[thm:exis\_fix\_point\]. Since $\varrho_{\kappa_2}$ is continuous in $\kappa_2$ and $\mathcal{F}[\varrho]$ is continuous in $\varrho$, then $\mathcal{F}$ is continuous in $\kappa_2$. By Lemma \[lem:minimisers\_always\_exist\] we know that a minimiser of $\mathcal{F}$ exists for each $\kappa_2$ and by Lemma \[lem:minimisers\_are\_positive\] the minimiser is strictly positive. Given $|\kappa_2|\geq \|V_2\|_{L^\infty(U)}^{-1}$ then by Proposition \[prop:F\_is\_strictly\_convex\], $\mathcal{F}$ is no longer convex and $\varrho_{\kappa_2}$ is either an inflection point or a local maximum of $\mathcal{F}$. Hence $\varrho_{\kappa_2}$ is unstable and by Lemma \[lem:minimisers\_always\_exist\] there exists $\varrho_{\ast}$ such that $\mathcal{F}[\varrho_\ast]<\mathcal{F}[\varrho_{\kappa_2}]$. Additionally by the self-adjointness and compactness of $\mathcal{R}$, one has that $\beta_n^{-1}\to 0$ as $n\to \infty$ and hence $\beta_n\to \infty$ as $n\to \infty$ and $\beta_1$ is the smallest of the $\{\beta_{n}\}_{n=1}^{\infty}$.

If $\mathcal{R}$ is positive definite, there are no negative $\beta_n$ and the only solution to is $u_n\equiv 0$ and $\varrho_{\kappa_2}$ will be stable for all $\kappa_2<\beta_1$. Similarly, if $\mathcal{R}$ is negative definite, there are no positive $\beta_n$ and the only solution to is $u_n\equiv 0$ and $\varrho_{\kappa_2}$ will be stable for all $\kappa_2>\beta_1$. If $\mathcal{R}$ is indefinite, by Remark \[rem:estimate\_betans\], for $|\beta_n|<\|V_2\|^{-1}_{L^\infty(U)}$ there are no solutions (other than $w_n\equiv 0$) to $\mathcal{R}[w_n] = \beta_n^{-1}w_n$, and once again for $|\kappa_2|<|\beta_1|$, $\varrho_{\kappa_2} = \varrho_{\infty}$ is stable. For $\kappa_2\geq \|V_2\|^{-1}_{L^\infty(U)}$ there are infinitely many non-trivial solutions to $\mathcal{R}[w_n] = \beta_n^{-1}w_n$ and $\kappa_2 = \beta_1$ is the first.

Hence if $|\kappa_2|\geq|\beta_1|$ then the unique stationary density $\varrho_\infty$ is unstable and by Lemma \[lem:minimisers\_always\_exist\] there must exist $\varrho_{\ast}$ such that $\mathcal{F}[\varrho_\ast]<\mathcal{F}[\varrho_{\kappa_2}]$.

We define the $\mathcal{W}:L^2(U)\to \mathbb{R}$ transform such that $$\begin{aligned}
\mathcal{W}[f](n) = \int_U\mathrm{d}\vec{r}'\varrho^{-1}_{\beta_n}w_n(\vec{r})f(r)\end{aligned}$$ where $\varrho_{\beta_n}$ solves with $\kappa_2 = \beta_n$. With this we may plot the bifurcation diagram for the stability of the unique equilibrium state $\varrho = \varrho_\infty$, see for example Figure \[fig:bifurcation\_diagram\].

Application To Nonlinear Diffusion Equations {#sec:manufactured_bif}
============================================

In this section we consider sufficient conditions for bifurcations under particular forms of nonlocal operators. We will show that, by use of numerical examples, there may be more than one stationary solution under additional assumptions on the two-body potential by making use of the bifurcation theory developed in Section \[sec:bifurcation\_theory\]. We fix $\kappa_1$ and consider boundary value problems where the nonlocal term is not of convolution type. Let $V_2(\vec{r},\vec{r}')$ be a two-body function and consider $$\begin{aligned}
\begin{cases}\label{eq:bif_stationary_eqn}
\mathcal{P}[\varrho]:=\nabla \cdot \Big[ \bm{D} \Big(\nabla\varrho + \varrho \kappa_1 \nabla V_1+\kappa_2 \varrho \nabla \int_U\mathrm{d}\vec{r}' V_2(\vec{r},\vec{r}') \varrho(\vec{r'}) \Big) \Big]= 0 & \text{ in } U,\\
\Omega [\varrho]\cdot\vec{n}  := \bm{D} \Big(\nabla\varrho + \varrho \kappa_1 \nabla V_1+\kappa_2 \varrho \nabla \int_U\mathrm{d}\vec{r}' V_2(\vec{r},\vec{r}') \varrho(\vec{r'}) \Big) \cdot \vec{n}= 0 & \text{ on } \partial U.
\end{cases}\end{aligned}$$ Solutions of $\mathcal{P}\varrho=0$ with $\Omega [\varrho]\cdot\vec{n} = 0 $ on the boundary are denoted by $\varrho=\varrho_{\kappa_2}$ and satisfy the self-consistency equation $$\begin{aligned}
\label{eq:bifur_self_con_eqn}
\varrho_{\kappa_2} = \frac{e^{-(\kappa_1V_1+\kappa_2\int_U\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\varrho_{\kappa_2}(\vec{r}'))}}{Z}.\end{aligned}$$

The linear stability of the steady state may be studied implicitly by examining the properties of the linearised self-consistency map. By linearising equation , by writing $\varrho_{\kappa_2} = \phi_0+\epsilon\phi_1$, for some small $\epsilon$, we obtain the original nonlinear problem $$\begin{aligned}
\label{eq:bifur_self_con_eqn_O_0}
\phi_{0} = \frac{e^{-(\kappa_1V_1+\kappa_2\int_U\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\phi_0(\vec{r}'))}}{Z_0} \quad \text{ s.t } \quad \Omega[\phi_{0}]\cdot \vec{n}=0\end{aligned}$$ where $Z_0 = \int_U\mathrm{d}\vec{r}\, e^{-(\kappa_1V_1+\kappa_2\int_U\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\phi_0(\vec{r}'))}$, along with the linearised equation $$\begin{aligned}
\label{eq:rho_pert_series_u1}
\phi_1 = - \kappa_2\,\phi_0\int_U\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\phi_1(\vec{r}') \quad \text{ s.t } \quad \int_U\mathrm{d}\vec{r}\, \phi_1(\vec{r}) = 0.\end{aligned}$$ The integral condition in comes from the fact that higher order perturbations to $\phi_0$ must possess zero mean to preserve the mass in the system.

We define the linear operator $\mathcal{T}$ in $L_1(U)$ by $$\begin{aligned}
\label{eq:def_of_T_operator}
\mathcal{T}\phi(\vec{r}) := \phi_0(\vec{r})\int_U\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\phi(\vec{r}').\end{aligned}$$ We also define the mapping from $\mathcal{G}:\left(L_1(U), \mathbb{R}\right) \to L_1(U)$ by $$\begin{aligned}
\mathcal{G}(v,\kappa) := \phi-f(\phi,\kappa)\end{aligned}$$ where $f(\phi,\kappa) := \tfrac{e^{-(\kappa_1V_1+\kappa\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\phi(\vec{r}'))}}{\int_U\mathrm{d}\vec{r}e^{-(\kappa_1V_1+\kappa\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\phi(\vec{r}'))}}$. To construct the bifurcation diagram, we will use the following result from [@tamura1984asymptotic Tamura (1984)], or [@greg_mckean_vlasov_torus Carrillo et al. 2019], which is a direct consequence of the Crandall-Rabinowitz theorem, see, e.g. [@crandall1971bifurcation].\

\[thm:tamura\_bifurcations\]

Let $V_2(x,y) = V_2(y,x)$. Also let $(\psi_0,\mu_0)$ be a fixed point in $L^1(U)\times \mathbb{R}$ such that:

1.  $\mathcal{G}(\psi_0,\mu_0)=0$,

2.  $\mu_0^{-1}$ is an eigenvalue of $\mathcal{T}$,

3.  $\int_U\mathrm{d}\vec{r}\,V_2(\vec{r},\vec{r}')\psi_0(\vec{r})=0$,

4.  $\dim \{\phi\in L^1(U)\,:\, v = \mu_0 \mathcal{T}\phi\}=1$.

Then $(\psi_0,\mu_0)$ is a bifurcation point of $\mathcal{G}=0$. That is, for any neighbourhood $B$ of $(\psi_0,\mu_0)$ in $L^1(U)\times \mathbb{R}$ there exists $(\psi_1,\mu_1)\in B$ such that $\psi_1\neq \psi_0$ and $\mathcal{G}(\psi_1,\mu_1)=0$.

\[thm:bifurcations\]

The proof relies on checking the conditions of the Crandall-Rabinowitz Theorem and is equivalent to Tamura’s proof [@tamura1984asymptotic].

[0.48]{} ![Stable densities bifurcating from (a). $\psi_0 = \exp\{-\kappa_1V_1(x)\}/Z$ and (b). $\psi_0 = \frac{N}{2\mathrm{L}}$ where $2L$ is the length of the interval, which solve for different two-body functions (a). $V_2(x,y) = xy$ and (b). $V_2(x,y) = -\cos\left(\frac{2\pi(x-y)}{\mathrm{L}}\right)$ . Insets show the shape of perturbation function.[]{data-label="fig:bifurcations"}](xy_bifurcation_perturb_2.pdf "fig:"){width="\textwidth"}

 

[0.48]{} ![Stable densities bifurcating from (a). $\psi_0 = \exp\{-\kappa_1V_1(x)\}/Z$ and (b). $\psi_0 = \frac{N}{2\mathrm{L}}$ where $2L$ is the length of the interval, which solve for different two-body functions (a). $V_2(x,y) = xy$ and (b). $V_2(x,y) = -\cos\left(\frac{2\pi(x-y)}{\mathrm{L}}\right)$ . Insets show the shape of perturbation function.[]{data-label="fig:bifurcations"}](stable_equilibria_V2_cosine.pdf "fig:"){width="\textwidth"}

Note that $\psi_0$ is, by construction, the background density given by $\psi_0 = \tfrac{e^{-V_1(\vec{r})}}{\int_U\mathrm{d}\vec{r}e^{-V_1(\vec{r})}}$. Theorem \[thm:tamura\_bifurcations\] presents sufficient conditions to permit bifurcations from $v_0$ with stationary equations of the form . In particular it will be sufficient that the two-body potential satisfies the normality condition (condition 3. of Theorem \[thm:tamura\_bifurcations\]). Then bifurcations occur at discrete eigenvalues of the nonlocal operator $\mathcal{T}$ as defined in . We remark that these conditions are consistent with Theorem \[thm:necessary\_bifurcation\_conditions\].

Numerical Experiments. {#subsec:numericalexperiments}
----------------------

In this section we compute the branches of solutions that may evolve in the DDFT-like example considered in Section \[sec:manufactured\_bif\] with nonlinear, nonlocal boundary conditions. Given simple interaction kernels we show that symmetry-breaking systems may be constructed quite easily given sufficiently high interaction strength. For the numerical examples presented here, $\varrho$ is a number density and hence $\int_U\mathrm{d}\vec{r}\,\varrho = 0$. We consider numerical solutions to $$\begin{aligned}
\begin{cases}\label{eq:bif_evolve_eqn}
\partial_t \varrho = \nabla \cdot \Big[ \bm{D} \Big(\nabla\varrho + \varrho \kappa_1 \nabla V_1+\kappa_2 \varrho \nabla \int_U\mathrm{d}\vec{r}' V_2(\vec{r},\vec{r}') \varrho(\vec{r'}) \Big) \Big] & \text{ in } U,\\
\Omega [\varrho]\cdot\vec{n}  := \bm{D} \Big(\nabla\varrho + \varrho \kappa_1 \nabla V_1+\kappa_2 \varrho \nabla \int_U\mathrm{d}\vec{r}' V_2(\vec{r},\vec{r}') \varrho(\vec{r'}) \Big) \cdot \vec{n}= 0 & \text{ on } \partial U, \\
\varrho(\vec{r},0) =  \tfrac{e^{-\kappa_1V_1(\vec{r})+\kappa_2\int_U\mathrm{d}\vec{r}'\,V_2(\vec{r},\vec{r}')\varrho(\vec{r}',0)}}{Z} \quad &\text{ at } t=0.
\end{cases}\end{aligned}$$\
The nonlocal terms in , both in the evolution equation and the boundary condition, mean that numerical implementations require efficient and accurate quadrature. We demonstrate the power with which the pseudo-spectral collocation scheme 2DChebClass [@DDFTCode] may compute solutions with such efficiency and accuracy. For a more detailed explanation of pseudospectral methods for DDFT problems, particularly the efficient computation of convolution integrals, see [@nold2017pseudospectral].

Some numerical experiments were performed by solving in 1D with the choice $V_2 = xy$ and $V_1 = \kappa_1 x^2$ on $U = [-1/2,1/2]$. Under this choice of confining and two-body potentials the normality condition (3) of Theorem \[thm:tamura\_bifurcations\] holds. Additionally, for $|\kappa_2|$ sufficiently small, the unique stationary density is $v_0 = e^{-\kappa_1V_1}/Z$. Upon increasing $\kappa_2$ and perturbing with a mean zero function $\eta(x,\theta)$ the stability of $v_0$ breaks and transitions may be observed to non symmetric equilibria. The asymmetry of the equilibria depends on the sign of $\eta$ as seen in Figure \[fig:bifurcations\].

Figure \[fig:bifurcations\] shows long time numerical solutions to the IBVP subject to a mean zero perturbation for different interaction strengths $\kappa_2$, with $V_2$ fixed. In Figure \[fig:bif\_fig\_right\], the symmetric solution $\psi_0 = \exp(-\kappa_1 V1)/Z$ was shown to be unstable as interaction strength $\kappa_2$ was made ever negative. In particular by perturbing with a sinusoidal function with positive or negative sign, the stationary density can be shown to adhere to one boundary, thereby bifurcating from the previously symmetric solution $\psi_0$. The skewness of the density is controlled by the sign of the perturbation function $\eta$ and $\eta\in \text{Span}{\mathcal{T}}$, hence densities which adhere to the left boundary may be obtained by changing the sign of $\eta$. We predict the stable and symmetric branch to bifurcate at the critical interaction energy $\kappa_2 = -2.4$ (to 1 decimal place) which is the negative inverse of smallest eigenvalue of $\psi_0^{-1}\mathcal{T}$ in $U = [-1/2,1/2]$. This is verified in Figure \[fig:bifurcation\_diagram\] and the transition between a stable symmetric density and a stable nonsymmetric one is observed in Figure \[fig:bif\_fig\_right\] for the curves labelled $\kappa_2 = -2$ and $\kappa_2 = -3$. In Figure \[fig:bif\_stable\_equilibria\_V2\_cosine\], we see how the uniform density may become unstable. Here $\psi_0 = N/(2\mathrm{L})$ where $2L$ is the length of the interval. We perturb with eigenvectors of $\mathcal{T}$ at increasing interaction strengths. The critical strength was $\kappa_{2\sharp} = 0.4$ (to 1 decimal place), the negative inverse of smallest eigenvalue of $\psi_0^{-1}\mathcal{T}$. This is verified in Figure \[fig:bifurcation\_diagram\] and the transition between a stable uniform density and a stable multi-modal one is observed in Figure \[fig:bif\_stable\_equilibria\_V2\_cosine\] for the curves labelled $\kappa_2 = 0$ and $\kappa_2 = 0.5$.

Existence & Uniqueness of Weak Solutions to Density with Full HI {#sec:existence_uniqueness_with_partial_HI}
================================================================

In this section we determine the existence and uniqueness of the weak density $\varrho(\vec{r},t)$ solving in the sense . To ease notation we suppress $\bm{A}[\vec{a}]$ as it may be trivially added (see Remark \[rem:final\_assumptions\]). We begin by determining some useful results: first, that $\varrho(\vec{r},t)$ is bounded above in $L^1(U)$ for all time by initial data $\varrho_0$ and second, the $L^1(U)$ norm of $\varrho$ is unity for all time and $\varrho(\vec{r},t)$ is non-negative. We will strengthen the non-negativity to strict positivity of $\varrho(\vec{r},t)$ in Section \[subsec:strict\_pos\_rho\]. The results in this section are analogous to those in [@chazelle2017well], [@greg_mckean_vlasov_torus] with the difference that the boundary conditions we consider are no-flux the diffusion tensor is non-constant.

Useful Results.
---------------

We identify the expansion of the absolute value function.

\[def:def\_of\_chi\_approx\_abs\]

Let $\epsilon>0$ and define the convex $C^2$ approximation of $|\cdot|$ by $$\begin{aligned}
\chi_\epsilon(\psi) = \begin{cases}
|\psi| \quad \text{ for } \quad \psi>\epsilon,\\
-\tfrac{\psi^4}{8 \epsilon^3 }+\tfrac{3 \psi^2}{4 \epsilon}+ \tfrac{3 \epsilon}{8} \quad \text{ for } \quad \psi\leq \epsilon.
\end{cases}\end{aligned}$$

We now present our first result concerning the boundedness of the the $L^1$ norm of $\varrho$ in terms of the initial data $\varrho_0$.\

If $ \varrho\in C^1([0,\infty);C^2(U))$ is a solution of with $\varrho_0\in L^1(U)$ then $\|\varrho(t)\|_{L^1(U)}\leq \|\varrho_0\|_{L^1(U)}$ for all time $t\geq 0$.

Multiplying by $\chi_\epsilon'(\varrho)$, integrating and using the divergence theorem and chain rule, we have $$\begin{aligned}
&\der[]{t}\int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho)+\|\bm{D}^{1/2}\nabla_{\vec{r}}\varrho\,[\chi_\epsilon''(\varrho)]^{1/2}\|_{L^2(U)}^2\nonumber\\
&=-\int\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\varrho\, \chi_\epsilon''(\varrho)\cdot[\varrho\,\bm{D}(\vec{r})\nabla_{\vec{r}}(\kappa _1V_1(\vec{r})+\kappa _2[V_2\star\varrho](\vec{r}) )].\end{aligned}$$ Now by H[ö]{}lder’s inequality and then Young’s inequality $$\begin{aligned}
&\der[]{t}\int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho)+\|\bm{D}^{1/2}\nabla_{\vec{r}}\varrho\,[\chi_\epsilon''(\varrho)]^{1/2}\|_{L^2(U)}^2\nonumber\\
&\leq \|\bm{D}^{1/2}\nabla_{\vec{r}}\varrho [\chi_\epsilon''(\varrho)]^{1/2}\|_{L^2(U)}
\times \| [\chi_\epsilon''(\varrho)]^{1/2}\varrho \,\bm{D}^{1/2}\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2[V_2\star \varrho] )\|_{L^2(U)}\nonumber\\
&\leq \tfrac{1}{2}\|\bm{D}^{1/2}\nabla_{\vec{r}}\varrho [\chi_\epsilon''(\varrho)]^{1/2}\|_{L^2(U)}^2
+\tfrac{1}{2}\| [\chi_\epsilon''(\varrho)]^{1/2}\varrho\, \bm{D}^{1/2}\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2[V_2\star \varrho] )\|_{L^2(U)}^2.\end{aligned}$$ Note there are no boundary terms due to the condition $\Pi[\varrho]\cdot\vec{n}=0$ on $\partial U$. All together this implies the inequality $$\begin{aligned}
&\der[]{t}\int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho)+\tfrac{1}{2}\|\bm{D}^{1/2}\nabla_{\vec{r}}\varrho\,\chi_\epsilon''(\varrho)^{1/2}\|_{L^2(U)}^2\nonumber\\
&\leq \tfrac{1}{2} \| [\chi_\epsilon''(\varrho)]^{1/2}\varrho\, \bm{D}^{1/2}\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2[V_2\star\varrho] )\|_{L^2(U)}^2\nonumber\\
&\leq \tfrac{1}{2}\| \bm{D}^{1/2}\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2[V_2\star\varrho] )\|_{L^\infty}^2
\| [\chi_\epsilon''(\varrho)]^{1/2} \varrho \|_{L^2}^2\nonumber\\
&\leq c_0 \| [\chi_\epsilon''(\varrho)]^{1/2} \varrho \|_{L^2}^2 (1+\|\varrho\|_{L^1(U)}^2)\label{eq:ddt_chi_rho_bound}\end{aligned}$$ for the constant $c_0 = 2\mu_{\max} \max\{|\kappa_1|^2 \|\nabla_{\vec{r}}V_1\|_{L^\infty}^2,|\kappa_2|^2 \|V_2\|_{L^\infty}^2\}$.

It is an elementary calculation to show that $$\begin{aligned}
\varrho^2\chi_\epsilon''(\varrho) = \tfrac{3\varrho^2}{2\epsilon} - \tfrac{3\varrho^4}{2\epsilon^3}\end{aligned}$$ for $\varrho\leq \epsilon$. With this, and the fact that $\chi''(\varrho) = 0$ for $\varrho>\epsilon$, we have $$\begin{aligned}
\| [\chi_\epsilon''(\varrho)]^{1/2} \varrho \|_{L^2}^2 &= \int_U  \mathrm{d}\vec{r}\, \varrho^2 \chi_\epsilon''(\varrho) \mathbb{I}_{\varrho\leq \epsilon}
+ \int_U  \mathrm{d}\vec{r}\, \varrho^2 \chi_\epsilon''(\varrho) \mathbb{I}_{\varrho > \epsilon} \nonumber \\
& = \int_U  \mathrm{d}\vec{r}\, \frac{3 \varrho^2(\epsilon^2 - \varrho^2)}{2 \epsilon^2} \mathbb{I}_{\varrho\leq \epsilon}
\leq \int_U  \mathrm{d}\vec{r}\, \frac{3 \epsilon}{2}\mathbb{I}_{\varrho\leq \epsilon} \leq c_1\epsilon \label{eq:sqrt_chi''_rho_L2}\end{aligned}$$ for some constant $c_1$ dependent on $U$. Applying Gr[ö]{}nwall’s lemma to $\eta(\cdot)$ a non-negative, absolutely continuous function on $[0,T]$ which satisfies for a.e. $t$ $$\begin{aligned}
\eta'(t)\leq \phi(t)\eta(t) + \psi(t)\end{aligned}$$ where $\phi$, $\psi$ non-negative and integrable functions on $[0,T]$ gives $$\begin{aligned}
\label{eq:statement_gronwall}
\eta(t)\leq e^{\int\mathrm{d}s_0^t\,\phi(s)}\eta(t)\Big[ \eta(0)+ \int_0^t\mathrm{d}s\,
\psi(s)\Big].\end{aligned}$$

Observe that $\|\varrho\|_{L^1(U)}\leq \int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho)$. Using this with , and with $\eta(t) = \phi(t) =c_1\epsilon \int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho)$ and $\psi(t) = c_1\epsilon$ we obtain $$\begin{aligned}
\int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho)\leq \left(\int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho_0)+c_1\epsilon \,t\right)\, e^{c_1\epsilon\int_0^t\mathrm{d}s\,\int_U\mathrm{d}\vec{r}\,\chi_\epsilon(\varrho(\vec{r},s))}.\end{aligned}$$ Now since $\varrho$ is assumed to be continuous in time on $[0,\infty)$ the integral in the exponential is finite. Therefore taking $\epsilon \to 0$ one obtains $$\begin{aligned}
\|\varrho\|_{L^1}\leq \|\varrho_0\|_{L^1}\end{aligned}$$ for every $t>0$.

\[cor:L\_1\_varrho\_is\_1\]

If $ \varrho\in C^1([0,\infty);C^2(U))$ is a solution of with $\varrho_0$ a probability density, that is $\varrho_0\geq 0$ and $\int_U \mathrm{d}\vec{r}\,\varrho_0(\vec{r})=1$, then $\|\varrho(t)\|_{L^1(U)}=1$ and $\varrho(t)\geq 0$ in $U$ for all time $t\geq 0$.

The argument is a standard one. Since, due to no-flux boundary conditions, $ \der[]{t} \int \mathrm{d}\vec{r} \, \rho(\vec{r},t) = 0$, we have $$\begin{aligned}
1  = \int_U\mathrm{d}\vec{r}\, \varrho_0(\vec{r}) = \int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t) \leq \|\varrho(t)\|_{L^1(U)} \leq \|\varrho(0)\|_{L^1(U)} = \int_U\mathrm{d}\vec{r}\, \varrho_0(\vec{r}) = 1,\end{aligned}$$ so $\|\varrho(t)\|_{L^1(U)} = 1$. Also observe the two equalities $$\begin{aligned}
&1 = \int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t) = \int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t)\mathbb{I}_{\varrho\geq 0} + \int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t)\mathbb{I}_{\varrho< 0},\nonumber\\
&1 = \int_U\mathrm{d}\vec{r}\, |\varrho(\vec{r},t)| = \int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t)\mathbb{I}_{\varrho\geq 0} - \int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t)\mathbb{I}_{\varrho< 0},\nonumber\end{aligned}$$ where in the second line we have used the definition of the absolute value function. Subtracting these equalities we obtain $$\begin{aligned}
2\int_U\mathrm{d}\vec{r}\, \varrho(\vec{r},t)\mathbb{I}_{\varrho< 0} = 0\end{aligned}$$ which implies $ \varrho(\vec{r},t)\geq 0$ almost everywhere in $U$. Non-negativity of $\varrho$ on all of $U$ follows from continuity.

With these results we may continue to determine the existence and uniqueness of weak densities solving in the sense . The method we use follows [@chazelle2017well] but here we must include calculations for the confining potential $V_1^{\text{eff}}$ (which for ease of notation is written $V_1$ for each $\vec{a}(\vec{r},t)$) and a much wider class of two-body potentials $V_2$ which are not necessarily step functions. To start we introduce , the frozen version of , indexed by $n\in \mathbb{N}$, by substituting $\varrho = u_n$ everywhere except in the convolution term where we substitute $\varrho = u_{n-1}$. Each equation is parametrised by $n$, a linear parabolic PDE for the unknown $u_n$ in terms of the solution $u_{n-1}$ at the previous index, for which we have existence and uniqueness of weak solutions for each $n$. The remainder of the argument is to show $\lim_{n\to\infty} u_{n}$ exists and is a limit point solving the weak problem . In this section we will make references to Appendix \[sec:classical\_paraboliv\_pde\] for results and definitions required for $u_{n}\in H^1(U)$, which differ slightly from the standard arguments found in textbooks for classical linear PDE theory (e.g. [@evans2002partial]).

Energy Estimates.
-----------------

The results of Appendix \[sec:classical\_paraboliv\_pde\] are that the initial boundary value problem $$\begin{aligned}
\label{eq:well_posed_un_ibvp}
\begin{cases}
\quad\partial_t u_n-\nabla_{\vec{r}}\cdot[\bm{D}\nabla_{\vec{r}}u_n]=\nabla_{\vec{r}}\cdot[u_n\,\bm{D}\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2V_2\star u_{n-1})],\\
\qquad \qquad \qquad \qquad \Xi [u_n]\cdot\vec{n} = 0 \quad \text{ on } \partial U \times [0, T],\\
\qquad\Xi[u_n]:= \bm{D}\,(\nabla_{\vec{r}}u_n+ 
u_n\, \nabla_{\vec{r}}(\kappa_1 V_1(\vec{r},t)+ \kappa_2 V_2\star u_{n-1})),\\
\qquad \qquad \qquad \qquad \quad  u_n = \varrho_0 \quad \text{ on } U\times \left\lbrace t=0\right\rbrace
\end{cases}\end{aligned}$$ is well posed, and there exists weak solutions $u_{n}$ for each $n\in\mathbb{N}$ in the sense . All that remains is to take the limit $n\to \infty$ to recover the original Smoluchowski equation . We start by deriving our first estimate on energy of $u_n$. To ease notation we derive all results with the time dependence on $\bm{D}$ suppressed since time may be trivially added to the exposition. Additionally, for a stationary density one has $$\begin{aligned}
\lim_{t\to\infty}\bm{D}(\vec{r},t) = \left(\bm{1} + \int \mathrm{d}\vec{r}' g(\vec{r},\vec{r}')\bm{Z}_{1}(\vec{r},\vec{r}')\varrho(\vec{r}')\right)^{-1}  \end{aligned}$$ which is a positive definite tensor and hence diagonalisable, and may be bounded by its smallest and largest eigenvalues which are positive and finite for $t\to\infty$. Hence energy estimates remain valid for $0<t\leq T$ when provided in terms of $\mu_{\min}$ and $\mu_{\max}$, both eigenvalues which depend on time but always remain positive and finite. It will be seen that a natural dual space to $H^1(U)$ is provided by the no-flux condition. In particular we denote by $H^{-1}(U)$ the dual space of $H^{1}(U)$, this is due to the divergence theorem and the boundary condition $\Xi[u_n]\cdot\vec{n}=0$ on $\partial U \times [0,T]$, there is no boundary term, and the normal characterisation of $H^{-1}=(H^1_0)^\ast$ carries over to $H^{1}(U)$.

We now obtain uniform estimates on $u_n$ in terms of the initial data $\varrho_0$ in all the required energy norms. The detailed calculations follow [@chazelle2017well] but take into account the confining potential and non-constant diffusion tensor $\boldsymbol{D}$. The explicit calculations can be found in RDMW’s PhD thesis [@rdmwthesisddft]. The first estimate is in $L^{\infty}([0,T];L^2(U))$ and $L^{2}([0,T];H^1(U))$ norms.\

\[prop:bound\_rhon\]

Let $T>0$ and suppose $\{u_n\}_{n\geq 1}$ satisfies with $\varrho_0\in C^{\infty}(U)$ a probability density. Then there exists a constant $C(T)$, dependent on time and $\mu_{\max}$, such that $$\begin{aligned}
\label{eq:L2_H_1_bound_for_rhon}
\|u_n\|_{L^{\infty}([0,T];L^2(U))}+\|u_n\|_{L^{2}([0,T];H^1(U))}\leq C(T,\mu_{\max}) \|\varrho_{0}\|_{L^2(U)}.\end{aligned}$$

The second estimate is for $L^{\infty}([0,T]; H^1(U))$ and $L^2([0,T];L^2(U))$ norms.\

\[prop:un\_H1\_rho0\_bound\]

Let $T>0$ and suppose $\{u_n\}_{n\geq 1}$ satisfies with $\varrho_0\in C^{\infty}(U)$ a probability density. Then there exists some constant dependent on time $C(T)$ such that $$\begin{aligned}
&\|u_n\|_{L^{\infty}([0,T]; H^1(U))}+\|\nabla_{\vec{r}}\cdot[\bm{D}\,\nabla_{\vec{r}}u_n]\|_{L^2([0,T];L^2(U))}^2\nonumber\\
&\leq C(T)(\|\varrho_0\|_{H^1(U)}^2+(1+\|\varrho_0\|_{L^2(U)})\|\varrho_0\|_{L^2(U)})^{1/2}.\end{aligned}$$

We now have strong convergence of $\left(u_n\right)_{n=1}^{\infty}$, by showing it is a Cauchy sequence in a complete metric space.\

\[lem:rhon\_is\_cauchy\]

Let $T>0$ and suppose $\{u_n\}_{n\geq1}$ satisfies with $\varrho_0\in C^\infty(U)$. Then there exists $\varrho\in L^1([0,T];L^1(U))$ such that $u_n\to \varrho$ in $L^1([0,T];L^1(U))$.

Lastly we have the uniform estimate on the limit point $\varrho(\vec{r},t)$ in terms of the initial data $\varrho_0$.\

\[lem:weak\_conv\_results\]

One has $\varrho\in L^2([0,T]; H^1(U))\cap L^\infty([0,T]; L^2(U))$ and that $\partial_t\varrho\in L^2([0,T]; H^{-1}(U))$ with the uniform bound $$\begin{aligned}
\label{eq:total_energy_bound}
\|\varrho\|_{L^{\infty}([0,T];L^2(U))}+\|\varrho\|_{L^{2}([0,T];H^1(U))}+\|\partial_t\varrho\|_{L^2([0,T]; H^{-1}(U))}\leq C(T) \|\varrho_{0}\|_{L^2(U)}.\end{aligned}$$ Additionally there exists a subsequence $\{u_{n_k}\}_{k\geq 1}$ such that $$\begin{aligned}
u_{n_k}\rightharpoonup \varrho& \quad \text{ in } L^2([0,T]; H^1(U)),\label{lim:rho_nk_in_H1_c}\\
\partial_tu_{n_k}\rightharpoonup \partial_t\varrho& \quad \text{ in } L^2([0,T]; H^{-1}(U)).\label{lim:partial_trho_nk_in_H_minus_1}\end{aligned}$$ where $\rightharpoonup$ denotes weak convergence.

The nature of convergence of the sequence $\{\varrho_n\}_{n\geq 1}$ as $n\to\infty$ are consolidated into the following result.\

\[cor:summary\_of\_convergence\]

There exists a subsequence $\{u_{n_k}\}_{k\geq 1}\subset \{u_{n}\}_{n\geq 1}$ and a function $\varrho\in L^2([0,T];H^1(U))$ with $\partial_t\varrho\in L^2([0,T];H^{-1}(U))$ such that $$\begin{aligned}
u_n\to \varrho& \quad \text{ in } L^1([0,T]; L^1(U)),\\
u_{n_k}\rightharpoonup \varrho& \quad \text{(weakly) in } L^2([0,T]; H^1(U)),\\
\partial_t u_{n_k}\rightharpoonup \partial_t\varrho& \quad \text{(weakly) in } L^2([0,T]; H^{-1}(U)).\end{aligned}$$

We are now in the position to obtain the existence and uniqueness of weak solutions $\varrho(\vec{r},t)$. First we state a calculus result which will be useful when working with the weak formulation .\

\[lem:calculus\_of\_inner\_product\]

Suppose $\varrho\in L^2([0,T]; H^1(U))$ and $\partial_t\varrho\in L^2([0,T]; H^{-1}(U))$ then the mapping $$\begin{aligned}
t\to \|\varrho(t)\|_{L^2(U)}^2\end{aligned}$$ is absolutely continuous with $$\begin{aligned}
\der[]{t}\|\varrho(t)\|_{L^2(U)}^2 = 2\langle \partial_t\varrho(t),\, \varrho(t)\rangle\end{aligned}$$ for a.e. $t\in [0,T]$.

Since the condition $\Pi[\varrho]\cdot \vec{n} = 0$ on $\partial U \times [0,T]$ guarantees integration by parts without extra terms the proof is identical to the textbook one [@evans2002partial].

We are now in the position to prove existence of the weak solution to .

Existence and Uniqueness.
-------------------------

By using Propositions \[prop:bound\_rhon\], \[prop:un\_H1\_rho0\_bound\] and Lemmas \[lem:rhon\_is\_cauchy\], \[lem:weak\_conv\_results\], \[lem:calculus\_of\_inner\_product\] we may obtain the following theorem.\

(Existence and Uniqueness of Weak Density)\[thm:existence\_and\_uniqueness\]

Let $\varrho_0 \in C^\infty(U)$, $\varrho\geq 0$ and $\int_U \mathrm{d}\vec{r} \varrho_0(\vec{r})=1$. Then there exists a unique weak solution $\varrho\in L^\infty([0,T];L^2(U))\cap L^2([0,T]; H^1(U))$, with $\partial_t\varrho\in L^2([0,T]; H^{-1}(U))$, to equation in the sense with the estimate .

Multiply by $\eta\in L^2([0,T]; H^1(U))$ after setting $n=n_k\in \mathbb{N}$ and integrate over $U_T$ to obtain $$\begin{aligned}
&\int_0^T\mathrm{d}t\, \langle \partial_tu_{n_k},\, \eta(t) \rangle\nonumber\\
& + \int_0^T\mathrm{d}t\,\int_U\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\eta \cdot\bm{D}\,[\nabla_{\vec{r}}u_{n_k}+  u_{n_k}\nabla_{\vec{r}}(\kappa_1V_1+\kappa_2\,V_2\star u_{n_k-1})]=0.\end{aligned}$$ For the transport term we write $$\begin{aligned}
&\int_0^T\mathrm{d}t\, \nabla_{\vec{r}}\eta \cdot u_{n_k}\bm{D}\,\nabla_{\vec{r}}[\kappa_1V_1+\kappa_2\,V_2\star u_{n_k-1}]\nonumber\\
& = \int_0^T\mathrm{d}t\, \nabla_{\vec{r}}\eta \cdot\,(u_{n_k}-\varrho)\bm{D}\,\nabla_{\vec{r}}[\kappa_1V_1+\kappa_2\,V_2\star u_{n_k-1}]\nonumber\\
 &\quad  +\int_0^T\mathrm{d}t\, \nabla_{\vec{r}}\eta \cdot \,\varrho \bm{D}\,\nabla_{\vec{r}}[\kappa_1V_1+\kappa_2\,V_2\star (u_{n_k-1}-\varrho)]
 +\int_0^T\mathrm{d}t\, \nabla_{\vec{r}}\eta\cdot \,\varrho \bm{D}\,\nabla_{\vec{r}}\kappa_2\,V_2\star \varrho.\end{aligned}$$ Note that $u_{n_k}\rightharpoonup\varrho$ in $L^2([0,T];H^1(U))\subset L^2([0,T];L^2(U))$ and $(\nabla_{\vec{r}}\cdot\bm{D})\cdot\nabla_{\vec{r}}[\kappa_1V_1(\vec{r})+\kappa_2\,V_2\star (\varrho_{n_k-1})]$ is uniformly bounded and so $$\begin{aligned}
\int_0^T\mathrm{d}t\, \int_U\mathrm{d}\vec{r}\,\nabla_{\vec{r}}^\top\eta \, (\varrho_{n_k}-\varrho)\bm{D}\,\nabla_{\vec{r}}[\kappa_1V_1+\kappa_2\,V_2\star \varrho_{n_k-1}]\to 0\end{aligned}$$ as $k\to \infty$.

Now by H[ö]{}lder’s inequality one has $$\begin{aligned}
&  \int_0^T\mathrm{d}t\, \nabla_{\vec{r}}\eta\cdot \,\varrho\,\bm{D}\,\nabla_{\vec{r}}\kappa_2(V_2\star (u_{n_k-1}-\varrho))
\leq \mu_{\max}\|\nabla_{\vec{r}}\eta \|_{L^2([0,T]; L^2(U))}\|\nabla_{\vec{r}}V_2\|_{L^{\infty}(U)}\nonumber\\
&\times \left( \int_0^T\mathrm{d}t\,\|u_{n_k-1}(t)-\varrho(t)\|_{L^1(U)}^2\right)^{1/2}\to 0.\end{aligned}$$ Now note that by Lemma \[lem:rhon\_is\_cauchy\], $\|\phi_n\|_{L^1(U)}$ is bounded and therefore $$\begin{aligned}
\int_0^T\mathrm{d}t \|u_{n_k-1}(t)-\varrho(t)\|_{L^1(U)}^2\leq C\int_0^T\mathrm{d}t \|u_{n_k-1}(t)-\varrho(t)\|_{L^1([0,T];L^1(U))}
\to 0.\end{aligned}$$

Therefore we have $$\begin{aligned}
\int_0^T\mathrm{d}t\, \nabla_{\vec{r}}\eta\, \cdot \,u_{n_k}\bm{D}\,\nabla_{\vec{r}}[\kappa_1V_1+\kappa_2\,V_2\star u_{n_k-1}]
\to\int_0^T\mathrm{d}t\, \nabla_{\vec{r}}\eta \cdot\varrho  \bm{D}\,\nabla_{\vec{r}}[\kappa_1V_1+\kappa_2\,V_2\star \varrho]\end{aligned}$$ as $k\to \infty$. By the weak convergence results of Lemma \[lem:weak\_conv\_results\] we have $$\begin{aligned}
\int_0^T\mathrm{d}t\, \langle \partial_t u_{n_k},\, u_{n_k} \rangle &\to \int_0^T\mathrm{d}t\, \langle \partial_t\varrho,\, \varrho \rangle,\nonumber\\
\int_0^T\mathrm{d}t\,\int_U\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\eta \cdot\bm{D}\,\nabla_{\vec{r}}u_{n_k} &\to \int_0^T\mathrm{d}t\,\int_U\mathrm{d}\vec{r}\, \nabla_{\vec{r}}\eta \cdot\bm{D}\,\nabla_{\vec{r}}\varrho\end{aligned}$$ as $k\to \infty$. This establishes existence of weak solution to in the sense . Establishing $\varrho(0)=\varrho_0$ is a routine argument (see [@evans2002partial]).

To prove uniqueness we set $\xi = \varrho_1-\varrho_2$ where $\varrho_1,\varrho_2$ are weak solutions then we have $$\begin{aligned}
&\int_0^T\mathrm{d}t\, \langle \partial_t \xi(t), \, \eta(t)  \rangle\nonumber\\
& +\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, \nabla_{\vec{r}}\eta\cdot \bm{D}\,[\nabla_{\vec{r}}\xi
 +\xi\,\nabla_{\vec{r}}\kappa _1V_1+\kappa _2\varrho_1\nabla_{\vec{r}}\,V_2\star\varrho_1-\kappa_2\varrho_1\nabla_{\vec{r}}\,V_2\star\varrho_2]=0\end{aligned}$$ Adding and subtracting $\int_0^T\mathrm{d}t\,\int_U\mathrm{d}\vec{r}'\,\nabla_{\vec{r}}\eta\cdot \kappa_1\varrho_2\nabla_{\vec{r}}V_2\star \varrho_1$ we find $$\begin{aligned}
&\int_0^T\mathrm{d}t\, \langle \partial_t \xi(t), \, \eta(t)  \rangle + \int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, \nabla_{\vec{r}}\eta\cdot \bm{D}\,\nabla_{\vec{r}}\xi \nonumber\\
&=-\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, \nabla_{\vec{r}}\eta\cdot \bm{D}\,[\xi\,\nabla_{\vec{r}}\kappa _1V_1+\kappa _2\xi\nabla_{\vec{r}}\,V_2\star\varrho_1-\kappa_2\varrho_2\nabla_{\vec{r}}\,V_2\star\xi]\nonumber\\
&\leq \int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, |\nabla_{\vec{r}}\eta\cdot \bm{D}^{1/2} \bm{D}^{1/2}\,[\xi\,\nabla_{\vec{r}}\kappa _1V_1+\kappa _2\xi\nabla_{\vec{r}}\,V_2\star\varrho_1-\kappa_2\varrho_2\nabla_{\vec{r}}\,V_2\star\xi]|.\label{eq:uniqueness_bound_1}\end{aligned}$$ By Young’s inequality we have $$\begin{aligned}
&\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, |\nabla_{\vec{r}}\eta\cdot \bm{D}^{1/2} \bm{D}^{1/2}\,[\xi\,\nabla_{\vec{r}}\kappa _1V_1+\kappa _2\xi\nabla_{\vec{r}}\,V_2\star\varrho_1-\kappa_2\varrho_2\nabla_{\vec{r}}\,V_2\star\xi]|\nonumber\\
&\leq \int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, |\bm{D}^{1/2} \nabla_{\vec{r}}\eta |^2\nonumber\\
& +\tfrac{1}{4}\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, | \bm{D}^{1/2} [\xi\,\nabla_{\vec{r}}\kappa _1V_1+\kappa _2\xi\nabla_{\vec{r}}\,V_2\star\varrho_1-\kappa_2\varrho_2\nabla_{\vec{r}}\,V_2\star\xi ] |^2.\end{aligned}$$ Using the triangle inequality and Young’s inequality we expand the absolute value inside the integral $$\begin{aligned}
&\tfrac{1}{4}\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, |\bm{D}^{1/2} [ \xi\,\nabla_{\vec{r}}\kappa _1V_1+\kappa _2\xi\nabla_{\vec{r}}\,V_2\star\varrho_1-\kappa_2\varrho_2\nabla_{\vec{r}}\,V_2\star\xi ]|^2\nonumber\\
&\leq \tfrac{1}{4}\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\, |\bm{D}^{1/2} \xi\,\nabla_{\vec{r}}\kappa _1V_1|^2
+\kappa _2^2|\bm{D}^{1/2} [ \xi\nabla_{\vec{r}}\,V_2\star\varrho_1-\varrho_2\nabla_{\vec{r}}\,V_2\star\xi ] |^2 \nonumber\\
&\leq \tfrac{1}{4}\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r} \left(  |\bm{D}^{1/2}\xi\,\nabla_{\vec{r}}\kappa _1V_1|^2
+2\kappa _2^2|\bm{D}^{1/2}\xi\nabla_{\vec{r}}\,V_2\star\varrho_1|^2
+2\kappa _2^2|\bm{D}^{1/2}\varrho_2\nabla_{\vec{r}}\,V_2\star\xi|^2\right)\nonumber\\
&\leq \tfrac{\mu_{\max}}{4}\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\,  |\xi\,\nabla_{\vec{r}}\kappa _1V_1|^2+2\kappa _2^2|\xi\nabla_{\vec{r}}\,V_2\star\varrho_1|^2+2\kappa _2^2|\varrho_2\nabla_{\vec{r}}\,V_2\star\xi|^2. \label{eq:uniqueness_bound_2}\end{aligned}$$ Estimating each of these terms, first $$\begin{aligned}
\label{eq:uniqueness_bound_3}
\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\,  |\xi\,\nabla_{\vec{r}}\kappa _1V_1|^2\leq \kappa_1^2\|\nabla_{\vec{r}}V_1\|_{L^\infty(U)}^2\|\xi\|_{L^2([0,T];L^2(U))}.\end{aligned}$$ Second, $$\begin{aligned}
\label{eq:uniqueness_bound_4}
2\kappa _2^2\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\,  |\xi\nabla_{\vec{r}}\,V_2\star\varrho_1|^2 \leq 2\kappa _2^2|U\||\nabla_{\vec{r}}V_2\|_{L^\infty(U)}^2\|\xi\|_{L^2([0,T];L^2(U))},\end{aligned}$$ and third $$\begin{aligned}
&2\kappa _2^2\int_0^T\mathrm{d}t\, \int_U \mathrm{d}\vec{r}\,|\varrho_2\nabla_{\vec{r}}\,V_2\star\xi|^2\nonumber\\
&\leq 2\kappa _2^2|U| \|\varrho_2\|_{L^\infty([0,T];L^2(U))} \|\nabla_{\vec{r}}V_2\|_{L^\infty(U)}^2\|\xi\|_{L^2([0,T];L^2(U))}.\label{eq:uniqueness_bound_5}\end{aligned}$$

Combining , , , , we obtain, after setting $\eta = \xi$, and using boundedness of $\varrho_2$ in terms of its initial data $$\begin{aligned}
\int_0^T\mathrm{d}t\, \langle \partial_t \xi(t),\, \xi(t)\rangle  \leq (C_1(T)+C_2(T)\|\varrho_0\|_{L^2(U)}^2)\|\xi \|_{L^2([0,T];L^2(U))}^2\end{aligned}$$ for some constants $C_1(T)$, $C_2(T)$ dependent on $U$. This holds for all $T$ so it must be the case that $$\begin{aligned}
\der[]{t}\|\xi(t)\|_{L^2(U)}^2\leq (C_1(T)+C_2(T)\|\varrho_0\|_{L^2(U)}^2)\|\xi(t)\|_{L^2(U)}^2\end{aligned}$$ implying by Gr[ö]{}nwall’s lemma that $$\begin{aligned}
\|\xi(t)\|_{L^2(U)} \leq (C_1(T)+C_2(T)\|\varrho_0\|_{L^2(U)}^2)\|\xi(0)\|_{L^2(U)}\end{aligned}$$ a.e. $t\in [0,T]$. However, $\xi(0)\equiv 0$ hence $\|\varrho_1(t)-\varrho_2(t)\|_L^2(U)=0$ for all $t\in [0,T]$.

Strict Positivity of $\varrho$. {#subsec:strict_pos_rho}
-------------------------------

With the existence of weak solutions we may establish positivity of $\varrho$ solving with reference to [@bogachev2015fokker]. In particular since $\boldsymbol{D}$ is positive definite and $\vec{b}$ is uniformly bounded and $$\begin{aligned}
\sup_{\vec{r}\in U}\varrho(\vec{r},t_1)<C \inf_{\vec{r}\in U}\varrho(\vec{r},t_2)\end{aligned}$$ for $0<t_1<t_2<\infty$ and $C$ is a constant depending on $d$ (the dimension) and $\mu_{\max}$. Since $\varrho$ is non-negative for all time we must have $\inf_{\vec{r}\in U}\varrho(\vec{r},t)$ is positive and hence $\varrho$ is positive.\

Discussion & Open Problems {#sec:discussion}
==========================

In this paper, the global asymptotic stability and well-posedness of overdamped DDFT with two-body HI was studied. It was shown that bifurcations occur in DDFT systems with no-flux boundary conditions at an infinite and discrete set of critical energies equal to eigenvalues of the two-body interaction integral operator $\mathcal{R}$. Additionally we have shown that a weak solution to the density with no-flux boundary conditions and strong solution to the flux equation exist and are unique under sensible assumptions on the confining and interaction potentials and initial data $V_1$, $V_2$ and $\varrho(\vec{r},0)$ respectively. Assuming a classical solution to the DDFT we also derived *a priori* convergence estimates in $L^2$ and relative entopy, the latter restricted to convex two-body potentials.

Well-posedness and global asymptotic stability of the phase space equation for the time evolution of $f(\vec{r},\vec{p},t)$ remains open (see [@goddard2012overdamped Proposition 2.1] for the evolution equation for $f(\vec{r},\vec{p},t)$). It is of similar form to the Vlasov equation considered by [@degond1986global] but with Hermite dissipative term and modified nonlocal term in the momentum variable $\vec{p}$ dependent on the HI tensors. To progress further some maximum principles on $f$ solving the linearised version of the phase space equation must be found. Additionally, the existence results on the overdamped equations considered here may be made more regular by routine arguments.

We also note that the present analysis is based on the Smoluchowski equation rigorously derived from the phase space Fokker-Planck equation using homogenisation methods [@goddard2012overdamped]. As an alternative to this, assuming inertia is small altogether, or if one is interested only in very short times to begin with, the system of interacting particles maybe considered solely in configuration space. Only the positions (and not the momenta) of a system of interacting Brownian particles are then taken into account with Smoluchowski equation as in [@rex2009dynamical], and, the underlying Langevin dynamics contain only velocity equations for each particle which are usually written down *a posteriori*. The justification for this is that the momentum distribution is assumed to have a minor role in the dynamical description of the fluid density, and indeed is taken to be irrelevant at the microscopic level. This Brownian approximation may also hold for highly dense suspensions, since in dense Newtonian systems there is a fast transfer of momentum and kinetic energy from the particle collisions, and this effect may be accounted for most efficiently by the bath in the Brownian dynamics with a non constant diffusion tensor. It is known however that the one-body Smoluchowski equation in [@rex2009dynamical] does not equate to equations - which are obtained in the rigorous overdamped limit starting from the Newtonian dynamics. Intuitively this is because the two-body assumption for the HI ($\boldsymbol{\Gamma}$) and mobility ($\boldsymbol{D}$) tensors and the matrix inversion $\boldsymbol{D} = \boldsymbol{\Gamma}^{-1}$ are not commutable operations; even if $\boldsymbol{D}$ is two-body then $\text{det}(\boldsymbol{D})$ is not. A flow chart demonstrating the permitted commutations between various formalisms is included in [@goddard2012overdamped]. The nonequivalence of the two Smoluchowski equations is not considered here, and therefore a natural extension for future work would be to determine the existence, uniqueness and regularity of of the density starting from [@rex2009dynamical] as well as the corresponding conditions for linear stability.

Finally we remark that a well-posedness analysis of DDFT equations of the form to include a hard-sphere contribution to the free energy by fundamental measure theory (FMT) e.g. Rosenfeld [@rosenfeld1989free] or Roth [@roth2010fundamental] would be very interesting.\
————————————-

Classical Linear Parabolic PDE {#sec:classical_paraboliv_pde}
==============================

The first goal is to derive a similar set of estimates as [@chazelle2017well Lemma 3.5, Lemma 3.7]. The standard argument is to set up a sequence of linear parabolic PDEs. Let $U$ be a bounded and open subset of $\mathbb{R}^d$ and set $U_T=U\times (0,\,T]$ for some time $T>0$. Now consider the linear parabolic equation $$\begin{aligned}
\label{eq:linear_classical_parabolic_eqn}
\partial_t u_n-\nabla_{\vec{r}}\cdot[\bm{D}\nabla_{\vec{r}}u_n]=\nabla_{\vec{r}}\cdot[u_n\,\bm{D}\nabla_{\vec{r}}(\kappa _1V_1+\kappa _2V_2\star u_{n-1})].\end{aligned}$$ In general $d$ dimensions we are in the divergence form of the parabolic PDE $$\begin{aligned}
\begin{cases}
\label{eq:ibvp_pde_for_rho_lin}
\qquad\partial_t u_n+Lu_n = 0 \quad \text{ in } U_T,\\
 \Xi [u_n]\cdot\vec{n} = 0 \quad \text{ on } \partial U \times [0, T],\\
\qquad u_n = \varrho_0 \quad \text{ on } U\times \left\lbrace t=0\right\rbrace
\end{cases}\end{aligned}$$ where $\partial U$ is a $C^1$ boundary with unit normal $\vec{n}$. We define $L$ to be the linear differential operator given by $$\begin{aligned}
&Lu_n := -\sum_{ij=1}^d\partial_{r_j}(\bm{D}_{ij}(\vec{r},t)\partial_{r_i}u_n)+\sum_{i=1}^d b_i(\vec{r})\partial_{r_i}u_n+c(\vec{r})u_n,\label{eq:def_of_linear_parabolic_op}\\
& \vec{b}(\vec{r}):=-\bm{D}(\vec{r},t)\nabla_{\vec{r}}(\kappa_1V_1(\vec{r})+\kappa _2 [V_2\star u_{n-1}]), \\
& c(\vec{r}):=-\nabla_{\vec{r}}\cdot (\bm{D}(\vec{r},t)\nabla_{\vec{r}}(\kappa_1V_1(\vec{r})+\kappa _2[V_2\star u_{n-1}]),\\
& \Xi[u_n]:= \bm{D}(\vec{r},t)\,(\nabla_{\vec{r}}u_n + 
u_n\, \nabla_{\vec{r}}(\kappa_1 V_1(\vec{r},t)+ \kappa_2 \,V_2\star u_{n-1})).\end{aligned}$$ Since $\bm{D}(\vec{r},t)$ is assumed to positive definite, there exists $\theta$ for every $\vec{r}$, $\xi$ such that $\xi^\top D(\vec{r},t) \xi\geq \theta |\xi |^2$, therefore the operator $\partial_t + L$ is uniformly parabolic. The Sobolev space of functions that permit the no-flux condition $\Xi [u_n]\cdot \vec{n}$ on $\partial U\times [0,T]$ is $H^1(U)$ which is reflexive, so that $\partial_t u$ interpreted as a bounded linear functional can be paired to an element in $H^1(U)$, and further by the Riez-Representation theorem there exists a unique element from $H^1(U)$ for the pairing. Additionally $H^1(U)$ is separable so that the (unique) weak solution may be approximated by a sequence of smooth functions coming from a countably dense subset.

Weak Formulation.
-----------------

Equation may be recast into weak form. We first introduce the bilinear operator, defined by $$\begin{aligned}
B[u,v;t] := \int_U\mathrm{d}\vec{r}\,   \nabla_{\vec{r}}v\cdot\bm{D}\,\nabla_{\vec{r}}u+\int_U\mathrm{d}\vec{r}\,\vec{b}(\vec{r})\cdot\nabla_{\vec{r}}u\, v+\int_U\mathrm{d}\vec{r}\,c(\vec{r})u\, v\end{aligned}$$ for $u,v\in H^1(U)$ and a.e.  $0\leq t \leq T$. We regard $u$ as a mapping $[\mathfrak{u}(t)](\vec{r}):=u(\vec{r},t)$ from the time interval $[0,T]$ to the function space $H^1(U)$. Now fixing $v \in H^{1}(U)$ we multiply by $v$ and integrate by parts to obtain the weak formulation $$\begin{aligned}
\label{eq:weak_form_for_un}
(\partial_t\mathfrak{u},\,v)+B[\mathfrak{u},\,v;t] = 0\end{aligned}$$ for each $0\leq t\leq T$ with $(\,,\,)$ denoting inner product in $L^2(U)$.

Existence.
----------

The method to establish weak solution for the indexed problem is a textbook one. The method is described as follows. Fix $n$ then the evolution equation is a uniformly parabolic PDE for the unknown $u_n = u$. One now expands $\mathfrak{u} = \mathfrak{u}^m$ in a linear combination of $m$ eigenvectors of the operator $-\nabla_{\vec{r}}\cdot( \bm{D}(\vec{r})\nabla_{\vec{r}}w_k)$ for finite dimensional approximation to $u$. Since $\bm{D}_{ij}(\cdot)$ is a compact and symmetric operator then the eigenfunctions $w_k$ form an orthonormal basis of $L^2(U)$ with $w_k\in H^1(U)$. Thus $\mathfrak{u}^m$ is projected onto the finite dimensional subspace spanned by $\left\lbrace w_k\right\rbrace_{k=1}^m$. The standard existence theory of ODEs (the Carath[é]{}odory conditions with the Cauchy–Picard theorem) gives existence of weak solutions $\mathfrak{u}^m$ as expanded in the functions $\left\lbrace w_k\right\rbrace_{k=1}^m$ on a finite dimensional subspace of $H^1(U)$. All that remains is to pass to the limit $m\to \infty$ to realise the result in $H^1(U)$. To do this energy estimates are required on $\mathfrak{u}^m$, these are routine calculations except in the textbooks they are done for simpler boundary condition choices (homogeneous Dirichlet or periodic) and make use of Poincare’s inequality (holding only for $H^1_0$ functions).

The calculations are similar and for the present boundary condition choice, a weaker Poincar[é]{}$-$Wirtinger inequality is used through out to obtain $$\begin{aligned}
&\max_{0\leq t\leq T}\|\mathfrak{u}^m(t)\|_{L^2(U)}+\|\mathfrak{u}^m\|_{L^2([0,T]; H^{1}    (U))}\nonumber\\
&+\|\pder[]{t}\mathfrak{u}^m\|_{L^2([0,T]; H^{-1}(U))}
\leq c_1\|u_{0}\|_{L^2(U)}+ c_2.\label{eq:uniform_bound_on_three_norms}\end{aligned}$$ where $c_1$, $c_2$ are constants dependent on $T$ and $U$ and $\mu_{\min},\mu_{\max}$. Note that the left hand side of forms a bounded sequence in $\mathbb{R}$ and by the Bolzano–Weierstrass theorem there exists a convergent subsequence $\{\mathfrak{u}^{ m_l}\}_{l\geq 1}\subset\{\mathfrak{u}^{ m}\}_{m\geq 1}$. In particular there exists $\mathfrak{u}$ such that $$\begin{aligned}
\begin{split}
\mathfrak{u}^{ m_l} \rightharpoonup \mathfrak{u}\quad &\text{ weakly in } L^2([0,T]; H^1(U)),\\
\partial_t\mathfrak{u}^{ m_l} \rightharpoonup \mathfrak{u}' \quad &\text{ weakly in } L^2([0,T]; H^{-1}(U)).
\end{split}\end{aligned}$$ Note of course that $\mathfrak{u} = \mathfrak{u}_n$, but we have not yet established existence of weak solution to the full nonlinear Smoluchowski equation . This result establishes existence of weak solution for the parabolic equation for every index $n$. Now since $L^2([0,T]; H^1(U))$ is separable, and weak solutions currently only exist in a finite dimensional subspace of $H^1(U)$, it makes sense to choose a test function $\bm{\phi}\in C^1([0,T]; H^1(U))\subset L^2([0,T]; H^1(U))$. We may therefore write $$\begin{aligned}
\int_0^T\mathrm{d}t\,\langle\partial_t\mathfrak{u}^{ m} ,\,\bm{\phi}^N\rangle +B[\mathfrak{u}^m,\,\bm{\phi}^N;t] =0\end{aligned}$$ for $\bm{\phi}^N=\sum_{k=1}^Nd^k(t)w_k$. Making the choice $N\leq m$ and letting $N\to\infty$ one obtains $$\begin{aligned}
\int_0^T\mathrm{d}t\,\langle\partial_t\mathfrak{u} ,\,\bm{\phi}^{\infty}\rangle +B[\mathfrak{u},\,\bm{\phi}^\infty;t] =0\end{aligned}$$ for any function $\bm{\phi}^\infty\in L^2([0,T]; H^1(U))$ since $\phi^N$ are dense in $L^2([0,T]; H^1(U))$. Now since $\bm{\phi}^\infty$ is arbitrary we obtain $$\begin{aligned}
\langle\partial_t\mathfrak{u} ,\,\phi\rangle +B[\mathfrak{u},\,\phi;t] = 0\end{aligned}$$ for an arbitrary $\phi\in H^1(U)$. Hence the criteria of weak solution is satisfied.

Uniqueness.
-----------

To show uniqueness we argue by contradiction that there exists two weak solutions solutions. By linearity, their difference $\bm{\chi}$ is a weak solution of with $\chi_0\equiv 0$, for $\chi_0$ initial data. Then as it is a weak solution, we may test $\bm{\chi}$ against itself $$\begin{aligned}
\langle\partial_t\bm{\chi} ,\,\bm{\chi}\rangle +B[\bm{\chi},\,\bm{\chi};t] \equiv 0\end{aligned}$$ giving $$\begin{aligned}
\tfrac{1}{2}\der[]{t}(\|\bm{\chi}(t)\|^2_{L^2(U)})+B[\bm{\chi},\,\bm{\chi};t] =0\end{aligned}$$ but $B[\bm{\chi},\,\bm{\chi};t]\geq -c_{7}\|\bm{\chi}(t)\|^2_{L^2(U)}$ which may be obtained by the following estimate $$\begin{aligned}
\label{eq:second_bound_on_B[u,u]}
c_5\|\mathfrak{u}^m-c\|_{H^{1}(U)}^2\leq B[\mathfrak{u}^m,\mathfrak{u}^m]+c_6\|\mathfrak{u}^m\|^2_{L^2(U)}.\end{aligned}$$ and hence by Gr[ö]{}nwall $$\|\bm{\chi}(t)\|^2_{L^2(U)}\leq c_{7}(t)\|\chi_0\|^2_{L^2(U)}=0$$ and $\bm{\chi} =0 $ for a.e. $\vec{r}\in U$ for every $0\leq t \leq T$. We have established the existence and uniqueness of the weak solution to the linear parabolic equation and may apply this to an iteration problem on .\

Nomenclature {#app:nomenclature}
============

  ------------------------------------------------- -- -----------------------------------------------------------------------------------------------------------------------------------
  $\alpha$                                             Mass diffusivity
  $\beta_i^{(\kappa_2)}$, $\theta_i^{(\kappa_2)}$      Expansion coefficients used in Theorem \[thm:epsilon\_as \_a\_fn\_of\_lambda\]
  $\beta_{n}^{-1}$                                     Eigenvalues of $\mathcal{R}$ from Definition \[def:def\_of\_R\_op\]
  $\gamma$                                             Friction coefficient
  $\gamma^{(\kappa_2)}_k$                              Eigenvalues of $\mathcal{A}_{\kappa_2}$ defined in
  $\delta_n$, $\phi_n$, $\psi_n$                       Expansion coefficients used in Theorem \[thm:eigenfn\_expansion\_of\_flux\]
  $\epsilon$                                           Small nondimensional parameter
  $\zeta^a_i$                                          Gaussian white noise process
  $\kappa_1$                                           Nondimensional external potential strength
  $\kappa_2$                                           Nondimensional two body potential strength
  $\kappa_{2_\sharp}$                                  Point of critical stability defined in
  $\lambda$                                            Scalar value of $\mathbb{C}$ used in Lemma \[lem:flux\_operator\_is\_compact\_self\_adjoint\]
  $\lambda(\kappa_2)$                                  Eigenvalue of $\mathcal{L}_{1}$ used in Theorem \[thm:epsilon\_as \_a\_fn\_of\_lambda\]
  $\mu$                                                Probability measure
  $\mu_c$                                              Chemical potential
  $\mu_i$                                              Eigenvalues of $\bm{D}$ defined in
  $\mu_{\max}$, $\mu_{\min}$                           Largest and smallest eigenvalues of $\bm{D}$ respectively
  $\varrho$                                            Density
  $\varrho_0$                                          Initial density
  $\varrho_n$                                          $n$-particle configuration space distribution function, for $n\geq 2$
  $\varrho_\infty$                                     Unique equilibrium density
  $\tau$                                               Characteristic time scale
  $\phi_n$                                             Eigenvalues of $\bm{1}+\mathcal{Z}_1^\varrho -\lambda \mathcal{Z}_2^\varrho$ used in Theorem \[thm:eigenfn\_expansion\_of\_flux\]
  $\chi_\epsilon$                                      Convex approximation to $|\cdot |$ in Definition \[def:def\_of\_chi\_approx\_abs\]
  ------------------------------------------------- -- -----------------------------------------------------------------------------------------------------------------------------------

  : [**[Upper Case Greek]{}**]{}

\
\

  -------------------- -- ---------------------------------------------------------
  $\bm{\Gamma}$           $3N\times 3N$ friction tensor
  $\bm{\Gamma}_{ij}$      $3\times 3$ block matrices of $\bm{\Gamma}$
  $\Pi$, $\Pi_1$          Nonlinear and linear boundary operators respectively, ,
  -------------------- -- ---------------------------------------------------------

  : [**[Upper Case Greek]{}**]{}

  --------------------------------- -- ------------------------------------------------------------------------------------------------------------------------------------
  $\vec{a}$                            Flux
  $c_{ls}$                             Log–Sobolev constant
  $c_{pw}$                             Poincar[é]{}–Wirtinger constant
  $d$                                  Dimension number
  $\vec{e}_i$                          Eigenvectors of $\bm{D}$ defined in
  $f(\vec{r},\vec{p},t)$               Phase space density
  $\vec{f}_i$                          Gaussian white noise vector
  $g(\vec{r},\vec{r}',[\varrho])$      Correlation function
  $k_B$                                Boltzmann constant
  $m$                                  Mass of particle $i$
  $\vec{p}_i$                          Momentum vector of particle $i$
  $\vec{r}_i$                          Position vector of particle $i$
  $t$                                  Time
  $u_n$                                Eigenfunctions of $\mathcal{R}$ from Definition \[def:def\_of\_R\_op\]
  $v^{(\kappa_2)}_k$                   Eigenfunctions of $\mathcal{A}_{\kappa_2}$ defined in
  $w^{(\kappa_2)}$                     Eigenfunction of $\mathcal{L}_1$ used in Theorem \[thm:epsilon\_as \_a\_fn\_of\_lambda\]
  $\vec{w}_n$                          Eigenfunctions of $\bm{1}+\mathcal{Z}_1^\varrho-\lambda\mathcal{Z}_2^\varrho$ used in Theorem \[thm:eigenfn\_expansion\_of\_flux\]
  --------------------------------- -- ------------------------------------------------------------------------------------------------------------------------------------

  : [**[Upper Case Roman]{}**]{}

\
\

  ---------------------------------- -- ---------------------------------------------------------------------------
  $\bm{1}$                              $3\times 3$ identity matrix
  $\bm{1} + \mathcal{Z}_1^\varrho$      Local operator on acting on $\vec{a}$ defined in
  $\mathrm{A}$                          Characteristic flux scale
  $\bm{A}([\vec{a}],t)$                 Advection tensor defined in
  $\mathcal{A}_\varrho$                 The operator $(\bm{1} + \mathcal{Z}_1^\varrho)^{-1}\mathcal{Z}_2^\varrho$
  $\mathcal{A}_{\kappa_2}$              Differential operator defined in
  $\boldsymbol{B}$                      $\left(m k_B T \boldsymbol{\Gamma}\right)^{1/2}$
  $\mathcal{B}$                         Nonlocal operator
  $C(T)$                                Constant dependent on the final time $T$
  $\bm{D}(\vec{r},[\varrho],t)$         Diffusion tensor defined in
  $F$                                   Nonlinear map defined in
  $\mathcal{F}$                         Free energy functional defined in
  $Fr$                                  Froude number
  $\mathcal{F}_H$                       Helmholtz free energy functional defined in
  $\mathcal{H}$                         Relative entropy functional defined in
  $\mathrm{L}$                          Characteristic length scale
  $\mathcal{L}_1$                       Linearised nonlocal differential operator defined in
  ---------------------------------- -- ---------------------------------------------------------------------------

  : [**[Upper Case Roman]{}**]{}

  --------------------------------- -- -------------------------------------------------------------------------------------------------------------------------------------
  $N$                                  Number of particles
  $Pe$                                 P[é]{}clet number
  $\mathcal{R}$                        Nonlocal operator from Definition \[def:def\_of\_R\_op\]
  $\mathrm{T}$                         Temperature
  $T$                                  Final time
  $\mathcal{T}$                        Nonlocal operator defined in
  $U$                                  Spatial domain
  $U_T$                                Space-time cylinder $U\times [0,T]$
  $\mathrm{U}$                         Characteristic velocity scale
  $V$                                  Potential
  $V_1$, $V_2$, $V_n$                  External, two body and $n$-body potentials respectively
  $\mathcal{X}_\varrho$                Inverse operator $(\bm{1}+\mathcal{Z}_1^\varrho+\mathcal{Z}_2^\varrho)^{-1}$ used in Theorem \[thm:association \_of\_free\_energy\]
  $Z$                                  Normalisation constant
  $\bm{Z}_1(\vec{r}_1,\vec{r}_2)$      Diagonal two body HI tensor
  $\bm{Z}_2(\vec{r}_1,\vec{r}_2)$      Off-diagonal two body HI tensor
  $\mathcal{Z}_2^\varrho$              Nonlocal operator acting on $\vec{a}$ defined in
  --------------------------------- -- -------------------------------------------------------------------------------------------------------------------------------------

  : [**[Sets and Mathematical Symbols]{}**]{}

\
\

  ------------------------- -- ----------------------------------------------------------------------------------------
  $L^2(U,\varrho^{-1})$        Weighted $L^2(U)$ space
  $P_{ac}(U)$                  Space of absolutely continuous probability densities supported on $U$
  $P_{ac}^{+}(U)$              $P_{ac}(U)$ restricted to strictly positive functions
  $\text{Tr}$                  Trace operator
  $\top$                       Transpose
  $u\star v$                   Convolution of two functions $\int_U\mathrm{d}\vec{r}\,u(\vec{r}-\vec{r}')v(\vec{r}')$
  $\vec{u}\otimes\vec{v}$      Outer product / dyadic of two vectors $\vec{u}\vec{v}^\top$
  ------------------------- -- ----------------------------------------------------------------------------------------

  : [**[Sets and Mathematical Symbols]{}**]{}

\[bib:References\]

[^1]: School of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, Edinburgh EH9 3FD, UK. Corresponding author: R. D. Mills-Williams (r.mills@ed.ac.uk).

[^2]: Department of Mathematics, Imperial College London, London SW7 2AZ, UK.

[^3]: BDG would like to acknowledge support from EPSRC EP/L025159/1. RDMW is grateful to EPSRC for PhD funding. GAP was supported by the EPSRC through grant numbers EP/P031587/1, EP/L024926/1, and EP/L020564/1. This research was funded in part by JPMorgan Chase & Co. Any views or opinions expressed herein are solely those of the authors listed, and may differ from the views and opinions expressed by JPMorgan Chase & Co. or its affiliates. This material is not a product of the Research Department of J.P. Morgan Securities LLC. This material does not constitute a solicitation or offer in any jurisdiction.
