# The conductivity of strong electrolytes from stochastic density functional theory

Vincent Démery, David S. Dean

<figure><img src="image_1.png" /></figure>

## To cite this version:

Vincent Démery, David S. Dean. The conductivity of strong electrolytes from stochastic density functional theory. Journal of Statistical Mechanics: Theory and Experiment, 2016, 2016 (2), pp.023106 (1-16). 10.1088/1742-5468/2016/02/023106 hal-01283991

HAL Id: hal-01283991

https://hal.science/hal-01283991

Submitted on 7 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# The conductivity of strong electrolytes from stochastic density functional theory

V Démery<fnref n="1" /> and D S Dean<fnref n="2" />

E-mail: vincent.demery@espci.fr and david.dean@u-bordeaux.fr

**Abstract.** Stochastic density functional theory is applied to analyze the conductivity of strong two species electrolytes at arbitrary field strengths. The corresponding stochastic equations for the density of the electrolyte species are solved by linearizing them about the mean density of ionic species, yielding an effective Gaussian theory. The non-equilibrium density–density correlation functions are computed and the conductivity of the electrolyte is deduced. In the bulk, our results give a simple derivation of the results of Onsager and coworkers who used very different methods. The method developed here can also be used to study electrolytes confined in one and two dimensions and interacting via either the three dimensional Coulomb interaction or the Coulomb interaction corresponding to that dimension of space.

**Keywords:** correlation functions, transport processes/heat transfer (theory), ionic liquids, electrolyte solutions, liquid metals and charged fluids (experiments)

<footnote n="1">¹ PCT, UMR Gulliver 7083, ESPCI and CNRS, 10 rue Vauquelin, 75005 Paris, France</footnote>

<footnote n="2">² Laboratoire Ondes et Matière d'Aquitaine (LOMA), CNRS, UMR 5798/ Université de Bordeaux, F-33400 Talence, France</footnote>
# Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Model</td><td>3</td></tr><tr><td>3. General computation using the linearized Stochastic Density Functional Theory</td><td>5</td></tr><tr><td>4. Zero field conductivity for purely electrostatic interactions</td><td>7</td></tr><tr><td>4.1. General results for a purely electrostatic interaction</td><td>7</td></tr><tr><td>4.2. Conductivity for ions moving in homogeneous space</td><td>8</td></tr><tr><td>4.3. Conductivity for ions confined to d = 1 or 2 and interacting with the 3d Green function</td><td>10</td></tr><tr><td>4.3.1. Confinement to a plane</td><td>10</td></tr><tr><td>4.3.2. Confinement to a line</td><td>11</td></tr><tr><td>5. Field dependence for monovalent salts</td><td>12</td></tr><tr><td>5.1. Correlation function</td><td>12</td></tr><tr><td>5.2. Conductivity</td><td>15</td></tr><tr><td>6. Conclusion</td><td>16</td></tr><tr><td>References</td><td>16</td></tr></table>

## 1. Introduction

Onsager's study of the conductivity of strong electrolytes is famous as one of his first major scientific contributions. In this work he showed how the theory developed by Debye and Hückel [1] should be modified to take into account the effect of Brownian motion on the ions [2]. For an electroneutral system with two ionic species $\alpha \in \{+, -\}$ with ionic charges $qz_\alpha$ ($z_\alpha$ denoting the signed valency), mobilities $\kappa_\alpha$ and average densities $\bar{\rho}_\alpha$, the naive formula for the conductivity, which neglects the interaction between the ions, is

$$ \sigma_0 = q^2 \sum_{\alpha} z_{\alpha}^2 \kappa_{\alpha} \bar{\rho}_{\alpha}. \qquad (1) $$

Debye and Hückel [1] computed a correction to this bare conductivity due to the flow induced on the solvent by the charge distribution about an ion. Onsager later pointed out how this result is changed (essentially additively at the level of approximation employed) by the Brownian motion of the ions and their mutual interaction [2]. Indeed it is clear that even in the absence of the solvent the conductivity will be modified by correlations between the ions: an ion moves in an electric field which is the sum of the uniform applied field and the field generated by the other ions in the system.
Onsager then went on to study the conductivity of strong and weak electrolytes in the presence of a finite applied field [3]. For strong electrolytes the modification due to the finite electric field is that the Debye screening of charges is killed off by the field and thus the conductivity increases—the so-called Wien effect [4]. For weak electrolytes, which consists of free anions and cations along with bound Bjerrum pairs, the external field has the effect of shifting the chemical equilibrium between ions and Bjerrum pairs (by pulling the pairs apart) [5]. This shift increases the density of charge carriers, leading to an increase in conductivity. The validity of Onsager’s calculation for weak electrolytes was verified recently, for a solvent free system, by numerical simulations of the second Wien effect [6]. Interestingly, Onsager returned regularly to this problem throughout his career, he compared it carefully with experimental results and also developed new mathematics to facilitate the study of these systems [3].

In this paper we revisit the problem of the conductivity of strong electrolytes for Brownian systems (we do not consider the effect of a solvent here) by analysing the stochastic density functional theory (SDFT) for the density field of the interacting particles. The SDFT in question was first written down by Kawasaki [7] on phenomenological grounds but was later shown formally to be the equation of evolution for the density field of interacting Brownian particles [8]. The SDFT is analytically intratable when the particles interact, and very difficult to analyse even when there are no interactions [9], however it can be used as the starting point for various approximate theories, for example mode coupling theory [10]. Recently it was shown that the linearised SDFT for charged particle systems yields a dynamical theory which has the Debye–Hückel theory as its static limit [11]. This has allowed the study of the dynamics of the thermal Casimir effect between adjacent plates containing Brownian charges. The same type of linearization approximation can also be used to compute the effective diffusion constant of a Brownian tracer particle in interaction with a bath of identical particles [12]. Remarkably, the result obtained from this computation is identical to that obtained by an arduous one-loop renormalisation group analysis of the full $N$ particle Fokker–Planck equation [13]. Finally we remark that the SDFT has even been used to derive integration results in random matrix theory [14].

Here we will show how Onsager’s results for electrolyte conductivity can be derived in a relatively straightforward fashion. As well as considering the case of a globally electroneutral electrolyte with two species of oppositely charged ions, we also deal with mobile ions that are not as an ensemble electroneutral. In this case, the mobile ion’s total charge is neutralized by a uniform immobile background charge—a so-called jelium model. We also generalize these results for the conductivity to lower dimensions $d$ with the corresponding $d$-dimensional Coulomb interaction. Finally, our results can be extended to systems with reduced dimension, for example ions confined to a two dimensional surface.

