Fluctuations of absorption of interacting diffusing particles by multiple absorbers

TL;DR

This paper analyzes the fluctuations in particle absorption by multiple patches on a domain with interacting diffusing particles, revealing how correlations depend on transport coefficients and linking absorption statistics to other physical systems.

Contribution

It extends previous work by calculating the covariance matrix of absorption fluctuations for multiple patches using macroscopic fluctuation theory, and identifies conditions for correlation signs independent of geometry.

Findings

Correlation sign determined by transport coefficients

Flux field develops vorticity during fluctuations

Absorption statistics linked to non-equilibrium current statistics

Abstract

We study fluctuations of particle absorption by a three-dimensional domain with multiple absorbing patches. The domain is in contact with a gas of interacting diffusing particles. This problem is motivated by living cell sensing via multiple receptors distributed over the cell surface. Employing the macroscopic fluctuation theory, we calculate the covariance matrix of the particle absorption by different patches, extending previous works which addressed fluctuations of a single current. We find a condition when the sign of correlations between different patches is fully determined by the transport coefficients of the gas and is independent of the problem's geometry. We show that the fluctuating particle flux field typically develops vorticity. We establish a simple connection between the statistics of particle absorption by all the patches combined and the statistics of current in a non-equilibrium steady state in one dimension. We also discuss connections between the absorption statistics and (i) statistics of electric currents in multi-terminal diffusive conductors and (ii) statistics of wave transmission through disordered media with multiple absorbers.

Equations

$${\partial }_t\rho =\nabla \cdot \left[D(\rho )\nabla \rho \right],$$

$${\partial }_t\rho =-\nabla \cdot \mathbf{J},\mathbf{J}=-D(\rho )\nabla \rho -\sqrt{\sigma (\rho )}\mathbf{\eta }(\mathbf{x},t),$$

$$\mathcal{P}\simeq \frac{T^{s/2}}{(2\pi {)}^{s/2}[\text{det}\mathbf{C}{]}^{1/2}}\exp \left(-\frac{T}{2}\sum _{i,j=1}^s\delta n_i{C}_{ij}^{-1}\delta n_j\right),$$

$$\rho (\mathbf{x}\in {\Omega }_i,t)=0,i=1,2,\dots ,s$$

$$\mathbf{J}(\mathbf{x}\in {\Omega }_r,t)\cdot \hat{n}=0$$

$$\rho (\mathbf{x},t=0)={\rho }_0.$$

$$\rho (|\mathbf{x}|\to \infty ,t)={\rho }_0.$$

$${\partial }_t\rho$$

$${\partial }_{t}p$$

$$\mathbf{J}=-D(\rho )\nabla \rho +\sigma (\rho )\nabla p.$$

$$\rho (\mathbf{x},t=0)$$

$$p(\mathbf{x},t=T)$$

$$\rho (|\mathbf{x}|\to \infty ,t)$$

$$p(|\mathbf{x}|\to \infty ,t)$$

$$\rho (\mathbf{x}\in {\Omega }_i,t)$$

$$p(\mathbf{x}\in {\Omega }_i,t)$$

$$\nabla \rho (\mathbf{x}\in {\Omega }_r,t)\cdot \hat{n}$$

$$-$$

$$\nabla \cdot \left[D(\bar{\rho })\nabla \bar{\rho }\right]=0$$

$$\varphi \left[\bar{\rho }\left(\mathbf{x}\right)\right]={\int }_{\bar{\rho }(\mathbf{x})}^{{\rho }_0}D(w)dw.$$

$${\nabla }^2\varphi =0.$$

$$\varphi (\mathbf{x}\in {\Omega }_i)=V({\rho }_0)\equiv {\int }_0^{{\rho }_0}D(w)dw,$$

$$\nabla \varphi (\mathbf{x}\in {\Omega }_r)\cdot \hat{n}=0,\varphi (\mathbf{x}\to \infty )=0.$$

$$\varphi \left(\bar{\rho }\right)=D_0\left({\rho }_0-\bar{\rho }\right)$$

$$\bar{\mathbf{J}}=-D(\bar{\rho })\nabla \bar{\rho }=\nabla \varphi .$$

$${\bar{n}}_i={\oint }_{{\Omega }_i}\bar{\mathbf{J}}\cdot \hat{n}dS={\oint }_{{\Omega }_i}\nabla \varphi \cdot \hat{n}dS,i=1,2,\dots ,s.$$

$$A_{ij}={\oint }_{{\Omega }_i}\nabla {\varphi }_j\cdot \hat{n}dS.$$

$$\varphi (\mathbf{x})=V({\rho }_0)\sum _{i=1}^s{\varphi }_i(\mathbf{x}).$$

$${\bar{n}}_i=V({\rho }_0)\sum _{j=1}^s{A}_{ij}.$$

$$A_{ij}={\oint }_{\Omega }{\varphi }_i\nabla {\varphi }_j\cdot \hat{n}dS$$

$$\bar{n}=\sum _{i=1}^s{\bar{n}}_i=V({\rho }_0)A,$$

$$A\equiv \sum _{i,j=1}^s{A}_{ij}$$

$$\nabla \cdot \mathbf{\delta J}$$

Full text

I Introduction

Fluctuations of currents of matter and energy is an important subject of nonequilibrium statistical mechanics. A prototypical model problem, which has attracted much attention, involves a diffusive lattice gas driven by two reservoirs of particles kept at different densities lebowitz ; Derrida2007 ; shpiel ; bd ; MFTreview . This simple setting has a direct relevance to experiment in at least two different contexts: statistics of electric current in mesoscopic conductors jor ; meso ; multi and statistics of wave transmission through disordered media penini ; shapiro2 . Here we consider a different but closely related problem: transport of diffusing molecules into the living cell through receptors distributed on its surface berg ; bergbook . The surrounding gas serves as a finite-density reservoir, whereas the cell receptors can be modeled as reservoirs kept at zero gas density. In their pioneering 1977 paper, Berg and Purcell berg evaluated the expected steady-state current of non-interacting diffusing particles into a single receptor. Their motivation was to assess physical limitations of the cell’s ability to sense changes in the environmental concentration berg ; bergbook . Building on their work, here we aim at (i) evaluating the current fluctuations at each receptor during a given time, and (ii) accounting for inter-particle interactions in the surrounding gas. These extensions are presently possible due to recent advances in the fluctuating hydrodynamics Spohn and the macroscopic fluctuation theory (MFT) MFTreview of diffusive lattice gases. These formalisms have already been successfully used in the simplest two-reservoir setting Derrida2007 ; MFTreview ; shpiel ; bd , and in several non-stationary settings which involved a single current. As we will show here, the presence of multiple absorbing patches, leading to multiple currents, allows one to ask new questions, and brings new effects. One new question concerns the joint probability distribution of, and correlations between, the absorption currents at different receptors. One new effect is that the most probable particle flux field, conditioned on a specified joint absorption statistics, exhibits a large-scale vorticity.

We will model the surrounding gas of interacting particles as a diffusive lattice gas Spohn ; Liggett . The large-scale long-time behavior of such gases can be described by fluctuating hydrodynamics Spohn ; KL . The average particle density $\rho(\mathbf{x},t)$ of a lattice gas obeys a diffusion equation