## 2. Model

We will consider a system containing two species of Brownian particles which interact via pairwise interactions and which are acted on by a uniform external field. The basic
formalism developed here can be applied to systems with an arbitrary number of species, however the reader will see below that on increasing the number of species the calculations become rather cumbersome. We will ultimately apply this model to the case of electrostatic interactions between charged particles in the presence of an applied field. The formalism developed however allows the possibility of including additional interactions between the ions. Let each species be denoted by an index $\alpha$ taking the label + or - depending if the ion is a cation or anion respectively. We denote by $\rho_{\alpha}(\mathbf{x}, t)$ the local density operator of the species $\alpha$ and $\kappa_{\alpha}$ its mobility. We denote by $V_{\alpha\beta}$ the pairwise interaction between species $\alpha$ and $\beta$. We will ultimately decompose the pairwise interaction $V_{\alpha\beta}$ into an electrostatic and non-electrostatic contribution as

$$V_{\alpha\beta}(\mathbf{x}) = q^2 z_{\alpha} z_{\beta} G_0(\mathbf{x}) + U_{\alpha\beta}(\mathbf{x}). \qquad (2)$$

Here $G_0(\mathbf{x})$ denotes the Green's function for the electrostatic interaction in the system, $q$ is the basic unit of charge and $z_{\alpha}$ denotes the valency of the ions of type $\alpha$. The interaction $U_{\alpha\beta}(\mathbf{x})$ denotes interactions other than the direct electrostatic interaction that may be present in the system, for example hard and soft core interactions. Note that in the case where there is a neutralizing uniform background charge (jellium model) there is no contribution to the electric field due to this uniform charge.

At temperature $T$, the density fields obey the SDFT [8]:

$$\partial_t \rho_{\alpha} = -\nabla \cdot \mathbf{j}_{\alpha}, \qquad (3)$$

$$\mathbf{j}_{\alpha} = -T\kappa_{\alpha}\nabla\rho_{\alpha} + \kappa_{\alpha}\rho_{\alpha}\mathbf{f}_{\alpha} + (\kappa_{\alpha}\rho_{\alpha})^{1/2}\eta_{\alpha}, \qquad (4)$$

where $\eta_{\alpha}(\mathbf{x}, t)$ is a Gaussian white noise with correlation function

$$\langle \eta_{\alpha}(\mathbf{x}, t) \eta_{\beta}(\mathbf{x}', t')^{\mathrm{T}} \rangle = 2 T \delta_{\alpha\beta} \delta(\mathbf{x} - \mathbf{x}') \delta(t - t'). \qquad (5)$$

In the above we have introduced the force $\mathbf{f}_{\alpha}(\mathbf{x}, t)$ generated by the applied external field $\mathbf{E}$ and the interactions between the ions:

$$\mathbf{f}_{\alpha} = z_{\alpha} q \mathbf{E} - \sum_{\beta} \nabla V_{\alpha\beta} * \rho_{\beta}, \qquad (6)$$

where * denotes the convolution over spatial variables.

The average electrical current is given by $\langle \mathbf{J} \rangle = q \sum_{\alpha} z_{\alpha} \langle \mathbf{j}_{\alpha} \rangle$ and from this the conductivity $\sigma$ is defined via $\langle \mathbf{J} \rangle = \sigma \mathbf{E}$. Inserting equations (4) and (6) in the definition of the current gives

$$\langle \mathbf{J} \rangle = q^2 \left( \sum_{\alpha} z_{\alpha}^2 \kappa_{\alpha} \bar{\rho}_{\alpha} \right) \mathbf{E} - q \sum_{\alpha, \beta} z_{\alpha} \kappa_{\alpha} \langle \rho_{\alpha} \nabla V_{\alpha\beta} * \rho_{\beta} \rangle. \qquad (7)$$

To rewrite the second term, we introduce the density fluctuations and their stationary equal time correlation as

$$n_{\alpha}(\mathbf{x}, t) = \rho_{\alpha}(\mathbf{x}, t) - \bar{\rho}_{\alpha}, \qquad (8)$$

$$C_{\alpha\beta}(\mathbf{x}) = \langle n_{\alpha}(\mathbf{x}, t) n_{\beta}(0, t) \rangle. \qquad (9)$$

The average electric current is thus
$$ \langle \mathbf{J} \rangle = \sigma_0 \mathbf{E} - q \sum_{\alpha, \beta} z_{\alpha} \kappa_{\alpha} \int C_{\alpha\beta}(\mathbf{x}) \nabla V_{\alpha\beta}(\mathbf{x}) d\mathbf{x}. \qquad (10) $$

Neglecting the correlations between the ions, only the electric field contribution in the second term of equation (4) does not average to zero and the conductivity is given by equation (1). In order to compute the correction to the bare conductivity $\sigma_0$, we need to evaluate the correlations of the density fluctuations.

### 3. General computation using the linearized Stochastic Density Functional Theory

We have to use the SDFT (equations (3)–(6)) to compute the stationary correlation function $C_{\alpha\beta}(\mathbf{x})$. Unfortunately these equations are not linear and contain multiplicative noise, thus the correlation function cannot be calculated in general. This difficulty is circumvented by assuming small density variations, $n_\alpha(\mathbf{x}, t) \ll \bar{\rho}_\alpha$, which allows one to linearize the SDFT [11, 12]. The dynamics of $n_\alpha(\mathbf{x}, t)$ become, to the lowest order

$$ \partial_t n_\alpha = T\kappa_\alpha \nabla^2 n_\alpha - \kappa_\alpha z_\alpha q \mathbf{E} \cdot \nabla n_\alpha + \kappa_\alpha \bar{\rho}_\alpha \nabla^2 \left[ \sum_\beta V_{\alpha\beta} * n_\beta \right] + (\kappa_\alpha \bar{\rho}_\alpha)^{1/2} \nabla \cdot \mathbf{\eta}_\alpha. \qquad (11) $$