[TABLE]

whereas macroscopic fluctuations are described by the conservative Langevin equation Spohn ; KL

[TABLE]

where $\boldsymbol{\eta}(\mathbf{x},t)$ is a zero-mean Gaussian noise, delta-correlated in space and in time. As one can see from Eq. (2), a diffusive lattice gas is completely specified by two transport coefficients: the diffusivity $D(\rho)\geq 0$ and the mobility $\sigma(\rho)\geq 0$. For lattice gases $\sigma(0)=0$.

The simplest example of a diffusive lattice gas is a gas of non-interacting Random Walkers (RWs). For the RWs one has $D(\rho)=D_{0}=\text{const}$, and $\sigma(\rho)=2D_{0}\rho$ Spohn . An example of interacting lattice gas is the Simple Symmetric Exclusion Process (SSEP). Here at each time step a particle can jump, with equal probability, to any empty neighboring site. The average behaviors of the RWs and the SSEP turn out to be identical: they share the same density-independent diffusivity $D_{0}$. The inter-particle interactions of the SSEP are manifested at the level of fluctuations, as the SSEP’s mobility $\sigma(\rho)=2D_{0}\rho(1-\rho)$ is a non-linear function of $\rho$ Spohn . The lattice gases are not the only systems describable by the Langevin equation (2). Important additional examples describe transport of noninteracting electrons in mesoscopic materials at zero temperature jor and wave transmission in disordered media penini ; shapiro2 .

Now we formulate the model which we will study in this paper. Consider a lattice gas of initially uniform density $\rho(\mathbf{x},t=0)=\rho_{0}$ which fills the whole space outside of a simple-connected three-dimensional domain (the cell) of the characteristic linear size $L$. The domain boundary $\Omega$ includes absorbing patches (the cell receptors) $\Omega_{i},\,i=1,2,\dots,s$. Whenever a particle hits any of these, it is immediately absorbed. Whenever a particle hits the rest of the surface, $\Omega_{r}$ (the cell wall), it is reflected, see Fig. 1.

At times $T$ much longer than the characteristic diffusion time $L^{2}/D(\rho_{0})$, the system reaches a non-equilibrium steady state. In the steady state the average gas density $\bar{\rho}(\mathbf{x})$ and the average absorption current $\bar{n}_{i}=\bar{N}_{i}/T$ into each absorbing patch are independent of time. Here $\bar{N}_{i}$ is the average number of particles absorbed by the $i$-th patch during the time interval $0<t<T$. In this work we will determine the joint probability distribution, ${\mathcal{P}}(\delta n_{1},\delta n_{2},\dots,\delta n_{s};\rho_{0},T)$, of observing small fluctuations, $\delta n_{i}=(N_{i}-\bar{N}_{i})/T$, of the absorption currents during the time $T$. As we show here, this multi-variate probability distribution is Gaussian and given by the expression

[TABLE]

where $\boldsymbol{C}$ is an $s\times s$ positive-definite symmetric matrix which depends on $\rho_{0}$ and on the geometry of the problem, but is independent of time. We obtain a general condition when the sign of the currents’ cross-correlation $\overline{\delta n_{i}\delta n_{j}}$ is independent of the geometry and completely specified by the transport coefficients $D(\rho)$ and $\sigma(\rho)$ of the gas.

Further, we show that the cross-correlation between the current into a single absorbing patch and the total current into all patches combined has a simple structure where the dependence on the problem’s geometry is factorized out. We apply this result to the Berg-Purcell model of a living cell berg . An important finding of Ref. berg was the dependence of the expected total current on the number of receptors distributed on the cell’s surface. We extend their result by finding how the number of receptors affects correlations, and also account for interactions.

Our calculations employ the MFT in conjunction with the additivity principle. The latter was proposed by Bodineau and Derrida bd in the context of statistics of current in a one-dimensional lattice gas driven by two boundaries. We determine the optimal (most probable) spatial profiles of the fluctuating gas density and flux fields, conditioned on the specified absorption current into each patch. By virtue of the additivity principle the optimal absorption current into the $i$-th patch is independent of time and equal to $n_{i}=N_{i}/T$. In order to calculate the variance of the probability distribution (3), we use the approach of Ref. KrMe and linearize the MFT equations around the deterministic solution. The additivity principle and linearization enable us to solve the problem in quite a general form and for an arbitrary diffusive lattice gas. We show that the optimal flux field, conditioned on the specified absorption current into each patch, exhibits a large-scale vortex structure. This feature is unique to multi-reservoir systems sustaining multiple currents, and it appears even when the surrounding gas is modeled as non-interacting RWs. The vorticity is absent in systems with a single current shpiel . Remarkably, it is also absent when the process is conditioned on the total absorption current into all patches combined.

Although the original motivation for this work came from a living cell sensing via multiple receptors, many of our results can be generalized to other diffusive systems which are driven by multiple reservoirs and therefore sustain multiple currents. The reservoirs can be disjoint (rather than placed on a common reflecting boundary), and the system can be finite. Still, for concreteness we focus here on the setting shown in Fig. 1.

Here is how the remainder of the paper is organized. In Sec. II we present the MFT formulation of the problem. The deterministic limit is discussed in Sec. III. The joint distribution of fluctuating absorption currents is obtained in Sec. IV.1 and analyzed in Sec. IV.2. In Sec. IV.3 we briefly discuss the shot-noise-driven fluctuations of current in multi-terminal diffusive conductors, previously studied in Ref. multi . In Sec. IV.4 we determine the statistics of total absorption current into all patches combined, and relate it to previous findings for single-current settings. We also find, in Sec. IV.4.1, the cross-correlation between the total current and the current into a single absorbing patch and apply this finding in Sec. IV.4.2 to the Berg-Purcell model. In Sec. V we study optimal fluctuations of the density field and of the flux field and uncover a large-scale vortex structure of the optimal flux field. Section V also discusses the particular case of the density and flux fields conditioned on the total absorption current. We discuss our main results in Sec. VI. Some of the technical derivations are relegated to appendixes.

II Macroscopic fluctuation theory of joint absorption statistics

The starting point of the MFT formulation is the Langevin equation (2) with the boundary conditions

[TABLE]

at the absorbing patches, and

[TABLE]

at the reflecting surface. The fluctuating flux field $\mathbf{J}$ is defined in Eq. (2). Here and in the following $\hat{n}$ denotes a local unit vector normal to the domain boundary and directed into the domain. At $t=0$ the gas has a uniform density $\rho_{0}$,

[TABLE]

The boundary condition at infinity is, therefore,

[TABLE]

In Appendix A we derive the MFT equations for the joint absorption statistics. The derivation, by now pretty standard MFTreview , yields the governing equations [Eqs. (8) and (9) below] and problem-specific boundary conditions. We repeat the derivation here in order to establish the previously unknown boundary conditions on the absorbing patches. The derivation starts from a path-integral formulation for Eq. (2) with specified numbers of absorbed particles $N_{i}$ by time $T$. The derivation exploits a large parameter – the typical number of particles in relevant regions of space – to perform a saddle-point evaluation of the path integral. The ensuing minimization procedure yields the Euler-Lagrange equations which can be cast into a Hamiltonian form for the optimal density history $\rho(\mathbf{x},t)$ (where $\rho$ plays the role of a “coordinate”) and the conjugate momentum density $p(\mathbf{x},t)$:

[TABLE]

The Hamiltonian $H$ is given in Eq. (82), and the prime denotes the derivative with respect to the single argument. The optimal flux field is given by:

[TABLE]

The boundary conditions in time are

[TABLE]

where here and in the following $\mathbf{x}$ is outside of the domain. The boundary conditions in space are the following. Far from the domain the gas is unperturbed, so

[TABLE]

On the domain boundary we have

[TABLE]

where $\lambda_{i}$ are a priori unknown Lagrange multipliers which are ultimately set by the $s$ constraints of having the specified numbers $N_{i}$ of particles absorbed by the patches by time $T$. The boundary conditions (15) for $p$ generalize their simple analogs in single-current settings bd ; MFTreview ; shpiel ; main ; b ; fullabsorb .

Having solved the coupled nonlinear partial differential equations (8) and (9) with the boundary conditions in space and time, one determines the optimal history of the system conditioned on $N_{i}$, $i=1,2,\dots,s$. With the solutions at hand, one can calculate the action $S$ which yields $-\ln{\mathcal{P}}$ up to a pre-exponential factor:

[TABLE]

Here and in the following the volume integral $\int d\mathbf{x}$ is performed over all space outside of the domain. For typical fluctuations the pre-exponential factor in $\mathcal{P}$, see Eq. (3), is determined from normalization to unity.

III Deterministic theory

The choice $\lambda_{i}=0$ sets $p$ to vanish at all times and describes the deterministic solution, where all $N_{i}$ are equal to their expected values $\bar{N}_{i}$. In this case Eq. (8) reduces to the deterministic diffusion equation (1) for the average density. At long times, $T\gg L^{2}/D(\rho_{0})$, the solution approaches a stationary one, $\bar{\rho}(\mathbf{x})$, obeying the time-independent equation

[TABLE]

and the boundary conditions in space. As a result, a steady-state particle flux field, and steady-state currents into the absorbing patches, set in. Essentially, this is the problem which Berg and Purcell berg solved for gases with constant diffusivity $D_{0}$. They did it using an illuminating electrostatic analogy shpitzer ; Rednerbook which can be easily generalized to a density-dependent diffusivity. Let us introduce the following function of the steady-state density $\bar{\rho}(\mathbf{x})$:

[TABLE]

Then Eq. (18) for $\bar{\rho}(\mathbf{x})$ becomes the Laplace’s equation for $\phi$:

[TABLE]

The boundary conditions for $\bar{\rho}(\mathbf{x})$ transform to the following boundary conditions for $\phi$:

[TABLE]

As in Ref. berg , $\phi$ can be interpreted as the electrostatic potential outside an insulating domain with boundary $\Omega$ over which $s$ conducting patches $\Omega_{i}$, held at voltage $V({\rho_{0}})$, are distributed. Having solved for this potential, one obtains the complete solution of the deterministic problem: $\bar{\rho}$ is obtained by inverting the relation for $\phi\left(\bar{\rho}\right)$ in Eq. (19). For gases with constant diffusivity $D_{0}$, such as the RWs and SSEP, Eq. (19) defines a linear relation

[TABLE]

which we will use in the following. The deterministic steady-state flux field $\bar{\mathbf{J}}$ is minus the electric field of this system:

[TABLE]

In their turn, the average steady-state currents $\bar{n}_{i}=\bar{N}_{i}/T$ are (up to a factor of $4\pi$) the electric charges accumulated on the conducting patches:

[TABLE]

These charges are linearly related to the voltage $V(\rho_{0})$ via an $s\times s$ symmetric capacitance matrix (which we will denote by $\mathbf{A}$), determined solely by the problem’s geometry. The $i,j$ element of $\mathbf{A}$ is the charge on the patch $i$ induced by the unit voltage applied to the patch $j$, the rest of the patches being grounded smythe . $\mathbf{A}$ can be expressed via a set of $s$ characteristic electrostatic potentials $\phi_{i}(\mathbf{x})$. Each of them appears when the corresponding conducting patch $\Omega_{i}$ is held at unit voltage, the rest of the conducting patches are grounded, and the Neumann boundary condition is specified at the reflecting part of the boundary $\Omega_{r}$. The capacitance matrix is given by

[TABLE]

In its turn, $\phi(\mathbf{x})$ from Eq. (19) can be expressed as

[TABLE]

Using Eqs. (25)-(27), we obtain

[TABLE]

To highlight the symmetry of the capacitance matrix $\mathbf{A}$, let us rewrite Eq. (26) as

[TABLE]

where we used the boundary conditions for $\phi_{i}$, the Gauss theorem, and the fact that $\phi_{i}$ are harmonic functions in the bulk. We will call the integrand of the last expression in Eq. (29), $\nabla\phi_{i}\cdot\nabla\phi_{j}$, the capacitance density of the system. Finally, the total absorption current into all patches in this interpretation is equal to the total charge on all the patches when they are held at voltage $V(\rho_{0})$:

[TABLE]

where

[TABLE]

is the total capacitance of the system, that is the total charge accumulated on all the patches when they are kept at unit voltage.

IV Absorption Statistics

IV.1 Solving linearized MFT equations

As already mentioned, we will solve the MFT problem under two simplifying assumptions. The first is the additivity principle bd : Being interested in the long-time limit, $T\gg L^{2}/D(\rho_{0})$, we look for stationary solutions of Eqs. (8) and (9) which satisfy the boundary conditions (12)-(16) blayers .

The second assumption involves linearization of the stationary MFT equations around the steady-state deterministic solution $\rho=\bar{\rho}(\mathbf{x})$ and $p=0$. This corresponds to typical, small fluctuations and suffices for the evaluation of the variance of the joint probability distribution KrMe . Note that, for the typical fluctuations, the additivity principle appears to be a safe assumption for all diffusive lattice gases Bertini2005 ; BD2005 ; Hurtado ; ZM2016 .

Let us denote small deviations of $\rho$ and $p$ from their average values as $\rho_{1}(\mathbf{x})=\rho(\mathbf{x})-\bar{\rho}(\mathbf{x})$ and $p_{1}(\mathbf{x})=p(\mathbf{x})$. They determine time-independent current deviations $\delta n_{i}=n_{i}-\bar{n}_{i}$. The linearized steady-state MFT equations read:

[TABLE]

The boundary conditions for $\rho_{1}$ and $p_{1}$ follow from Eqs. (12)-(16):

[TABLE]

As one can see, the Laplace’s equation (33) for $p_{1}$ is decoupled from Eq. (32). Subject to the boundary conditions (35)-(37), it has a unique solution which can be expressed in terms of the auxiliary potentials $\phi_{i}(\mathbf{x})$, introduced in the previous section:

[TABLE]