The dynamics takes a simpler form in Fourier space, where $\tilde{n}_\alpha(\mathbf{k}, t) = \int e^{-i\mathbf{k} \cdot \mathbf{x}} n_\alpha(\mathbf{x}) d\mathbf{x}$:
$$ \partial_t \tilde{n}_\alpha(\mathbf{k}) = -\kappa_\alpha (Tk^2 + iz_\alpha q \mathbf{E} \cdot \mathbf{k}) \tilde{n}_\alpha(\mathbf{k}) - \kappa_\alpha \bar{\rho}_\alpha k^2 \sum_\beta \tilde{V}_{\alpha\beta}(\mathbf{k}) \tilde{n}_\beta(\mathbf{k}) + \tilde{\xi}_\alpha(\mathbf{k}, t), \qquad (12) $$

where the Gaussian noise $\tilde{\xi}_\alpha(\mathbf{k}, t)$ has a correlation function

$$ \langle \tilde{\xi}_{\alpha}(\mathbf{k}, t) \tilde{\xi}_{\beta}(\mathbf{k}', t') \rangle = 2(2\pi)^d T \kappa_{\alpha} \bar{\rho}_{\alpha} k^2 \delta_{\alpha\beta} \delta(t-t') \delta(\mathbf{k}+\mathbf{k}'). \qquad (13) $$

We can now rewrite the dynamics in terms of the two component vector

$$ N(\mathbf{x}, t) = \begin{pmatrix} n_+(\mathbf{x}, t) \\ n_-(\mathbf{x}, t) \end{pmatrix}, \qquad (14) $$

using the matrices

$$ \tilde{R}(\mathbf{k}) = k^2 \begin{pmatrix} \bar{\rho}_+ \kappa_+ & 0 \\ 0 & \bar{\rho}_- \kappa_- \end{pmatrix}, \qquad (15) $$

$$ \tilde{A}(\mathbf{k}) = T \begin{pmatrix} \frac{1}{\bar{\rho}_+} \left( 1 + i \frac{z_+ q \mathbf{E} \cdot \mathbf{k}}{Tk^2} \right) + \frac{\tilde{V}_{++}(k)}{T} & \frac{\tilde{V}_{+-}(k)}{T} \\ \frac{\tilde{V}_{+-}(k)}{T} & \frac{1}{\bar{\rho}_-} \left( 1 + i \frac{z_- q \mathbf{E} \cdot \mathbf{k}}{Tk^2} \right) + \frac{\tilde{V}_{--}(k)}{T} \end{pmatrix}. \qquad (16) $$

The evolution equation (12) now reads

$$ \partial_t \tilde{N} = -\tilde{R} \tilde{A} \tilde{N} + \Xi, \qquad (17) $$
where the correlation function of the Gaussian noise $\tilde{\Xi}(\mathbf{k}, t)$ is given by

$$ \langle \tilde{\Xi}(\mathbf{k}, t) \tilde{\Xi}(\mathbf{k}', t')^{\mathrm{T}} \rangle = 2(2\pi)^d T\tilde{R}(\mathbf{k}) \delta(\mathbf{k} + \mathbf{k}') \delta(t - t'). \quad (18) $$

The density fluctuation correlation matrix $C_{\alpha\beta}(\mathbf{x})$ reads

$$ C(\mathbf{x}) = \langle N(\mathbf{x})N(0)^{\mathrm{T}} \rangle, \quad (19) $$

and in Fourier space, the correlation is given by

$$ \langle \tilde{N}(\mathbf{k}) \tilde{N}(\mathbf{k}')^{\mathrm{T}} \rangle = \int e^{-i\mathbf{k}\cdot\mathbf{x}-i\mathbf{k}'\cdot\mathbf{x}'} C(\mathbf{x}-\mathbf{x}') d\mathbf{x} d\mathbf{x}' \quad (20) $$