This solution suffices for determining the action (17) in terms of $\lambda_{i}$-s. Indeed, in the leading order Eq. (17) yields

[TABLE]

Plugging Eq. (38) into Eq. (39), we obtain the action in terms of a bilinear form in $\lambda_{i}$-s:

[TABLE]

where $\boldsymbol{\Lambda}$ is a vector with components $\lambda_{i}$, and $\boldsymbol{C}$ is a $s\times s$ symmetric matrix,

[TABLE]

given in terms of the volume integral of the product of the capacitance density of the system $\nabla\phi_{i}\cdot\nabla\phi_{j}$ and the effective local noise magnitude $\sigma[\bar{\rho}(\mathbf{x})]$. As we will see shortly, this is the covariance matrix. It is fully determined by the deterministic solution, and it plays a crucial role in our results.

What is left is to express the $\lambda_{i}$-s in Eq. (40) in terms of the current deviations $\delta n_{i}$. This requires solving Eq. (32) which, once $p_{1}$ is known, is a Poisson’s equation for $\rho_{1}$:

[TABLE]

In Appendix F we present an explicit solution of this equation for the RWs. For a general lattice gas an explicit solution is unavailable. Still, we were able to derive an explicit relation for $\delta n_{i}$ vs. $\lambda_{i}$ by using a Green’s function identity, see Appendix B. This relation is given by the same matrix $\boldsymbol{C}$ defined in Eq. (41):

[TABLE]

where $\boldsymbol{\delta}\mathbf{n}$ is the vector with components $\delta n_{i}$. As shown in Appendix C, the symmetric matrix $\boldsymbol{C}$ is positive definite and therefore invertible, enabling one to solve Eq. (43) for $\boldsymbol{\Lambda}$. Plugging this inverse relation in the bilinear form (40) and using Eq. (39), we obtain

[TABLE]

This probability distribution describes Gaussian fluctuations. Normalizing it to unity, we arrive at the announced result (3).

IV.2 Covariance matrix

As is clear from Eq. (44), the joint statistics of absorption is encoded in the covariance matrix $\boldsymbol{C}$. The diagonal elements of $\boldsymbol{C}$ describe the variance of the current into the patch $i$:

[TABLE]

where the over-line denote averaging with respect to the Gaussian distribution (3). The off-diagonal elements of $\boldsymbol{C}$ describe cross-correlations between the currents into different patches:

[TABLE]

A necessary condition for nonzero cross-correlations is inter-particle interactions. This can be seen from an alternative expression for $\boldsymbol{C}$ which we derived in Appendix D:

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta. For the non-interacting RWs one has $\sigma^{\prime}(\rho)/D(\rho)=2$, so the integrand vanishes, and one is left with $C_{ij}=\bar{n}_{i}\delta_{i,j}$, no correlations. From Eq. (45) one can see that, for the RWs, the variance of the number of absorbed particles in each patch, $\overline{\delta N_{i}^{2}}$, is equal to the mean value $\overline{N}_{i}$, regardless of the system’s geometry. In fact, for the RWs the steady-state MFT problem can be solved exactly, without linearization. We performed these calculations and found the complete joint distribution of the number of absorbed particles. As to be expected, this distribution is equal to a product of independent Poisson distributions with the mean $\overline{N}_{i}$.

A simple interacting lattice gas with zero cross-correlations of absorption by different patches is the zero-range process, where a particle can hop to a neighboring site with a rate which increases with the number of particles on the departure site but is independent of the number of particle on the target site zrp . For the zero-range process one has, in the hydrodynamic limit, $D(\rho)=(1/2)\,\sigma^{\prime}(\rho)$ void , and the integrand in Eq. (47) vanishes.

Going back to general $D(\rho)$ and $\sigma(\rho)$, we notice that the expression $\phi_{i}\phi_{j}\left(\nabla\phi\right)^{2}/D(\bar{\rho})$ under the integral in Eq. (47) is everywhere positive decay . Therefore, if

[TABLE]

for any value of $\bar{\rho}\in[0,\rho_{0}]$, then the currents into different patches $i\neq j$ are all anti-correlated, $\overline{\delta n_{i}\delta n_{j}}<0$, regardless of the system’s geometry. In particular, this is true for the SSEP, where $\left[\sigma^{\prime}\left(\rho\right)/D(\rho)\right]^{\prime}=-4$.

If the inequality in Eq. (48) is reversed, the absorption currents are all positively correlated. This happens for the Kipnis-Marchioro-Presutti (KMP) model KMP ; MFTreview which describes an ensemble of agents on a lattice which randomly redistribute energy among neighbors. Using the KMP model one can study fluctuations of energy absorption by absorbing patches located on the domain boundary. For the KMP model $D(\rho)=D_{0}=\text{const}$ and $\sigma\left(\rho\right)=2D_{0}\rho^{2}$, and so $\left[\sigma^{\prime}\left(\rho\right)/D(\rho)\right]^{\prime}=4$.

Remarkably, the same inequality (48) guarantees the validity of the additivity principle for arbitrary currents ber , and also determines the sign of the two-point density correlation function bertinistat ; tridib , in single-current systems.

IV.3 A comparison with Sukhorukov and Loss multi

Sukhorukov and Loss multi studied shot-noise-driven fluctuations of current in multi-terminal diffusive conductors. Although theirs and our geometries are different, their expression (3.16) for the zero-frequency mode of the power spectrum of correlations of the current has a mathematical structure which resembles that of our covariance matrix (41). In their case it is a volume integral over the “conductance density” – an analog of our capacitance density – times the local noise magnitude. If we closely examine the governing equations of both systems, this resemblance should not come as a surprise. Sukhorukov and Loss started from a Boltzmann-Langevin description of the distribution function of electrons (where noise comes from electron scattering on static impurities). Then, applying a diffusion approximation, they derived a time-independent linear Langevin equation for the electrical potential $V(\mathbf{x})$ inside the conductor: Eqs. (2.12) and (2.13) of Ref. multi . For an isotropic medium their equation can be written as

[TABLE]

where $\boldsymbol{\eta}$ is a Gaussian noise, delta-correlated in $\mathbf{x}$ and $t$. The transport coefficients $D(\mathbf{x})$ and $\Pi\left(\mathbf{x}\right)$ depend on $\mathbf{x}$ but are independent of the fluctuating potential $V(\mathbf{x})$ and of time. Equation (49), therefore, is very different from the nonlinear time-dependent Langevin equation (2) which was the starting point of our analysis. However, in order to make a progress within the MFT formalism, we employed the additivity principle and linearization. In fact, we could have derived the same final results from the following Langevin equation which describes small quasi-stationary density fluctuations $\rho_{1}(\mathbf{x})$ around the steady-state average density profile $\bar{\rho}(\mathbf{x})$:

[TABLE]

The mathematical difference between the linear Langevin equations (49) and (50) comes from the fact that the transport coefficients $D(\rho)$ and $\sigma(\rho)$ in Eq. (50) are $\rho$-dependent. In particular, the $\rho$-dependence of $D$ causes an additional contribution to the fluctuating flux field, $D^{\prime}\left(\bar{\rho}\right)\rho_{1}\nabla\bar{\rho}$, which is absent in Eq. (49). Of course, this difference reflects different physical problems which Ref. multi and this work address.

Notwithstanding these differences, there is a close similarity between the two problems. An extensively discussed regime in Ref. multi is that of zero-temperature elastic scattering, manifested in a particular form of the spatial mobility profile $\Pi(\mathbf{x})$ in Eq. (49). An important finding of Ref. multi for this case is the strictly negative sign of the cross-correlations of the current, regardless of the system’s geometry. Remarkably, the profile $\Pi(\mathbf{x})$ for such a system turns out to be equivalent, up to irrelevant factors, to the mobility coefficient $\sigma\left[\bar{\rho}(\mathbf{x})\right]$ for the SSEP. As we saw in the previous section, the cross-correlations of the absorption currents for the SSEP are indeed strictly negative.

Here we deal with arbitrary $D(\rho)$ and $\sigma(\rho)$, where the currents’ cross-correlations can behave differently. In particular, the sign of the cross-correlations can be strictly positive, as it happens for the KMP model. In addition, the MFT formalism, which we employ here, makes it possible to determine the optimal fluctuating profiles, as we report in Sec. V.

IV.4 Statistics of total absorption current

Given the joint absorption statistics, Eq. (3) or (44), what is the probability density that the total current into all the patches combined, $\sum_{i=1}^{s}{n_{i}}$, is equal to a prescribed value $n$? To answer this question, we can minimize the action (44) under the constraint $\sum_{i=1}^{s}n_{i}=n$ or, equivalently, $\sum_{i=1}^{s}\delta n_{i}=n-\bar{n}$. The minimization, performed in Appendix E, yields the following optimal values of the individual currents $n_{i}$:

[TABLE]

The corresponding Lagrange multipliers $\lambda_{i}$-s [see Eq. (15)] turn out to be all equal:

[TABLE]

where

[TABLE]

Note that $I_{1}(\rho_{0})=V(\rho_{0})$, see Eq. (21). Now we can compute the action by using either Eqs. (44) and (51) or Eqs. (40) and (52). An explicit result can be obtained with the help of the relation

[TABLE]

derived in Appendix E. In this way we arrive at a Gaussian distribution of the total absorption current $n$:

[TABLE]

The first two distribution cumulants of the number of absorbed particles $N=Tn$, following from Eq. (55), are

[TABLE]

where we have used Eq. (30). Remarkably, these cumulants are equal to the total capacitance $A$ multiplied by the cumulants of the integrated current, obtained for a one-dimensional lattice gas driven by two reservoirs, at $\rho_{a}=\rho_{0}$ and $\rho_{b}=0$, see Eq. (3) of Ref. bd .

Similar relations were established, for all cumulants of the current, in Ref. shpiel which dealt with the SSEP driven by two reservoirs in a finite system in any spatial dimension. There too the cumulants of the current are equal to the corresponding cumulants of the one-dimensional system multiplied by a constant geometrical factor shpiel . As was noticed in Ref. MFTreview , this geometric factor is the electric capacitance.

Our setting involves an infinite system with a single reservoir at infinity. Still, as we have just shown, the first two cumulants of the total absorption current have the same structure as those in the finite two-reservoir settings. (For a spherically symmetric absorber this property was established, for the complete absorption statistics, in Ref. fullabsorb .) Extending the results of Refs. shpiel and fullabsorb , the relations (56) show that the dependence on the problem’s geometry, through the electric capacitance $A$, is factorized out from the expressions for the current’s mean and variance. We now show that the same feature also holds for the cross-correlations.

IV.4.1 Cross-correlation between the total current and the current into a single absorbing patch

This cross-correlation is obtained by summing over $j$ in Eq. (46) and using Eqs. (54), and (28):

[TABLE]

Note that the variance in (56) can be obtained from the cross-correlation by summing over $i$ and using the definition (31). We see that the dependence on the geometry in Eq. (57) is factorized out as in Eq. (56), although via a different geometrical factor. This geometrical factor is the same as the factor which factorizes the expression for the average current into the $i$-th patch in Eq. (28). This result is remarkable for two reasons.

First, it means that the geometry affects different gases in the same way for the purpose of calculating the average current (28) and the cross-correlation (57). For instance, changing the geometry for one gas model so as to increase the average current and the cross-correlation will have the same effect on any other gas model, in spite of their different microscopic dynamics. This property should be contrasted with the cross-correlations of currents into different patches, Eqs. (46) and (41). There the geometry is not factorized out, so that a change in the geometry affects different gas models differently.

Second, the geometry dependence of the average current (28) and the cross-correlation (57) is factorized out via the same geometrical factor $\sum_{j=1}^{s}A_{ij}$. That is, a given gas in two different geometries but sharing the same average current into the $i$-th patch will have the same cross-correlation between this specific current and the total current. This should again be contrasted with the cross-correlation which in general is not the same for a specified gas in two different geometries sharing the same average currents into the patches.

IV.4.2 Cell sensing by multiple receptors

Let us apply our results to the Berg-Purcell model berg modified by an account of interactions. In this model the cell is a sphere of radius $a$ covered with $s\gg 1$ disk-shaped absorbing receptors with a small radius $b\ll a$, so that only a small fraction of the cell surface is covered by the disks: $sb^{2}/4a^{2}<<1$. The cell is immersed in a gas of diffusing molecules at density $\rho_{0}$. Berg and Purcell approximated the total capacitance of the system (31) as $A\simeq 4\pi sab/(sb+\pi a)$. Then, assuming that the diffusing molecules are non-interacting, they evaluated the average total current (30) of molecules into the receptors:

[TABLE]

where $\bar{n}_{\text{full}}$ is the average steady state current into a fully absorbing sphere. Based on Eq. (58), Berg and Purcell concluded the following: “For large $s$ the intake approaches that of a completely absorbing cell, as it ought to. But it can become almost that large before more than a small fraction of the cell’s surface is occupied by absorbent patches”. Berg and Purcell gave an elegant explanation of Eq. (58) in terms of a trajectory of a single diffusing molecule. Although this single-particle picture breaks down for interacting molecules, Eq. (58) remains intact. Indeed, as Eq. (30) shows, the geometry is factorized out. It is also factorized out, for any gas model, in Eq. (56) for the variance of the total current. Therefore, if the number of absorbing patches is increased, but the area fraction of the absorbers is kept constant, fluctuations in the total current grow. Furthermore, our expression (57) for the cross-correlation shows how the number of patches affects cross-correlations:

[TABLE]

where we have used Eq. (57), and $\overline{\delta n^{2}}/\bar{n}^{2}|_{\text{full}}$ is the (normalized) variance of the total current of a fully absorbing sphere. Equation (59) shows that, as the number of patches increases, the cross-correlations, normalized by the mean current, decrease. This result is to be expected on physical grounds. Importantly, it holds for any (in general, interacting) diffusive lattice gas.

V Optimal density and flux fields