$$ = \int e^{-i(\mathbf{k}+\mathbf{k}')\cdot\mathbf{x}+i\mathbf{k}'\cdot\mathbf{u}} C(\mathbf{u}) d\mathbf{x} d\mathbf{u} \quad (21) $$

$$ = (2\pi)^d \delta(\mathbf{k} + \mathbf{k}') \tilde{C}(\mathbf{k}). \quad (22) $$

In the stationary regime, $\tilde{C}$ satisfies [15]

$$ \tilde{R}\tilde{A}\tilde{C} + \tilde{C}\tilde{A}^*\tilde{R} = 2T\tilde{R}. \quad (23) $$

where we have used that $\tilde{R}(-\mathbf{k})^{\mathrm{T}} = \tilde{R}(\mathbf{k})$ and $\tilde{A}(-\mathbf{k})^{\mathrm{T}} = \tilde{A}(\mathbf{k})^*$. Without an electric field, the system would be in thermal equilibrium and $\tilde{A}$ would be self-adjoint; in this case one finds $\tilde{C} = T\tilde{A}^{-1}$, and one recovers the standard Debye-Hückel approximation for the density fluctuation correlation function. In the presence of an electric field, we have to solve the system of equations (23). The number of equations to solve is $M(M+1)/2$, where $M$ is the number of species; in our case of two species we have only three linear equations to solve.

Writing $\rho_{\pm}\kappa_{\pm} = r_{\pm}$ and $\tilde{A} = T\begin{pmatrix} a & b \\ b & c \end{pmatrix}$, solving for the components of the correlation function leads to

$$ \tilde{C} = \frac{2}{(a + a^*)(c + c^*)|r_+ a + r_- c^*|^2 - b^2[r_+(a + a^*) + r_-(c + c^*)]^2} \times \begin{pmatrix} (c + c^*)|r_+ a + r_- c^*|^2 & -b(r_+ a^* + r_- c)[r_+(a + a^*) + r_-(c + c^*)] \\ -b(r_+ a + r_- c^*)[r_+(a + a^*) + r_-(c + c^*)] & (a + a^*)|r_+ a + r_- c^*|^2 \end{pmatrix}. \quad (24) $$

Notice that the correlation function's off diagonal components satisfy $\tilde{C}_{+-}(\mathbf{k}) = \tilde{C}_{-+}^*(\mathbf{k})$, or equivalently $\tilde{C}_{+-}(\mathbf{k}) = \tilde{C}_{-+}(-\mathbf{k})$, which in real space corresponds to $C_{+-}(\mathbf{x}) = C_{-+}(-\mathbf{x})$, which is a symmetry condition pointed out by Onsager. This is because the symmetry $\mathbf{x} \to -\mathbf{x}$ is broken by the electric field (in the direction of the field), however reversing the direction of the field but at the same time swapping the charges generates a physically identical situation.

The average electrical current, which is given by equation (10), reads

$$ \langle \mathbf{J} \rangle = \sigma_0 \mathbf{E} + q \sum_{\alpha, \beta} \kappa_{\alpha} z_{\alpha} \int i \mathbf{k} \tilde{V}_{\alpha\beta}(-\mathbf{k}) \tilde{C}_{\alpha\beta}(\mathbf{k}) \frac{d\mathbf{k}}{(2\pi)^d}. \quad (25) $$
Recall here that the correction to the bare current is given by the average of the inter-action term which was neglected upon linearizing the full SDFT, thus our expansion is only valid when the computed correction is small. From the expression for the current we can define the field dependent conductivity via

$$ \langle \mathbf{J} \rangle = \sigma(E) \mathbf{E} = [\sigma_0 + \Delta\sigma(E)] \mathbf{E}, \qquad (26) $$

where we assumed that the system is isotropic. We find that the correction to the bare conductivity is given by

$$ \Delta\sigma(E) = -\frac{q^2 \bar{\rho}_+ \bar{\rho}_- (\kappa_+ z_+ - \kappa_- z_-)^2}{(\kappa_+ + \kappa_-) T^2} \times \int \frac{\frac{k_{\parallel}^2}{k^2} \tilde{V}_{+-}^2 \left(1 + \frac{\kappa_+ \bar{\rho}_+ \tilde{V}_{++} + \kappa_- \bar{\rho} \tilde{V}_{--}}{(\kappa_+ + \kappa_-)T}\right)}{\left(1 + \frac{\bar{\rho}_+ \tilde{V}_{++}}{T}\right) \left(1 + \frac{\bar{\rho} \tilde{V}_{--}}{T}\right) \left[\left[1 + \frac{\kappa_+ \bar{\rho}_+ \tilde{V}_{++} + \kappa_- \bar{\rho} \tilde{V}_{--}}{(\kappa_+ + \kappa_-)T}\right]^2 + \left[\frac{(\kappa_+ z_+ - \kappa_- z_-)qE k_{\parallel}}{(\kappa_+ + \kappa_-)T k^2}\right]^2\right]} (2\pi)^d, \\ -\frac{\bar{\rho}_+ \bar{\rho} \tilde{V}_{+-}^2}{T^2} \left(1 + \frac{\kappa_+ \bar{\rho}_+ \tilde{V}_{++} + \kappa_- \bar{\rho} \tilde{V}_{--}}{(\kappa_+ + \kappa_-)T}\right)^2 \qquad (27) $$

where $k_{\parallel}$ denotes the component of the vector $\mathbf{k}$ in the direction of $\mathbf{E}$ (the third line of the equation belongs to the denominator of the integrand).

We notice from equation (27) that the correction $\Delta\sigma(E)$ is zero when $\kappa_+ z_+ = \kappa_- z_-$. In this case, the two ionic types move with the same average velocity $\mathbf{v} = q\kappa_{\pm}z_{\pm}\mathbf{E}$ due to an applied uniform field and the density fluctuation correlation keeps its equilibrium form, which is isotropic and does not modify the average electric field felt by an ion.

At zero field, the correction to the conductivity is given by

$$ \Delta\sigma(0) = -\frac{q^2 \bar{\rho}_+ \bar{\rho}_- (\kappa_+ z_+ - \kappa_- z_-)^2}{d(\kappa_+ + \kappa_-) T^2} \times \int \frac{\tilde{V}_{+-}(k)^2}{\left[1 + \frac{\bar{\rho}_+ \tilde{V}_{++}(k) + \bar{\rho} \tilde{V}_{--}(k)}{T} + \frac{\bar{\rho}_+ \bar{\rho} (\tilde{V}_{++}(k) \tilde{V}_{--}(k) - \tilde{V}_{+-}(k)^2)}{T^2}\right] \left[1 + \frac{\kappa_+ \bar{\rho}_+ \tilde{V}_{++}(k) + \kappa_- \bar{\rho} \tilde{V}_{--}(k)}{(\kappa_+ + \kappa_-)T}\right]} \frac{\mathrm{d}\mathbf{k}}{(2\pi)^d}. \qquad (28) $$

# 4. Zero field conductivity for purely electrostatic interactions

## 4.1. General results for a purely electrostatic interaction

If the interaction is purely electrostatic, i.e. only the first term is present in equation (2), our result at zero field equation (28) reduces to

$$ \Delta\sigma(0) = -\frac{q^6 \bar{\rho}_+ \bar{\rho}_- z_+^2 z_-^2 (\kappa_+ z_+ - \kappa_- z_-)^2}{d(\kappa_+ + \kappa_-) T^2} \int \frac{\tilde{G}_0(\mathbf{k})^2}{\left(1 + \frac{q^2(z_+^2 \bar{\rho}_+ + z_-^2 \bar{\rho}_-) \tilde{G}_0(\mathbf{k})}{T}\right) \left(1 + \frac{\sigma_0 \tilde{G}_0(\mathbf{k})}{(\kappa_+ + \kappa_-)T}\right)} \frac{\mathrm{d}\mathbf{k}}{(2\pi)^d}. \qquad (29) $$
At this point we should note that theories with a $k$ dependent (non-local in space) dielectric function [16] can also be treated within this formalism. We now define, following the notation of Onsager [2],

$$m_{\pm}^2 = \frac{\bar{\rho}_{\pm} z_{\pm}^2 q^2}{\epsilon T} \qquad (30)$$

$$m^2 = m_+^2 + m_-^2 \qquad (31)$$

$$m'^2 = \frac{\kappa_+ m_+^2 + \kappa_- m_-^2}{\kappa_+ + \kappa_-} = \frac{\sigma_0}{\epsilon T (\kappa_+ + \kappa_-)} \qquad (32)$$

Using this, we can write

$$\frac{\Delta\sigma(0)}{\sigma_0} = -\frac{q^2}{d\epsilon T} \frac{m_+^2 m_-^2 (\kappa_+ z_+ - \kappa_- z_-)^2}{m'^2 (\kappa_+ + \kappa_-)^2} \int \frac{1}{([\epsilon \tilde{G}_0(\mathbf{k})]^{-1} + m^2)([\epsilon \tilde{G}_0(\mathbf{k})]^{-1} + m'^2)} \frac{d\mathbf{k}}{(2\pi)^d} \qquad (33)$$

In the case where all the charges are mobile in an electroneutral system, $\sum_{\alpha} z_{\alpha} \bar{\rho}_{\alpha} = 0$ and thus

$$\frac{m_+^2 m_-^2 (\kappa_+ z_+ - \kappa_- z_-)^2}{m'^2 (\kappa_+ + \kappa_-)^2} = -z_+ z_- m'^2, \qquad (34)$$

leading to a simpler form for the correction:

$$\frac{\Delta\sigma(0)}{\sigma_0} = \frac{q^2 z_+ z_- m'^2}{d\epsilon T} \int \frac{1}{([\epsilon \tilde{G}_0(\mathbf{k})]^{-1} + m^2)([\epsilon \tilde{G}_0(\mathbf{k})]^{-1} + m'^2)} \frac{d\mathbf{k}}{(2\pi)^d} \qquad (35)$$

## 4.2. Conductivity for ions moving in homogeneous space

In a homogeneous $d$-dimensional space with solvent dielectric permittivity $\epsilon$, the electrostatic Green function is given by $\tilde{G}_0(\mathbf{k}) = (\epsilon k^2)^{-1}$ and the conductivity correction (equation (33)) becomes

$$\frac{\Delta\sigma(0)}{\sigma_0} = -\frac{q^2}{d\epsilon T} \frac{m_+^2 m_-^2 (\kappa_+ z_+ - \kappa_- z_-)^2}{m'^2 (\kappa_+ + \kappa_-)^2} \int \frac{1}{(k^2 + m^2)(k^2 + m'^2)} \frac{d\mathbf{k}}{(2\pi)^d} \qquad (36)$$