Here we determine the optimal profiles of the gas density and flux fields conditioned on specified absorption currents into each patch. Let us start with the gas of non-interacting RWs. Here the calculations are straightforward, see Appendix F. The resulting stationary optimal density profile, $\rho_{\text{RWs}}(\mathbf{x})=\bar{\rho}_{\text{RWs}}(\mathbf{x})+\rho_{1_{\text{RWs}}}(\mathbf{x})$ can be written as

[TABLE]

where, using Eqs. (23) and (27), we have:

[TABLE]

The optimal flux field is (see Appendix F):

[TABLE]

Let us calculate the vorticity $\nabla\times\boldsymbol{J_{\text{RWs}}}$. Taking the curl of both sides of Eq. (62) we obtain, after cancellation of two terms,

[TABLE]

This quantity is, in general, non-zero. That is, the optimal fluctuating flux field is, in general, not a potential vector field, in contrast to the average flux field (24).

Equation (63) can be generalized to an arbitrary lattice gas. For small fluctuations the flux field is equal to

[TABLE]

see Eq. (32). The first two terms on this right-hand-side are potential vector fields, so a non-zero contribution to $\nabla\times\boldsymbol{J}$ can come only from the last term. Substituting $p_{1}$ from Eq. (38) and using Eqs. (19) and $(\ref{phii})$, we obtain

[TABLE]

The explicit result (63) for the RWs is a particular limit of this relation. As is clear from Eq. (64), for a fluctuating system (that is, when $\lambda_{i}\neq 0$ for some $i$) the vorticity vanishes if and only if the two vector fields, $\sum_{i=1}^{s}\lambda_{i}\nabla\phi_{i}$ and $\sum_{i=1}^{s}\nabla\phi_{i}$, entering Eq. (64), are parallel. This happens in a single-current system, that is if there is only one absorbing patch. (This also happens in the extensively studied finite two-reservoir system, also sustaining a single current.) The vorticity also vanishes if all $\lambda_{i}$ are equal to each other. As shown in the previous Sec. IV.4, this special case appears when the process is conditioned on the total absorption current into all patches combined, $n=\sum_{i=1}^{s}{n_{i}}$. We will return to this special case in Sec. V.2.

Let us examine the large-$|\mathbf{x}|$ asymptotic of the optimal vorticity field. The fields $\nabla\phi_{i}$, comprising the cross product in Eq. (64), can be expanded in multipoles, see e.g. Ref. jec , Chap. 4. The two monopole terms, each decaying as $|\mathbf{x}|^{-2}$, do not contribute to the cross product, as they are directed in the radial direction. The leading contribution to the vorticity comes from the cross product of a monopole field of one potential and a dipole field of the other. This product decays as $|\mathbf{x}|^{-5}$ regardless of the problem’s geometry or specific lattice gas.

Bodineau et al. partial considered a finite two-dimensional lattice gas in contact with two particle reservoirs. They studied fluctuations of the partial current through an imaginary slit located in the bulk and showed that the fluctuations are dominated by point-like flux vortices localized at the edges of the slit. The action, evaluated over these solutions, can be made arbitrarily small indicating a breakdown of the MFT. In our system the flux-field vorticity is macroscopic, and is generated in the bulk, while point-like vortices are forbidden by the boundary conditions on the absorbing-reflecting surface.

V.1 Two hemispheres

To illustrate our results on the optimal density and flux fields, we consider a gas of RWs in contact with a sphere of radius $R$ placed in the origin. The whole sphere is absorbing, and we are interested in the absorption statistics of each of the two hemispheres: the northern one, $\theta\in\left[0,\pi/2\right]$, which we denote by $\Omega_{1}$ and the southern one, $\theta\in\left[\pi/2,\pi\right]$, denoted by $\Omega_{2}$. Here $\theta$ is the polar angle of the spherical coordinates. The problem possesses cylindrical symmetry. The effective potential, defined in Eq. (19), is $\phi=D_{0}\rho_{0}R/r$, where $r$ is the radial coordinate. The average absorption currents (25) to each of the two hemispheres are identical and equal to $\bar{n}_{1}=\bar{n}_{2}=2\pi D_{0}\rho_{0}R$. The set of effective potentials $\phi_{i}$ can be found explicitly via expansion in spherical harmonics:

[TABLE]

where the plus and minus signs refer to the northern and southern hemispheres, respectively. The function $f(r,\theta)$ is given by

[TABLE]

$P_{2k-1}(\dots)$ are the Legendre polynomials, and $\Gamma(\dots)$ is the gamma function. Plugging expressions (65) and (66) in Eq. (63) we obtain, after some algebra, the vorticity:

[TABLE]

where $\hat{\varphi}$ is the unit vector in the azimuthal direction $\varphi$ in spherical coordinates. Using Eq. (66), we can obtain the large-$r$ behavior of the vorticity:

[TABLE]

The $r^{-5}$ dependence is expected from the general argument given above. The vorticity exhibits a delta-function singularity at the equator of the sphere:

[TABLE]

We can also find the optimal flux and optimal density fields in this example by plugging expressions (65) and (66) in Eqs. (60)-(62). Figure 2 shows a vertical cross-section of the optimal density field $\rho_{\text{RWs}}$ for $n_{1}/\bar{n}_{1}=2$ and $n_{2}/\bar{n}_{2}=0.5$. As intuitively expected, the optimal density is higher/lower than the mean density close to the northern/southern hemisphere, respectively. The optimal flux field (62) and the vorticity field (67) in this case are shown in Figs. 3 and 4, respectively. Evident in Fig. 4 is the singularity of the vorticity at $r\rightarrow R$, displayed by Eq. (69). The flux field (62) exhibits a milder singularity at $r\rightarrow R$, where the radial component of the flux is discontinuous as a function of $\theta$ at $\theta=\pi/2$. These singularities are due to the discontinuous boundary conditions for the potentials $\phi_{i}$’s, and they are generic when considering several patches (or reservoirs) in direct contact with each other. The density field (60) is continuous at $r=R$.

V.2 Optimal profiles conditioned on the total absorption current

The optimal density and flux fields, conditioned on the total absorption current deviation $\delta n$, are given by the solution of the linearized MFT equations (32) and (33) alongside with Eqs. (51) and (52). In this case the equations can be explicitly solved for any lattice gas, see Appendix G, and we obtain

[TABLE]

where

[TABLE]

Note that $\phi\left(\mathbf{x}\right)=\tilde{I}_{1}\left[\bar{\rho}\left(\mathbf{x}\right),\rho_{0}\right]$, see Eq. (19). As Eqs. (70) and (71) show, both $\rho(\mathbf{x})=\bar{\rho}(\mathbf{x})+\rho_{1}(\mathbf{x})$, and $p(\mathbf{x})=p_{1}(\mathbf{x})$ depend on $\bar{\rho}(\mathbf{x})$ alone. Remarkably, these $\bar{\rho}(\mathbf{x})$-dependencies are described by (the small-fluctuation limit of) the optimal profiles $\rho^{(1)}(x)$ and $p^{(1)}(x)$ of the one-dimensional problem for the reservoir density values $\rho_{a}=\rho_{0}$ and $\rho_{b}=0$, found in Ref. bd :

[TABLE]

where $\phi\left(\bar{\rho}\right)$ is defined in Eq. (19). As one can see from Eq. (LABEL:q1d), the average density $\bar{\rho}(\mathbf{x})$ plays the role of the natural “spatial coordinate” of the problem. (The same property was uncovered and exploited in a different geometry in Ref. shpiel .) With respect to this spatial coordinate the problem is effectively one-dimensional, hence the relation to the current fluctuations of the one-dimensional system, described at the end of Sec. IV.4.

Finally, using Eqs. (70) and (71), we obtain the optimal fluctuating flux field (32), conditioned on the total absorption current:

[TABLE]

where $\bar{\mathbf{J}}$ is the average steady state flux field (24). As Eq. (74) shows, in order to generate a total absorption current which is larger by a factor $n/\bar{n}$ than the average one, the system has to increase the optimal flux by the same factor everywhere in space (see also Refs. fullabsorb ; shpiel ). It is not surprising, therefore, that the optimal flux field, conditioned on the total absorption current, is vortex-free.

VI Discussion

In this work we employed the Macroscopic Fluctuation Theory (MFT) to determine the statistics of particle absorption by several patches located on the surface of a domain immersed in a gas composed of interacting diffusing particles. Essentially, this work extends a number of previous results Derrida2007 ; MFTreview ; bd ; fullabsorb ; shpiel , obtained for different single-current settings, to a multiple-current setting. Among central results of this work are Eqs. (41) and (47) for the covariance matrix of the absorption currents into different patches, and Eq. (48) which establishes, independently of the system’s geometry, whether the absorption currents correlate or anti-correlate. The same condition (48) has recently gained much attention in a different context – as a sufficient condition for the validity of the additivity hypothesis for arbitrary currents main ; ber ; main2 . This coincidence hints at a possible relation between correlations and dynamical phase transitions via which the additivity property breaks down.

One particular example of diffusive transport, extensively studied in the past, describes electronic transport in mesoscopic wires meso . In this case Eq. (2), with $D$ and $\sigma$ corresponding to the SSEP, serves as a suitable mathematical description jor . As follows from Eq. (48), the cross-correlations in this system are strictly negative. This result has been known for some time: it was obtained both by the scattering matrix method but and by the effective Langevin description multi . It has been verified in experiment and gained much attention, see e.g. f . It is often referred to as the electronic Hanbury Brown and Twiss effect meso , and is intimately related to the Fermi statistics of charge carriers. The MFT framework which we developed here enables one to address a much broader class of diffusive systems by considering more general diffusivity and mobility. As we have seen, this may lead to qualitatively different fluctuational behaviors, such as different signs of the correlation terms.

An important example of positive correlations in the (energy) absorption appears in the context of diffusive wave propagation in disordered media penini ; main . Indeed, mesoscopic wave transport can often be described by Eq. (2) which relates the wave energy flux to the local wave intensity $I$: an equivalent of particle density $\rho$. Remarkably, this “lattice gas” can be described in terms of the KMP model. Indeed, the corresponding diffusion coefficient is constant here, but the mobility behaves as $\sigma(I)\propto I^{2}$ penini . In this case our Eq. (48) predicts strictly positive cross-correlations in the photon absorption. It would be very interesting to test this prediction in experiment. The experimental setting can be similar to that of Ref. shapiro2 , where a waveguide filled with a disordered material was used to study the light-intensity correlations between two distant points along the waveguide.

Another important advantage of the MFT is that it predicts, for all settings, the optimal (most probable) spatial profiles of the density and of the flux, conditioned on a given multiple-current statistics. As we have shown here, the corresponding optimal flux field generically exhibits vorticity. The vorticity disappears if the process is conditioned on the total absorption current into all the patches combined. It would be interesting to see whether the vorticity still disappears for the total absorption current if one goes beyond typical, small fluctuations which we addressed here. We found that, under the assumption that the optimal flux field is vortex-free, the simple solution (LABEL:q1d) (which solves the time-independent nonlinear MFT equations exactly), is unique. As a result all the cumulants [and not only the first two as in Eq. (56)], are equal to the corresponding one-dimensional cumulants times the capacitance of the system. The irrotational character of the flux field is, however, a strong assumption. Its breaking, at a critical value of the total absorption current, may have a character of phase transition.

The MFT formalism can be readily extended to other geometries sustaining multiple currents in systems of interacting diffusing particles. It can also be extended to finite two-dimensional systems where a set of characteristic potentials $\phi_{i}$ and the capacitance matrix can always be defined. An infinite two-dimensional system does not reach a true steady state, and logarithmic corrections to the linear scaling of the action with time $T$ are to be expected b . Last but not least, the formalism can also be extended, with some modifications, to the case when different reservoirs are kept at different densities.

ACKNOWLEDGMENTS

We are grateful to Eric Akkermans, Yaron Bromberg, Giovanni Jona-Lasinio and Naftali Smith for useful discussions of different parts of this work. We acknowledge financial support from the Israel Science Foundation (grant No. 807/16) and the United States-Israel Binational Science Foundation (BSF) (grant No. 2012145).

Appendix A Derivation of the MFT equations and boundary conditions

There are several methods for derivation of the MFT equations MFTreview . Here we use the Martin-Siggia-Rose formalism MFTreview ; MSR ; deridamft ; map . We start from Eq. (2) and represent the probability of observing a joint density and flux history $\rho(\mathbf{x},t),\mathbf{J}(\mathbf{x},t)$, constrained by the conservation law $\partial_{t}{\rho}+\nabla\cdot\mathbf{J}=0$, as a path integral:

[TABLE]

Using an integral representation for the $\delta$-function with the help of an auxiliary field $p(\mathbf{x},t)$, we can rewrite Eq. (A) as a path integral over three unconstrained fields map :

[TABLE]

We wish to evaluate the path integral over only those histories which led to $N_{i}$ particles being absorbed, by time $T$, by the $i$-th patch, $i=1,2,\dots,s$. Assuming that all characteristic length scales are macroscopic and include large numbers of particles, we can evaluate this path integral via a saddle-point approximation. The dominant contribution comes from the optimal fluctuation: the most probable history $(\rho,\mathbf{J},p)$ leading to a specified number of absorbed particles $N_{i}$. The problem, therefore, reduces to finding the minimum of $\mathcal{L}$ under the $s$ constraints:

[TABLE]

The latter can be incorporated via $s$ Lagrange multipliers $\lambda_{i}$, $i=1,2,\dots,s$:

[TABLE]

where the a priori unknown Lagrange multipliers are ultimately set by the $s$ constraints (77). Taking the first variation of $S$ we obtain:

[TABLE]

Here we have used the conditions

[TABLE]

following from the boundary conditions (4)-(7). We have also used Gauss’s theorem to transform volume integrals to surface integrals. There are no contributions from the surface integral at infinity.

Setting $\delta S$ to zero and using the fact that $\delta\rho(\mathbf{x}\in\Omega_{r},t)$, $\delta\rho(\mathbf{x},t=T)$, and $\boldsymbol{\delta J}(\mathbf{x}\in\Omega_{i},t)$ can be arbitrary, we see that the the optimal flux field is given by