We note that while $m$ is the usual inverse Debye screening length, which only depends on static quantities, the inverse length scale $m'$ depends generally explicitly on the dynamical properties of the system via the mobilities $\kappa_{\alpha}$. However, for a monovalent electrolyte ($z_{\pm} = \pm 1$), $m' = m/\sqrt{2}$ and is thus independent of the ionic mobilities. Furthermore, we note that the correction to the bare conductivity is always negative. Explicit evaluation for dimensions $d=1$ to 3 yields

$$\left( \frac{\Delta\sigma(0)}{\sigma_0} \right)_{d=3} = - \frac{q^2}{12\pi\epsilon T} \frac{m_+^2 m_-^2 (\kappa_+ z_+ - \kappa_- z_-)^2}{m'^2 (\kappa_+ + \kappa_-)^2} \frac{1}{m + m'} \qquad (37)$$
$$ \left( \frac{\Delta\sigma(0)}{\sigma_0} \right)_{d=2} = -\frac{q^2}{4\pi\epsilon T} \frac{m_+^2 m_-^2 (\kappa_+ z_+ - \kappa_- z_-)^2}{m'^2 (\kappa_+ + \kappa_-)^2} \frac{\log\left(\frac{m}{m'}\right)}{m^2 - m'^2}, \quad (38) $$

$$ \left( \frac{\Delta\sigma(0)}{\sigma_0} \right)_{d=1} = -\frac{q^2}{2\epsilon T} \frac{m_+^2 m_-^2 (\kappa_+ z_+ - \kappa_- z_-)^2}{m'^2 (\kappa_+ + \kappa_-)^2} \frac{1}{mm'(m+m')}. \quad (39) $$

For an electroneutral system where all charges are mobile, the correction is given by

$$ \left( \frac{\Delta\sigma}{\sigma_0} \right)_{d=3} = \frac{q^2 z_+ z_-}{12\pi\epsilon T} \frac{m'^2}{m+m'}, \quad (40) $$

$$ \left( \frac{\Delta\sigma}{\sigma_0} \right)_{d=2} = \frac{q^2 z_+ z_-}{4\pi\epsilon T} \frac{\log\left(\frac{m}{m'}\right)}{\left(\frac{m}{m'}\right)^2 - 1}, \quad (41) $$

$$ \left( \frac{\Delta\sigma}{\sigma_0} \right)_{d=1} = \frac{q^2 z_+ z_-}{2\epsilon T} \frac{m'}{m(m+m')}. \quad (42) $$

The result in dimension three agrees with that obtained by Onsager [2].

We now recall that the validity of the linearization of the SDFT can be checked a posteriori by verifying that the correction to the average current or conductivity is small. In all dimensions this condition is satisfied if the electrostatic coupling constant $\Gamma = \frac{q^2 z_+ z_-}{2\pi\epsilon T}$ is small, but the dependence on the density depends on the dimension. For example, for monovalent ions:

* For $d=3$, the correction is proportional to $\sqrt{\bar{\rho}}$ and is small if $\bar{\rho}$ is small.

* For $d=2$, the magnitude of the correction is purely controlled by the coupling constant $\Gamma$, which is dimensionless here, and does not depend on the density. In $d=2$ it is well known that there is a transition from a weak coupling conducting phase to a strong coupling dielectric phase (where the charge carriers exist only in bound pairs)—the Kosterlitz-Thouless transition [17]. The weak coupling approach applied here is clearly only valid in the conducting phase.

* For $d=1$, the correction is proportional to $1/\sqrt{\bar{\rho}}$ and is thus small at high densities, however the long range nature of the electric field generated by one dimensional charges means that the one dimensional problem is not very realistic.

If, instead of the Coulomb interaction, the particles interact with a Yukawa or screened Coulomb interaction, with screening length $\xi$ (i.e. the relevant Green's function becomes $\tilde{G}_0(k) = [\epsilon(k^2 + \xi^{-2})]^{-1}$), in our precedent calculations this simply amounts to replacing the inverse Debye length $m$ by $\sqrt{m^2 + \xi^{-2}}$, and likewise for $m'$, in the integral over $k$ in equation (36).
### 4.3. Conductivity for ions confined to $d=1$ or $2$ and interacting with the 3d Green function

We now turn to the case of charges constrained to move in a low dimensional space ($d=1$ or $2$) but that still interact with the 3-dimensional interaction $G_0(\mathbf{x}) = 1/(4\pi\epsilon|\mathbf{x}|)$. In this case we find that the expressions for the conductivity corrections for systems confined in both two and one dimensions exhibit ultraviolet (large $k$) divergences. These divergences arise due to the singular nature of the three dimensional Coulomb interaction at short distances. The Coulomb interaction between the charges can be regularized by using a Gaussian distribution of charge about each ion type to give a local charge density

$$ \rho_{c\alpha}(\mathbf{x}) = \frac{qz_{\alpha}}{(2\pi\ell_{\alpha}^2)^{d/2}} \exp\left(-\frac{\mathbf{x}^2}{2\ell_{\alpha}^2}\right), \qquad (43) $$

where $\ell_\alpha$ is the size of the region on the ion over which the net charge is localized. In Fourier space this gives an effective electrostatic interaction

$$ \tilde{V}_{\alpha\beta}(\mathbf{k}) = \frac{q^2 z_{\alpha} z_{\beta}}{\epsilon k^2} \exp\left(-\frac{k^2[\ell_{\alpha}^2 + \ell_{\beta}^2]}{2}\right). \qquad (44) $$

#### 4.3.1. Confinement to a plane.
Now if we consider a 3 dimensional interaction between charges restricted to the plane $z=0$, the two dimensional Fourier transform of the effective electrostatic interaction in the $(x, y)$ plane is given by

$$ \tilde{V}_{\alpha\beta}^{d=2}(\mathbf{k}) = \int \tilde{V}_{\alpha\beta}(\mathbf{k}, k_z) \frac{dk_z}{2\pi} \qquad (45) $$

where here $\mathbf{k}$ is a two dimensional vector. This then gives

$$ \tilde{V}_{\alpha\beta}^{d=2}(\mathbf{k}) = \frac{q^2 z_{\alpha} z_{\beta}}{\pi \epsilon k} \int_0^{\infty} \frac{\exp\left(-\frac{[\ell_{\alpha}^2 + \ell_{\beta}^2]k^2[u^2 + 1]}{2}\right)}{u^2 + 1} du. \qquad (46) $$

In two dimensions on taking the limit $\ell_\alpha \to 0$ for both ions we recover the standard result $\tilde{G}_0^{d=2}(\mathbf{k}) = 1/(2\epsilon k)$; clearly this same behavior is recovered in the small $k$ limit and the long distance properties of the Coulomb interaction are thus unaffected by the short distance regularization. For large $k$,

$$ \tilde{V}_{\alpha\beta}^{d=2}(\mathbf{k}) \underset{k \to \infty}{\sim} \frac{q^2 z_{\alpha} z_{\beta}}{\sqrt{2\pi} \epsilon k^2 \sqrt{\ell_{\alpha}^2 + \ell_{\beta}^2}} \exp\left(-\frac{k^2[\ell_{\alpha}^2 + \ell_{\beta}^2]}{2}\right); \qquad (47) $$

subsequently all the formula for the conductivities derived in two dimensions are now rendered finite.

For simplicity if we take $\ell_+ = \ell_- = \ell$ and define $s = \sqrt{2}\ell$, the interaction potential can be written in the scaling form
$$ \tilde{V}_{\alpha\beta}^{d=2}(\mathbf{k}) = \frac{q^2 z_{\alpha} z_{\beta}}{\epsilon} sf(ks). \qquad (48) $$

The resulting integrals for the conductivity correction cannot be carried out analytically but the leading divergence as $s \to 0$ can be extracted by writing $\tilde{V}_{\alpha\beta}^{d=2}(\mathbf{k}) = q^2 z_{\alpha} z_{\beta} / (2\epsilon k)$ for $ks \ll 1$ and cutting off the resulting integral at $k = 1/s$; the regularisation is thus equivalent to the Pauli-Villars regularisation [18]. The result of the regularisation is that the correction to the bare conductivity behaves as

$$ \left( \frac{\Delta\sigma(0)}{\sigma_0} \right)_{d=2, 3\text{d int.}} \sim_{s \to 0} -\frac{q^2 z_+ z_- m'^2 \log(sm^2)}{16\pi\epsilon T}. \qquad (49) $$

The prefactor to the logarithmic correction is small at small coupling constant $\Gamma$ and at small density; the argument of the logarithmic term behaves as $s\Gamma\bar{\rho}$. For charges restricted to a plane we thus see that the change in the conductivity due to interactions behaves as $\bar{\rho}\log(\Gamma\bar{\rho}s)$ as opposed to the $\sqrt{\bar{\rho}}$ correction seen when they are free to move in three dimensions.

### 4.3.2. Confinement to a line.
In one dimension, for charges restricted to the line $y = z = 0$, the effective interaction is

$$ \tilde{V}_{\alpha\beta}^{d=1}(k) = \int \tilde{V}_{\alpha\beta}(k, k_y, k_z) \frac{dk_z dk_y}{(2\pi)^2}, \qquad (50) $$

and thus

$$ \tilde{V}_{\alpha\beta}^{d=1}(k) = \frac{q^2 z_{\alpha} z_{\beta}}{2\pi\epsilon} \int_0^{\infty} \frac{\exp\left(-\frac{[\ell_{\alpha}^2 + \ell_{\beta}^2]k^2[u^2 + 1]}{2}\right)}{u^2 + 1} u du. \qquad (51) $$

The asymptotic behavior of the potential is given by

$$ \tilde{V}_{\alpha\beta}^{d=1}(k) \sim_{k \to 0} -\frac{q^2 z_{\alpha} z_{\beta}}{2\pi\epsilon} \log\left(k\sqrt{\ell_{\alpha}^2 + \ell_{\beta}^2}\right), \qquad (52) $$

$$ \tilde{V}_{\alpha\beta}^{d=1}(k) \sim_{k \to \infty} \frac{1}{2\pi\epsilon k^2 (\ell_{\alpha}^2 + \ell_{\beta}^2)} \exp\left(-\frac{k^2[\ell_{\alpha}^2 + \ell_{\beta}^2]}{2}\right). \qquad (53) $$

Again this regularises all integrals appearing in the formulas for the conductivity both at zero and finite field.

In the simple case where $\ell_+ = \ell_- = s/\sqrt{2}$, the scaling form of $\tilde{V}_{\alpha\beta}^{d=1}(k)$ is given by

$$ \tilde{V}_{\alpha\beta}^{d=1}(k) = \frac{q^2 z_{\alpha} z_{\beta}}{\epsilon} f(ks). \qquad (54) $$

The regularized correction to the bare conductivity thus behaves as
$$ \left( \frac{\Delta\sigma(0)}{\sigma_0} \right)_{d=1, 3d \text{ int.}} \stackrel{s\to 0}{\sim} -\frac{q^2 z_+ z_- m'^2}{\pi \epsilon T s} g(m, m') \quad (55) $$

where $g(m, m')$ is a finite function of $m$ and $m'$ given by

$$ g(m, m') = \int_0^\infty \frac{dp}{[f(p)^{-1} + m^2][f(p)^{-1} + m'^2]} \quad (56) $$

Interestingly, the integral defining $g$ turns out to be convergent when $m = m' = 0$. This means that, for small densities, the change in the conductivity due to interactions is proportional to $\bar{\rho}$.

# 5. Field dependence for monovalent salts

## 5.1. Correlation function

Here we examine the case of monovalent salts (with no background charge, $z_+ = -z_- = 1$ and $\bar{\rho}_+ = \bar{\rho}_- = \bar{\rho}$), which turns out to be the simplest case where one can obtain completely analytical results for the conductivity at any applied field. We define the inverse Debye length $m$ and dimensionless electric field $\mathbf{F}$ by

$$ m^2 = \frac{2\bar{\rho}q^2}{\epsilon T}, \quad (57) $$

$$ \mathbf{F} = \frac{q}{mT} \mathbf{E}. \quad (58) $$

For this case the correlation function is given by

$$ \tilde{C}(\mathbf{k}) = \frac{\bar{\rho}}{\left(1 + \frac{m^2}{2k^2}\right) \left[1 + \frac{m^2}{k^2} + \left(\frac{mk \cdot \mathbf{F}}{k^2}\right)^2\right]} \left( \frac{\left(1 + \frac{m^2}{2k^2}\right)^2 + \left(\frac{mk \cdot \mathbf{F}}{k^2}\right)^2}{2k^2} \left(1 + \frac{m^2}{2k^2} - i \frac{mk \cdot \mathbf{F}}{k^2}\right) \right) \frac{m^2}{2k^2} \left(1 + \frac{m^2}{2k^2} + i \frac{mk \cdot \mathbf{F}}{k^2}\right) \left(1 + \frac{m^2}{2k^2}\right)^2 + \left(\frac{mk \cdot \mathbf{F}}{k^2}\right)^2 \quad (59) $$

The pair correlation functions $h_{++}(\mathbf{x})$ and $h_{-+}(\mathbf{x})$, defined by

$$ h_{\alpha\beta}(\mathbf{x}) = \frac{1}{\bar{\rho}^2} C_{\alpha\beta}(\mathbf{x}) - \frac{\delta_{\alpha\beta}\delta(\mathbf{x})}{\bar{\rho}}, \quad (60) $$

are plotted in figures 1 and 2 for different values of the electric field; they give the change in density of cations and anions, respectively, around a cation. These figures are obtained by inverting the Fourier transforms $\tilde{h}_{\alpha\beta}(k)$. Despite the fact that the expressions for the conductivity that we have found are finite, the integral of $\tilde{h}_{\alpha\beta}(k)$ is divergent at large $k$, meaning that the inverse Fourier transform needs to be regularized (the corresponding integrals at large $k$ decay as $1/k^2$ indicating a $1/|\mathbf{x}|$ divergence at small $|\mathbf{x}|$ in three dimensions). We have thus regularized the resulting integrals using a
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Pair correlation function $h_{++}(\mathbf{x})$ (equation (59)) for different values of the dimensionless electric field $F = qE/(mT)$ (equation (58)); the unit of length is the Debye screening length $m^{-1}$ defined in equation (57). The coordinates $x_\parallel$ and $x_\perp$ denote the directions parallel and perpendicular, respectively, to the field direction. Multiplied by $\bar{\rho}$, the pair correlation function does not depend on $\bar{\rho}$. As explained in the main text, $h_{++}(\mathbf{x})$ is regularized by convolution with a Gaussian whose width is set to the pixel size.</figcaption></figure>

Gaussian cut-off function $\exp(-\ell^2 k^2/2)$, the pair correlation function shown is thus the convolution of the corresponding pair correlation functions with a Gaussian of width $\ell$. In figures 1 and 2, the cut-off length is set to the pixel size, $\ell = 0.2m^{-1}$.

In figure 1, at zero field, $h_{++}(\mathbf{x})$ is negative about the origin and then decays monotonically to zero at large distances. When $F \neq 0$ we see that the correlations have a longer range as screening is destroyed by the field. At larger values of $F$, the probability of an other positive ion being present around a given positive ion is largest in the direction of the field as indicated by the lobes on the figure. In figure 2, without external field, the correlation function is positive about the origin and decays to zero as the distance from the origin increases, indicating screening of positive charges (at the origin) by negative charges. As $F$ is turned on, the negative ions are again more likely to be found in the direction of the applied field but most likely in front of the positive charge than behind, which of course makes physical sense. There is also a region of depletion of negative ions further from the cations, an effect which as also been seen in the pair correlation function in recent simulations of a lattice based ionic model [6] (the correlation function $g_{-+}(\mathbf{x}) = h_{-+}(\mathbf{x}) + 1$ becomes less than 1 in their figure 4). Both figures thus indicate a tendency for the positive and negative ions to form chains in the direction of the field.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Pair correlation function $h_+(x)$ (equation (59)) in monovalent electrolytes for different values of the dimensionless electric field $F = qE/(mT)$ (equation (58)); the unit of length is the Debye screening length $m^{-1}$ (equation (57)). The coordinates $x_\|$ and $x_\perp$ denote the directions parallel and perpendicular, respectively, to the field direction. Multiplied by $\bar{\rho}$, the pair correlation function does not depend on $\bar{\rho}$. As explained in the main text, $h_+(x)$ is regularized by convolution with a Gaussian whose width is set to the pixel size.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Conductivity correction for $d=3$ as a function of the (dimensionless) applied electric field (solid line), and asymptotic expression $\Delta\sigma \sim 1/F$ (dashed).</figcaption></figure>
## 5.2. Conductivity

The field dependent conductivity is given by

$$ \frac{\sigma(F)}{\sigma_0} = 1 - \frac{m^d}{2\bar{\rho}} \int \frac{u_{\parallel}^2}{(2u^2+1)(u^4+u^2+F^2u_{\parallel}^2)} \frac{du}{(2\pi)^d}, \qquad (61) $$

where, again, the subscript ‘‖’ denotes the component parallel to the electric field **E**. For a finite field and a dimension *d* = 3, we need to evaluate the integral

$$ I = \int \frac{u_{\parallel}^2}{(2u^2+1)(u^4+u^2+F^2u_{\parallel}^2)} \frac{du}{(2\pi)^3} = \frac{1}{2\pi^2} \int_0^{\infty} du \int_0^1 dv \frac{u^2 v^2}{(2u^2+1)(u^2+1+F^2v^2)}. \qquad (62) $$

The integral over *u* gives

$$ I = \frac{1}{8\pi} \int_0^1 \frac{2(1 + F^2 v^2) - \sqrt{2} \sqrt{1 + F^2 v^2}}{(1 + 2F^2 v^2) \sqrt{1 + F^2 v^2}} v^2 dv \qquad (63) $$

$$ = \frac{1}{16\pi F^3} \left[ F\sqrt{1+F^2} - \arctan\left(\frac{F}{\sqrt{1+F^2}}\right) - \sqrt{2}F + \arctan(\sqrt{2}F) \right]. \qquad (64) $$

This gives the final result

$$ \frac{\sigma(F)}{\sigma_0} = 1 - \frac{m^3}{32\pi\bar{\rho}F^3} \left[ F\sqrt{1+F^2} - \arctan\left(\frac{F}{\sqrt{1+F^2}}\right) - \sqrt{2}F + \arctan(\sqrt{2}F) \right]. \qquad (65) $$

The first term in brackets is dominant at large fields and the correction decays as $1/F$. This result agrees with that of Onsager and Kim [4]. The expression (65) is plotted on figure 3. The correction decays as the applied field increases, which is the so called Wien effect in simple strong electrolytes.

In lower dimensions, with the corresponding Coulomb interaction in that dimension, the conductivity at finite field is given by

$$ \left(\frac{\sigma(F)}{\sigma_0}\right)_{d=1} = 1 - \frac{m}{4\bar{\rho}} \frac{1}{\sqrt{F^2+1}\left(\sqrt{2(F^2+1)}+1\right)}, \qquad (66) $$

$$ \left(\frac{\sigma(F)}{\sigma_0}\right)_{d=2} = 1 - \frac{m^2}{16\pi\bar{\rho}} \left( \frac{1}{F^2\sqrt{1+2F^2}} \left[ \log(1+F^2+\sqrt{1+2F^2}) - \log(2+3F^2+2\sqrt{(1+F^2)(1+2F^2)}) \right] + \frac{1}{F^2}\log\left(1+\frac{F^2}{2}+\sqrt{1+F^2}\right) \right). \qquad (67) $$
## 6. Conclusion

We have used SDFT to compute the conductivity of strong electrolytes. Our approach is rather simple and compact and reproduces classical results of Onsager and collaborators which were derived by very different methods. We have furthermore seen how the calculations of Onsager are modified when there is a static uniform background charge. In this situation, even a system where all ions have the same charge but have different mobilities exhibits a negative correction to the bare conductivity. Our formalism has also allowed us to generalize the results of Onsager to Coulomb systems in arbitrary dimensions.

Clearly our results are applicable to more general systems; they can be applied to more complex electrostatic models, for instance models with non-local dielectric constants [16]. The general formulas given here can also be applied to the motion of ions confined to quasi-1d geometries, e.g. in carbon or boron-nitride nanotubes [19–21]. Finally, our formalism can treat general non-electrostatic interactions like steric interactions, which enters in the dynamics of ionic liquids [22].

Further extensions of this work are clearly possible, one could for example analyze the influence of a solvent on the conductivity. The effect of a solvent in Onsager's calculations is in fact additive at the level of approximation he employed, the additional term is generated by the solvent induced flow due to the movement of the ions. Using the formalism here it is possible that the perturbation analysis could be carried out in a self consistent manner thus potentially improving its results at higher concentrations in three dimensions.

## References

[1] Debye P and Hückel E 1923 *Phys. Z.* **24** 305–25

[2] Onsager L 1927 *Trans. Faraday Soc.* **23** 341–9

[3] Hemmer P C, Holden H and Kjelstrup Ratkje S 1996 *The Collected Works of Lars Onsager: with Comments vol 17* (Singapore: World Scientific)

[4] Onsager L and Kim S K 1957 *J. Phys. Chem.* **61** 215–29

[5] Onsager L 1934 *J. Chem. Phys.* **2** 599–615

[6] Kaiser V, Bramwell S T, Holdsworth P C W and Moessner R 2013 *Nat. Mater.* **12** 1033–7

[7] Kawasaki K 1994 *Phys. A: Stat. Mech. Appl.* **208** 35–64

[8] Dean D S 1996 *J. Phys. A: Math. Gen.* **29** L613–7

[9] Velenich A, Chamon C, Cugliandolo L F and Kreimer D 2008 *J. Phys. A: Math. Theor.* **41** 235002

[10] Kim B, Kawasaki K, Jacquin H and van Wijland F 2014 *Phys. Rev. E* **89** 012150

[11] Dean D S and Podgornik R 2014 *Phys. Rev. E* **89** 032117

[12] Démery V, Bénichou O and Jacquin H 2014 *New J. Phys.* **16** 053032

[13] Dean D S and Lefèvre A 2004 *Phys. Rev. E* **69** 061111

[14] Bun J, Bouchaud J P, Majumdar S N and Potters M 2014 *Phys. Rev. Lett.* **113** 070201

[15] Zwanzig R 2001 *Nonequilibrium Statistical Mechanics* (Oxford: Oxford University Press)

[16] Bopp P A, Kornyshev A A and Sutmann G 1996 *Phys. Rev. Lett.* **76** 1280–3

[17] Kosterlitz J M and Thouless D J 1973 *J. Phys. C: Solid State Phys.* **6** 1181

[18] Itzykson C and Zuber J B 2006 *Quantum Field Theory* (New York: Dover) (Dover books on physics)

[19] Liu H, He J, Tang J, Liu H, Pang P, Cao D, Krstic P, Joseph S, Lindsay S and Nuckolls C 2010 *Science* **327** 64–7

[20] Siria A, Poncharal P, Biance A L, Fulcrand R, Blase X, Purcell S T and Bocquet L 2013 *Nature* **494** 455–8

[21] Secchi E, Nigues A, Jubin L, Siria A and Bocquet L 2015 arXiv e-prints (preprint 1509.01026)

[22] Lee A A, Kondrat S, Vella D and Goriely A 2015 *Phys. Rev. Lett.* **115** 106101