[TABLE]

and the optimal profiles $\rho$ and $p$ satisfy the MFT equations (8) and (9) in the bulk and the boundary conditions (10)-(16). The MFT equations (8) and (9) are Hamiltonian, where $\rho$ and $p$ play the role of the conjugate “coordinate” and “momentum” density fields. The Hamiltonian is

[TABLE]

where the Hamiltonian density is given by MFTreview

[TABLE]

Finally, the action (78) can be simplified to (17) by using Eqs. (8), (77) and (81).

Appendix B Derivation of Eq. (43)

The formal solution to Eq. (42) is given by Green’s function $G(\mathbf{x},\mathbf{x^{\prime}})$ which satisfies the equations

[TABLE]

$G(\mathbf{x},\mathbf{x^{\prime}})$ can be interpreted as the electrostatic potential at point $\mathbf{x}$, induced by a point charge $q=-1/4\pi$ at point $\mathbf{x^{\prime}}$ outside the domain, when all the patches $\Omega_{i}$ are conducting and grounded. With $G(\mathbf{x},\mathbf{x^{\prime}})$, the solution to Eq. (42) can be written as

[TABLE]

Now let us evaluate the flux deviation

[TABLE]

whose surface integrals over the patches give the corresponding current deviations:

[TABLE]

Using a Green’s function identity

[TABLE]

jec , we obtain

[TABLE]

As the last step, we use Gauss’s theorem to transform the surface integral in Eq. (B) to a volume integral with the help of the potential $\phi_{i}$ defined in Sec. III:

[TABLE]

where in the last equality we have substituted $p_{1}$ from Eq. (38). This linear relation can be written in a matrix form as in Eq. (43).

Appendix C Matrix $\boldsymbol{C}$ is positive definite

To prove this statement one needs to show that, for every nontrivial $s$-dimensional vector $\boldsymbol{V}$, $\boldsymbol{V}^{T}\cdot\boldsymbol{C}\cdot\boldsymbol{V}>0$. We have

[TABLE]

which is clearly non-negative. Furthermore, since the fields $\nabla\phi_{i}$ are linearly independent, it is positive definite.

Appendix D Derivation of Eq. (47)

We start with rewriting (41):

[TABLE]

where we have used the fact that $\phi_{i}$-s are harmonic functions. Using Gauss’s theorem, we can transform the last integral into a surface integral. There is no contribution from the surface integral at infinity since both $\nabla\left(\phi_{i}\phi_{j}\right)$, and $\phi_{i}\phi_{j}\nabla\sigma(\bar{\rho})$ decay as $\mathbf{x}^{-3}$ as $\mathbf{x}\rightarrow\infty$. Therefore, the last integral in Eq. (LABEL:2) becomes

[TABLE]

As

[TABLE]

the surface integral is taken over the absorbing patches $\Omega_{i}$ alone. Moreover, since

[TABLE]

the surface integral of the first term vanishes. Finally, since $\phi_{i}\left(\mathbf{x}\in\Omega_{j}\right)=\delta_{i,j}$ we obtain

[TABLE]

where in the last equality we have used the fact that, on the absorbing patches, $\bar{\rho}=const=0$, and identified the average currents (25) expressed through the average flux field (24). What is left to arrive at Eq. (47) is to rewrite the integrand of the second integral in Eq. (LABEL:2) as

[TABLE]

where we have used Eq. (19) and the fact that $\phi$ is a harmonic function.

Appendix E Conditioning on the total absorption current

We seek the minimum of the bilinear form

[TABLE]

subject to the constraint

[TABLE]

Introducing an additional Lagrange multiplier $\gamma$, we minimize the function

[TABLE]

with respect to all $\delta n_{i}$. The minimization yields

[TABLE]

Comparing this expression with Eq. (43), we see that

[TABLE]

Now we can find the currents and all $\lambda_{i}$-s in terms of the total absorption current $n$. We plug the identity $\sum_{j}\phi_{j}=\phi/V(\rho_{0})$ in Eq. (41) and obtain

[TABLE]

Now we introduce a single-variable function $\bar{\rho}(\phi)$ as the inverse function to $\phi(\bar{\rho})$ defined in Eq. (19). Using $\bar{\rho}(\phi)$, we can rewrite the integrand as

[TABLE]

where we have also used the fact that $\phi_{i}$ is a harmonic function. Applying Gauss’s theorem for this form of the integrand in Eq. (97) we have (the surface integral at infinity vanishes due to our choice of the lower bound of the integral on $\sigma$, and that $\nabla\phi_{i}$ decays as $\mathbf{x}^{-2}$):

[TABLE]

where in the second and third equalities we used the boundary conditions:

[TABLE]

respectively. In the forth equality we employed the identity

[TABLE]

which reflects the symmetry of the capacitance matrix (26). As a result,

[TABLE]

where $\alpha(\rho_{0})$ is defined in (53). It is left to impose the constraint $\sum_{i}\delta n_{i}=\delta n$ in order to set the value of $\gamma$. This yields

[TABLE]

Equations (99) and (100) yield Eqs. (51) and (52).

Appendix F Optimal density and flux fields for the RWs

For the non-interacting RWs one can solve Eq. (42) explicitly. Since both $\phi_{i}$ and $\bar{\rho}_{\text{RWs}}$ are harmonic functions [see Eq. (23)], we can rewrite the r.h.s. of Eq. (42) as

[TABLE]

Therefore, Eq. (42) can be rewritten as a Laplace’s equation:

[TABLE]

By virtue of the boundary conditions for $\bar{\rho}$, $\phi_{i}$ and $\rho_{1}$, the function under the Laplacian vanishes on all the patches and at infinity, whereas its normal derivative vanishes on the reflecting part of the boundary $\Omega_{r}$. As a result, this function vanishes everywhere, and we obtain

[TABLE]

Now we use Eq. (43) alongside with $C_{ij}=\delta_{i,j}\bar{n}_{i}$ (see Sec. IV.2) and obtain:

[TABLE]

Plugging this relation into Eq. (102) we arrive at Eq. (60). Also, using Eq. (102), we can evaluate the flux deviation

[TABLE]

Substituting here $p_{1}$ from Eq. (38) and using (103), we arrive at Eq. (62).

Appendix G Optimal profiles for the total absorption current

We start from $p_{1}$ from Eq. (38). All $\lambda_{i}$ are given by Eq. (52). Then, using Eqs. (27) and (53), we obtain:

[TABLE]

where $\tilde{I}_{1}(\bar{\rho},\rho_{0})$ is defined by Eq. (72). Now we turn to Eq. (42) for $\rho_{1}$. Using

[TABLE]

we obtain

[TABLE]

As a result, Eq. (42) for $\rho_{1}$ becomes Laplace’s equation,

[TABLE]

with inhomogeneous boundary conditions on the absorbing patches [as can be deduced from the boundary conditions (34)-(37) for $\rho_{1}$ and the definition (72) of $\tilde{I}_{2}$]. The solution gives $\rho_{1}$ as a function of the average density field $\bar{\rho}(\mathbf{x})$ in Eq. (70).