*Article* **Joint Fluctuation Theorems for Sequential Heat Exchange**

#### **Jader Santos 1, André Timpanaro <sup>2</sup> and Gabriel Landi 1,\***


Received: 5 March 2020; Accepted: 8 July 2020; Published: 12 July 2020

**Abstract:** We study the statistics of heat exchange of a quantum system that collides sequentially with an arbitrary number of ancillas. This can describe, for instance, an accelerated particle going through a bubble chamber. Unlike other approaches in the literature, our focus is on the *joint* probability distribution that heat *Q*<sup>1</sup> is exchanged with ancilla 1, heat *Q*<sup>2</sup> is exchanged with ancilla 2, and so on. This allows us to address questions concerning the correlations between the collisional events. For instance, if in a given realization a large amount of heat is exchanged with the first ancilla, then there is a natural tendency for the second exchange to be smaller. The joint distribution is found to satisfy a Fluctuation theorem of the Jarzynski–Wójcik type. Rather surprisingly, this fluctuation theorem links the statistics of multiple collisions with that of independent single collisions, even though the heat exchanges are statistically correlated.

**Keywords:** fluctuation theorems; collisional models

#### **1. Introduction**

Fluctuations of thermodynamic quantities, which are usually negligible in macroscopic systems, are known to play a dominant role in the micro- and mesoscopic domain. These fluctuations are embodied in the so-called fluctuation theorems (FT) [1–4], a collection of predictions for systems evolving under nonequilibrium conditions valid beyond linear response. They can be summarized as [5,6]

$$\frac{P(+\Sigma)}{\overline{P}(-\Sigma)} = \epsilon^{\Sigma},\tag{1}$$

where *P*(Σ) denotes the probability that an amount of entropy Σ is produced in a certain process and *P*˜(Σ) denotes the corresponding probability for the time-reversed process.

Of the many scenarios which present FTs, one which is particularly interesting is that of heat exchange between a system *S*, prepared in equilibrium with a temperature *Ts*, and an environment *E*, prepared in a different temperature *Te*. In this case, as first shown by Jarzynski and Wójcik in Ref. [7], the distribution *P*(*Q*) of the heat exchanged between them, satisfies

$$\frac{P(+Q)}{\overline{P}(-Q)} = e^{\Lambda \beta Q},\tag{2}$$

where Δ*β* = *β<sup>e</sup>* − *β<sup>s</sup>* (with *β* = 1/*T* and *kB* = 1). Here, and throughout the paper, *Q* denotes the net heat transfer from the system to the environment. Quite surprising, in this case it turns out that *P*˜(*Q*) = *P*(*Q*), meaning the statistics of the forward and backward processes are the same. Equation (2) was subsequently generalized to allow for the exchange of both energy and particles between several interacting systems initially at different temperatures and chemical potentials [6,8,9].

Here we consider a generalization of this scenario, where the system interacts sequentially with multiple parts of the environment, exchanging heat with each part. One can imagine, for instance, an accelerated particle crossing a bubble chamber. In this case, the system will leave a trail on *E*, represented by the heat exchanged in each point. In the microscopic domain this process will be stochastic, with a random amount of heat exchanged in each interaction.

The key idea that we will explore in this paper is to look at the joint probability distribution for the heat exchanged with each part, *P*(*Q*1, *Q*2, *Q*3, ...). This allows us to understand the correlations between the different heat exchanges.

For instance, in a situation where all the ancillas have the same temperature, from a stochastic perspective a large exchange in the first collision increases the probability that the second collision exchanges less. This feature is fully captured by the joint distribution. This happens because thermal operations have the property of bringing the system closer to its thermal equilibrium state, *σ*eq, i.e., [10]

$$D\left(\sigma\_{0}\|\sigma\_{\text{eq}}\right) \ge D\left(\sigma\_{1}\|\sigma\_{\text{eq}}\right) \ge D\left(\sigma\_{2}\|\sigma\_{\text{eq}}\right) \ge \dots \ge D\left(\sigma\_{N}\|\sigma\_{\text{eq}}\right) \tag{3}$$

where *D*(*ρ* ||*ρ*) = Tr(*ρ* ln *ρ* − *ρ* ln *ρ*) is the quantum relative entropy. If in the first interaction the system exchange a large quantity of heat, the system gets a lot closer to its steady state. So in the next interaction, the system should exchange less heat.

To formalize this idea, we split the environment into a set of ancillas *Ai*, with which the system interacts sequentially, producing a collisional model [11–14]. The process is schematically illustrated in Figure 1 and the formal framework is developed in Section 2. In Section 3 we then show that *P*(*Q*1, *Q*2, *Q*3, ...) satisfies a fluctuation theorem that generalizes (2). Moreover, we show how this fluctuation theorem relates the joint distribution to the statistics of a single collision, even though the events are statistically correlated.

**Figure 1.** Schematic representation of a system *S* interacting sequentially with a series of ancillas. The system starts in the state *σ*<sup>0</sup> and the ancillas in an initial states *ρi*, which are assumed to be thermal but at possibly different temperatures. Each *SAi* interaction is also governed by a possibly different unitary *Ui*.

#### **2. Formal Framework**

We consider a quantum system *S*, with Hamiltonian *H<sup>s</sup>* , prepared in a thermal state *σ*<sup>0</sup> = *e*−*βsH<sup>s</sup>* /*Zs*, with temperature *Ts*. The system is put to interact sequentially with a series of *N* ancillas *Ai*, as depicted in Figure 1. The ancillas are not necessarily identical. Each has Hamiltonian *H<sup>i</sup>* and is prepared in a thermal state *<sup>ρ</sup><sup>i</sup>* = *<sup>e</sup>*−*βiH<sup>i</sup>* /*Zi*, with possibly different temperatures *Ti*. Each collision is described by a unitary operator *Ui* acting only on *SAi*, which may also differ from one interaction to another.

In order to comply with the scenario of Ref. [7], we assume that the *Ui* satisfy the strong energy-preservation condition

$$\left[\mathcal{U}\_{\dot{\iota}\prime}H^{s} + H^{\dot{\iota}}\right] = 0.\tag{4}$$

Or, what is equivalent, that each collision is a thermal operation [10,15]. This implies that all energy that leaves *S* enters *Ai*, so nothing is stuck in the interaction. As a consequence, there is no work involved and all the change in energy of the system can be unambiguously identified as heat flowing to the ancillas [13].

We label the eigenvalues and eigenvectors of the system as *<sup>H</sup>s*|*α* = *<sup>E</sup><sup>s</sup> <sup>α</sup>*|*α*. For concreteness, we assume these levels are non-degenerate. Time is labeled discretely by *i* = 1, 2, 3, ..., representing which collisions already took place. For instance, the initial state is decomposed as *σ*<sup>0</sup> = ∑*α*<sup>0</sup> *p*0(*α*0)|*α*0*α*0|, with *p*0(*α*0) = *e* <sup>−</sup>*βsE<sup>s</sup> <sup>α</sup>*<sup>0</sup> /*Zs* and we use *α*<sup>0</sup> to emphasize that this is before the first collision. Similarly, the eigenvalues and eigenvectors of the ancillas are labeled as *Hi*|*ni* = *<sup>E</sup><sup>e</sup> ni* |*ni*. The initial state of each *Ai* is thus decomposed as *ρ<sup>i</sup>* = ∑*ni qi*(*ni*)|*nini*| where *qi*(*ni*) = *e* <sup>−</sup>*βiE<sup>e</sup> ni* /*Zi*.

The dynamics depicted in Figure 1 generates a stroboscopic map for the system. The joint state of *SAi* after the interaction is given by

$$
\varrho\_i = \mathcal{U}\_i \big( \sigma\_{i-1} \otimes \rho\_i \big) \mathcal{U}\_i^\dagger. \tag{5}
$$

Taking the partial trace over *Ai* then leads to the updated state *σi*. Conversely, tracing over the system leads to the reduced state *ρ <sup>i</sup>* of the ancilla after the interaction (Figure 1).

The fact that the unitary is energy preserving (Equation (4)), together with the assumption that the energy levels are non-degenerate, means that it is possible to construct quantum trajectories for the system in two equivalent ways. The first is to assume a two-point measurement scheme in *S* at each step [16,17]. Equation (4) implies that the system will remain diagonal in the energy basis, so that measurements in this basis are non-invasive (that is, have no additional entropy production associated to it). Measuring *S* in the energy basis after each collision then leads to the trajectory

$$\gamma\_s = \{\mathfrak{a}\_0, \mathfrak{a}\_1, \dots, \mathfrak{a}\_N\}.\tag{6}$$

The heat associated with each collision is then readily defined as

$$Q\_i[\gamma\_\delta] = -E\_{a\_i}^s + E\_{a\_{i-1}\prime}^s \tag{7}$$

Alternatively, one can construct a quantum trajectory by measuring the ancillas, before and after each collision, plus a single measurement of the system before the process starts. That is, one can consider instead a quantum trajectory of the form

$$\gamma\_{\varepsilon} = \{n\_0, n\_1, n\_1', n\_2, n\_2', \dots, n\_N, n\_N'\}. \tag{8}$$

This, in a sense, is much more natural since the ancillas are only used once and thus may be experimentally more easily accessible. Furthermore, as far as heat exchange is concerned, this turns out to be equivalent to the trajectory (6). The reason is that Equation (4) implies the restriction

$$\delta\left \propto \delta\left( (E\_{n\_i}^s + E\_{n\_i'}^c) - (E\_{n\_{i-1}}^s + E\_{n\_i}^c) \right) \tag{9}$$

where *δ*(*x*) is the Kronecker delta. In addition, since the energy values are taken to be non-degenerate, energies uniquely label states. Thus, for instance, if we know *α*0, *n*1, *n* <sup>1</sup> we can uniquely determine *α*1, and so on. The converse, however, is not true: from *α*<sup>0</sup> and *α*<sup>1</sup> we cannot specify *n*<sup>1</sup> and *n* <sup>1</sup> (which is somewhat evident given that the number of points in Equation (6) is smaller than that in Equation (8)). This, however, is not a problem if one is interested only in the heat exchanged, which can also be defined from the trajectory (8) as

$$\mathbb{Q}Q\_i[\gamma\_t] = E\_{n\_i'}^{\varepsilon} - E\_{n\_i}^{\varepsilon}. \tag{10}$$

Due to Equation (9), this must coincide with Equation (7); i.e., *Qi*[*γe*] ≡ *Qi*[*γs*].

The assumption in Equation (4) may at first seem somewhat artificial. However, this is not the case. This assumption is a way to bypass the idea of weak coupling, which is one of the conditions used in [7]. Moreover, the interesting thing about the present analysis is that it establishes under which conditions Equations (6) and (8) are equivalent. Naively one would expect that this is often the case. However, as the above arguments show, several assumptions are necessary for this to be the case. This reflects some of the challenges that appear in describing thermodynamics in the quantum regime.

#### *2.1. Path Probabilities from Measurements in S*

Thermal operations imply that the probability that, after the *i*-th collision, the system is in a given eigenstate |*αi* depends only on the probabilities in the previous time. That is, the dynamics of populations and coherences completely decouple [18]. Indeed, Equation (5) together with Equation (4) imply that

$$p\_i(a\_i) = \langle a\_i | \sigma\_i | a\_i \rangle = \sum\_{a\_{i-1}} M\_i(a\_i | a\_{i-1}) p\_{i-1}(a\_{i-1}) \, . \tag{11}$$

where

$$M\_i(a\_i|\alpha\_{i-1}) = \sum\_{n\_i,n\_i'} |\langle a\_i,n\_i'|\boldsymbol{l}I\_i|a\_{i-1},n\_i\rangle|^2 q\_i(n\_i). \tag{12}$$

The populations therefore evolve as a classical Markov chain, with *Mi*(*αi*|*αi*−1) representing the transition probability of going from *αi*−<sup>1</sup> to *αi*. Moreover, Equation (9) together with the fact that the ancillas are initially thermal, imply that *Mi*(*αi*|*αi*−1) satisfies detailed balance

$$M\_i(a\_i|a\_{i-1})e^{-\beta\_i E\_{a\_{i-1}}^\*} = M\_i(a\_{i-1}|a\_i)e^{-\beta\_i E\_{a\_{i\_s}}^\*} \tag{13}$$

where, notice, what appears here is the temperature *β<sup>i</sup>* of ancilla *Ai*.

The path probability associated with *γ<sup>s</sup>* in Equation (6) will then be

$$\mathcal{P}[\gamma\_s] = M\_N(\mathfrak{a}\_N|\mathfrak{a}\_{N-1})\dots M\_2(\mathfrak{a}\_2|\mathfrak{a}\_1)M\_1(\mathfrak{a}\_1|\mathfrak{a}\_0)p\_0(\mathfrak{a}\_0),\tag{14}$$

which is nothing but the joint distribution of a Markov chain. We call attention to the clear causal structure of this expression: marginalizing over future events has no influence on past ones. For instance, summing over *α<sup>N</sup>* leads to a distribution of the exact same form. Conversely, marginalizing over past variables completely changes the distribution.

The joint distribution of heat can then be constructed from Equation (14) in the usual way:

$$P(Q\_1, \ldots, Q\_N) = \sum\_{\gamma\_s} \mathcal{P}[\gamma\_s] \left( \prod\_{i=1}^N \delta \left( Q\_i - Q\_i[\gamma\_s] \right) \right). \tag{15}$$

This is the basic object that we will explore in this paper.

#### *2.2. Path Probabilities from Measurements in the Ai*

Alternatively, we also wish to show how Equation (15) can be constructed from the trajectory *γ<sup>e</sup>* in Equation (8). The easiest way to accomplish this is to first consider the augmented trajectory

$$\gamma\_{st} = \{a\_0, n\_1, n\_1', a\_1, n\_2, n\_2', a\_2, \dots, n\_{N'}, n\_{N'}', a\_N\} \tag{16}$$

Introducing the transition probabilities *Ri*(*αi*, *n i* |*αi*−1, *ni*) = |*αi*, *n i* <sup>|</sup>*Ui*|*αi*−1, *ni*|2, the path distribution associated with the augmented trajectory *γse* will be

$$\mathcal{P}[\gamma\_{sc}] = R\_N(a\_N, n\_N'|a\_{N-1}, n\_N) \dots R\_1(a\_1, n\_1'|a\_0, n\_1) q\_N(n\_N) \dots q\_1(n\_1) p\_0(a\_0) \dots$$

As a sanity check, if we marginalize this over *ni* and *n <sup>i</sup>* we find

$$\mathcal{P}[\gamma\_{\mathcal{S}}] = \sum\_{\substack{n\_1,\ldots,n\_N\\n'\_1,\ldots,n'\_N}} R\_N(a\_{N'}n'\_N|a\_{N-1},n\_N)\ldots R\_1(a\_1,n'\_1|a\_0,n\_1)q\_N(n\_N)\ldots q\_1(n\_1)p\_0(a\_0)$$

$$= M\_N(a\_N|a\_{N-1})\ldots M\_2(a\_2|a\_1)M\_1(a\_1|a\_0)p\_0(a\_0),$$

where we used Equation (12). This is therefore precisely P[*γs*] in Equation (14), as expected.

Instead, from P[*γse*] one can now obtain P[*γe*] by marginalizing over *α*1,..., *αN*; viz.,

$$\mathcal{P}[\gamma\_{\varepsilon}] = \sum\_{a\_1, \dots, a\_N} R\_N(a\_{N\prime} n\_N^{\prime} | a\_{N-1\prime} n\_N) \dots R\_1(a\_1, n\_1^{\prime} | a\_0, n\_1) q\_N(n\_N) \dots q\_1(n\_1) p\_0(a\_0) \dots \tag{17}$$

The above analysis puts in evidence the Hidden Markov nature of the dynamics in Figure 1. When measurements are done in the ancilla, the system plays the role of the hidden layer, which is not directly accessible. Instead, predictions about the system must be made from the visible layer (i.e., the ancillas).

This Hidden Markov nature manifests itself on the fact that even though the system obeys a Markov chain [Equation (14)], the same is not true for the ancillas. In symbols, this is manifested by the fact that *n <sup>i</sup>* depends not only on *ni* and *n <sup>i</sup>*−1, but on the entire past history (*n*1, *<sup>n</sup>* <sup>1</sup>, ... , *ni*−1, *n <sup>i</sup>*−1, *ni*). This is intuitive in a certain sense: the amount of heat exchanged at the *i*-th collision will depend on the heat exchanged in all past events.

With P[*γe*], the distribution of heat, Equation (15) can be equivalently defined using Equation (10). One then finds

$$P(Q\_1, \ldots, Q\_N) = \sum\_{\gamma\_\epsilon} \mathcal{P}[\gamma\_\epsilon] \left( \prod\_{i=1}^N \delta \left( Q\_i - Q\_i[\gamma\_\epsilon] \right) \right). \tag{18}$$

The reason why this is equivalent to Equation (15) becomes clear from the way we derived P[*γe*] above: we can expand the summation to *γse* and then use the fact that *Qi*[*γs*] = *Qi*[*γe*].

#### *2.3. Backward Process*

To construct the fluctuation theorem, we must now establish the backward process. As shown in [19], however, there is an arbitrariness in the choice of the initial state of the backward process; different choices lead to different definitions of the entropy production. Here we are interested specifically in heat and the generalization of the Jarzynski–Wójcik fluctuation theorem [7]. Hence, we assume that in the backward process both system and ancillas are fully reset back to their thermal states. As usual, the time-reversed interaction between *SAi* now takes place by means of the unitary *U*† *i* . However, the order of the interactions must now be flipped around, as shown in Figure 2. More about the choice of backward process can be found in [20,21] and its relation to the notion of recovery maps is discussed in [22].

**Figure 2.** Schematic representation of the backward process.

In the backward process, the system will therefore evolve according to

$$\vec{p}\_i(\mathfrak{a}\_{N-i}) = \sum\_{\mathfrak{a}\_{N-i+1}} M\_{N-i+1}(\mathfrak{a}\_{N-i}|\mathfrak{a}\_{N-i+1}) \vec{p}\_{i-1}(\mathfrak{a}\_{N-i+1}),$$

where we index the states as *αN*−*<sup>i</sup>* instead of *α<sup>i</sup>* just so that the trajectory *γ<sup>s</sup>* can remain the same as in the forward process. The path probability <sup>P</sup>˜ [*γs*] associated to this process will then be

$$\mathcal{P}[\gamma\_s] = M\_1(\mathfrak{a}\_0|\mathfrak{a}\_1)\dots M\_N(\mathfrak{a}\_{N-1}|\mathfrak{a}\_N)p\_0(\mathfrak{a}\_N)\,. \tag{19}$$

which is similar to that used in the original Crooks fluctuation theorem [23]. The corresponding heat distribution is

$$\mathcal{P}(Q\_N, \dots, Q\_1) = \sum\_{\gamma\_s} \mathcal{P}[\gamma\_s] \prod\_{i=1}^N \delta \left( Q\_i + Q\_i[\gamma\_s] \right), \tag{20}$$

where *Qi* continues to be the heat exchanged with *Ai* (which is now different from the heat exchanged at collision *i*).

#### **3. Joint Fluctuation Theorem for Heat Exchange**

We are now ready to construct the fluctuation theorem. The detailed balance condition (13) immediately implies that Equations (15) and (20) will be related by

$$\frac{P(Q\_1, \ldots, Q\_N)}{\tilde{P}(-Q\_N, \ldots, -Q\_1)} = e^{\sum\_{i=1}^N (\beta\_i - \beta\_s) Q\_i}.\tag{21}$$

This is a theorem for the joint distribution of the heat exchanged between multiple ancillas. It thus represents a generalization of Ref. [7] to the case where the system interacts sequentially with multiple reservoirs. This result has several features which are noteworthy. First, note that the temperature *β<sup>i</sup>* of the ancillas are not necessarily the same. Second, note how after the first collision the state of the system is no longer thermal. However, still, this does not affect the fluctuation theorem. All that matters is that before the first collision the system is in equilibrium.

It is also important to point out that any Markov chain satisfying the detailed balance relation also satisfies a fluctuation theorem [24]. This fact can be used to obtain Equation (21) when properly choosing the rates of the Markovian evolution. Beyond that, a generalization of the detailed FT to multiple reservoirs has also being obtained before, e.g., in Ref. [25].

#### *3.1. Causal Order and Relation to Single Collisions*

The causal order of the process plays a crucial role here. Marginalizing over future events has no effect on the fluctuation theorem. That is, from (21) one could very well construct a similar relation for *P*(*Q*1, ... , *QN*−1), by simply summing over *QN*. This is not possible, however, for marginalization over past events. That is, *P*(*Q*2,..., *QN*), for instance, does not satisfy a fluctuation theorem.

The right-hand side of Equation (21) is very similar to what appears in the original FT (2). We can make this more rigorous as follows. Let us consider a different process, consisting of a single collision between the system thermalized in *β<sup>s</sup>* and an ancilla thermalized in *β<sup>i</sup>* (Figure 3). The associated heat distribution *P*sc(*Qi*) will then satisfy Equation (2); viz.,

$$\frac{P\_{\rm sc}(Q\_i)}{P\_{\rm sc}(-Q\_i)} = e^{(\beta\_i - \beta\_s)Q\_i},\tag{22}$$

where, recall, in this case of a single collision the backward process coincides with the forward one, so that the distribution *P*˜ sc in the denominator is simply *P*sc. It is very important to emphasize, however, that *P*sc(*Qi*) is not the marginal of *P*(*Q*1, ... , *QN*) (with the exception of *Q*1). Notwithstanding, comparing with Equation (21), we see that the full process in Figure 1 is related to the single-collision processes according to

$$\frac{P(Q\_1, \ldots, Q\_N)}{P(-Q\_N, \ldots, -Q\_1)} = \frac{P\_{\rm sc}(Q\_1)}{P\_{\rm sc}(-Q\_1)} \cdots \frac{P\_{\rm sc}(Q\_N)}{P\_{\rm sc}(-Q\_N)}.\tag{23}$$

*Entropy* **2020**, *22*, 763

This result is noteworthy, for the right-hand side is a product whereas the left-hand side is not. The full distribution *P*(*Q*1, ... , *QN*) cannot be expressed as a product because the heat exchanges are, in general, not statistically independent. Notwithstanding, the ratio on the left-hand side of (23) does factor into a product. The point, though, is that this is not the product of the marginals, but of another distribution *P*sc.

One can also write a formula of the form (23), but for only some of the heat exchanges. For instance, it is true that

$$\frac{P(Q\_1, \ldots, Q\_N)}{\bar{P}(-Q\_N, \ldots, -Q\_1)} = \frac{P(Q\_1, \ldots, Q\_{N-1})}{\bar{P}(-Q\_{N-1}, \ldots, -Q\_1)} \frac{P\_{\text{sc}}(Q\_N)}{P\_{\text{sc}}(-Q\_N)}.\tag{24}$$

This kind of decomposition, however, depends crucially on the causal structure since it can only be done for future exchanges. For instance, we cannot write something involving *P*(*Q*2, ... , *QN*). The reason is that *P*(*Q*1, ... , *QN*−1) satisfies the fluctuation theorem (21), but *P*(*Q*2, ... , *QN*) does not (since, after the first collision the system is no longer in a thermal state).

**Figure 3.** Schematic representation of a single collision event.

#### *3.2. Information-Theoretic Formulation of the Entropy Production*

We define the entropy production associated with Equation (21) as

$$\Sigma[\gamma\_s] = \ln \frac{\mathcal{P}[\gamma\_s]}{\bar{\mathcal{P}}[\gamma\_s]} = \sum\_{i=1}^{N} (\beta\_i - \beta\_s) Q\_i[\gamma\_s]. \tag{25}$$

The second equality is obtained using the detailed balance relation (13). We emphasize that this is the entropy production associated with the choice of backward protocol used in Section 2.3, which may differ from other definitions in the literature [18,26]. As discussed in [19], the interpretation of the entropy production depends on the choice of the initial state of the backwards process. For instance, if we have chosen the initial state as the final state of the forward process, i.e., the state *<sup>N</sup>* (see Equation (5)), we would have a a contribution related to the correlations between the system and the ancillas. This type of entropy production was called the inclusive entropy production in Ref. [19]. This happens because this state carries the information about the correlations. Here we have choose a initial state for the backward process that does not have this contributions.

In [7], Jarzynski and Wójcik calculated an upper bound on the probability of observing a violation of the second law, i.e., the passage of heat from a colder to a hotter body. We can apply the same reasoning to Equation (25). Let us assume that all ancillas start in the same thermal state with temperature *Ta* and *β<sup>a</sup>* − *β<sup>s</sup>* > 0. The probability that the heat transfer from the system to *i*-th ancilla will fall below a specified value *qi* in each interaction through the whole process, obeys the inequality

$$\int\_{-\infty}^{q\_1} \mathrm{d}Q\_1 \cdots \int\_{-\infty}^{q\_N} \mathrm{d}Q\_N \ P(Q\_1, \dots, Q\_N) \le e^{(\beta\_4 - \beta\_s)(q\_1 + \dots + q\_N)}\tag{26}$$

which is the multiple-exchange extension of the result obtained in [7]. Equation (26) shows that observing a positive total transference of heat from the hot system to the cold ancillas dies exponentially with *q*<sup>1</sup> + ··· + *qN*.

Alternatively, we can consider the entropy production from the perspective of the global trajectory *γse* in Equation (16). Using also that *Qi*[*γs*] = *Qi*[*γe*], we can then write Σ[*γse*] as

$$\Sigma[\gamma\_{s\varepsilon}] = \sum\_{i=1}^{N} \beta\_i Q\_i[\gamma\_\varepsilon] - \beta\_s (E\_{a\_N}^s - E\_{a\_0}^s) = \sum\_{i=1}^{N} \ln \frac{q\_i(n\_i)}{q\_i(n\_i')} + \ln \frac{p\_0(a\_0)}{p\_0(a\_N)}.\tag{27}$$

The average entropy production may then be written as

$$
\langle \Sigma[\gamma\_{s\varepsilon}] \rangle = \mathcal{S}(\sigma\_N) - \mathcal{S}(\sigma\_0) + D(\sigma\_N||\sigma\_0) + \sum\_{i=1}^N \left\{ \mathcal{S}(\rho\_i') - \mathcal{S}(\rho\_i) + D(\rho\_i'||\rho\_i) \right\}, \tag{28}
$$

where *S*(*ρ*) = −Tr(*ρ* ln *ρ*) is the von Neumann entropy. Here *σ<sup>N</sup>* is the final state of the system after the *N* collisions. In the Equation (28), we can identify

$$S(\sigma\_N) - S(\sigma\_0) + \sum\_{i=1}^{N} S(\rho\_i') - S(\rho\_i) = \Delta I\_{\mathfrak{sl}} \tag{29}$$

where Δ*Ise* is the change in the mutual information between the system and the ancillas. This way we can have a more clear meaning of the expression (28). One term is proportional to the total correlations built between system and ancillas and the other two relative entropy terms measure the disturbance on the environment and the system during the process.

The important aspect of this result is that it depends only on local changes in the ancillas. That is, all quantities refer to the local states *ρ <sup>i</sup>* of each ancilla after the interaction. In reality, because the ancillas all interact with the system, they actually become indirectly correlated. These correlations are still represented indirectly in Σ[*γse*], but they do not appear explicitly. This, ultimately, is a consequence of the choice of backward process that is used in the Jarzynski–Wójcik scenario [7].

#### *3.3. Initially Correlated Ancillas*

One possible extension of our formalism is to consider the case of initially correlated system-ancillas. In this case, we could explore how the correlation between the system and the ancillas affect the XFT. This problem was studied for a single heat exchange in [27] and in our case, the same approach yields

$$\frac{\partial \left[\gamma\_{\mathcal{S}\varepsilon}\right]}{\partial \overline{\boldsymbol{\Im}} \{\gamma\_{\mathcal{S}\varepsilon}\}} = \varepsilon^{-\Delta T \left(\gamma\_{\mathcal{S}\varepsilon}\right) + \sum\_{i=1}^{N} \left(\beta\_{\varepsilon} - \beta\_{i}\right) \left(E\_{\boldsymbol{u}\_{i}}^{\boldsymbol{s}} - E\_{\boldsymbol{u}\_{i-1}}^{\boldsymbol{s}}\right) + \sum\_{i=1}^{N} \beta\_{i} \left[E\_{\boldsymbol{u}\_{i}}^{\boldsymbol{s}} + E\_{\boldsymbol{u}\_{i}^{\boldsymbol{s}}}^{\boldsymbol{s}} - \left(E\_{\boldsymbol{u}\_{i-1}}^{\boldsymbol{s}} + E\_{\boldsymbol{u}\_{i}}^{\boldsymbol{s}}\right)\right] \tag{30}$$

where the ΔI(*γse*) = *I*<sup>∗</sup> − *I* with

$$I^\* = \ln \left[ \frac{p(\alpha\_{n\_1} n'\_{1'}, \dots, n'\_N)}{p\_0(\alpha\_N) q\_1(n'\_1) \dots q\_N(n'\_N)} \right] \tag{31}$$

$$I = \ln\left[\frac{p(\boldsymbol{a}\_0, \boldsymbol{n}\_1, \dots, \boldsymbol{n}\_N)}{p\_0(\boldsymbol{a}\_0)q\_1(\boldsymbol{n}\_1)\dots q\_N(\boldsymbol{n}\_N)}\right] \tag{32}$$

where we define *p*(*α*0, *n*1, ... , *nN*) = *α*0, *n*1, ... , *nN*|*ρSE*|*α*0, *n*1, ... , *nN*. Here *ρSE* is the initial state for the system-ancillas. This result is similar to the one found in [27]. Because in our case, we are working with thermal operations, we can write Equation (30) as

$$\frac{P[\gamma\_{sc}]}{\bar{P}[\gamma\_{sc}]} = \varepsilon^{-\Delta \mathcal{T}(\gamma\_{sc}) + \sum\_{i=1}^{N} (\beta\_s - \beta\_i)(E\_{a\_i}^s - E\_{a\_{i-1}}^s)}\tag{33}$$

By taking the above equation and sum over all trajectories, to obtain the nonequilibrium equality for an initially correlated state

$$\langle e^{\Delta \mathcal{X} + \sum\_{i=1}^{N} (\beta\_i - \beta\_i) Q\_i} \rangle = 1 \tag{34}$$

and then using Jensen's inequality we have that

$$\sum\_{i=1}^{N} (\beta\_i - \beta\_s) \langle Q\_i \rangle \ge \langle \Delta \mathcal{Z} \rangle \tag{35}$$

So it is possible to obtain a type of Clausius relation where now the entropy production has a new lower bound.

#### **4. Conclusions**

To summarize, we have considered here the sequential heat exchange between a system and a series of ancillas. We assume all entities start in thermal state, but at possibly different temperatures. Moreover, all interactions are assumed to be described by thermal operations, which makes the identification of heat unambiguous. The main object of our study was the joint probability of heat exchange *P*(*Q*1, ... , *QN*) for a set of *N* collisions. This object contemplates the correlations between heat exchange, a concept which to the best of our knowledge, has not been explored in the quantum thermodynamics community. We showed that *P*(*Q*1, ... , *QN*) satisfies a fluctuation theorem, which relates this joint distribution with single collision events. This result, we believe, could serve to highlight the interesting prospect of analyzing thermodynamic quantities in time-series and other sequential models.

**Author Contributions:** J.S., A.T. and G.L. equally contributed to conceptualization, investigation and writing of the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** G.L. acknowledges the financial support from the São Paulo funding agency FAPESP and from the Instituto Nacional de Ciência e Tecnologia em Informação Quântica. J.S. would like to acknowledge the financial support from the CAPES/PNPD (process No. 88882.315481/2013-01).

**Acknowledgments:** The authors acknowledge E. Lutz and M. Paternostro for fruitful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Probability Distributions with Singularities**

#### **Federico Corberi 1,2,\* and Alessandro Sarracino <sup>3</sup>**


Received: 28 February 2019; Accepted: 20 March 2019; Published: 21 March 2019

**Abstract:** In this paper we review some general properties of probability distributions which exhibit a singular behavior. After introducing the matter with several examples based on various models of statistical mechanics, we discuss, with the help of such paradigms, the underlying mathematical mechanism producing the singularity and other topics such as the condensation of fluctuations, the relationships with ordinary phase-transitions, the giant response associated to anomalous fluctuations, and the interplay with fluctuation relations.

**Keywords:** large deviations; phase transitions; condensation of fluctuations; fluctuation relations

#### **1. Introduction**

Quantitative predictions on the occurrence of rare events can be very useful particularly when these events can produce macroscopic effects on the system. This occurs, for instance, when a large fluctuation triggers the decay of a metastable state [1] leading the system to a completely different thermodynamic condition. Other examples with rare deviations producing important effects are found in many other contexts, as in information theory [2] and finance [3].

For a collective variable *N*, namely a quantity formed by the addition of many microscopic contributions, such as the energy of a perfect gas or the mass of an aggregate, typical fluctuations are regulated by the central limit theorem. Rare events, instead, may go beyond the theorem's validity and are described by large deviations theory [4,5] which, in principle, aims at describing the whole spectrum of possible fluctuations, no matter how large or rare they are, by means of their full probability distribution *P*(*N*).

It has been found that, in many cases, *P*(*N*) exhibits a singular behavior, in that it is non-differentiable around some value (or values) *Nc* of the fluctuating variable [3,6–39]. Such singularities have an origin akin to those observed in the thermodynamic potentials of systems at criticality. Indeed, a correspondence can be shown between *P*(*N*) and the free energy of a companion system, related to the one under study by a duality map [4,34–36], which is interested by a phase-transition.

Recently, a great effort has been devoted to the characterization of these singular behaviors in the large deviation functions of different models where analytical results can be obtained. This has unveiled a rich phenomenology which shares common features. In most cases non-analyticities are a consequence of a particular condensation phenomenon denoted as condensation of fluctuations.

It occurs when a significant contribution to the fluctuations is built within a limited part of phase-space, or is provided by just one of the degrees of freedom of the system. This is analogous to what happens, for instance, in the usual condensation of a gas when it concentrates in a liquid drop, or in the well-known Bose–Einstein condensation, where the mode with vanishing wavevector contributes macroscopically. However, while usual condensation represents the typical behavior of the system, the condensation of fluctuations can only be observed when certain rare events take place.

Another interesting feature of systems with singular probability distributions can be their extreme sensibility to small perturbations. Usually, the properties of a system made of many constituents or degrees of freedom do not change much if some features of a single particle are slightly changed. This is true both for the average properties and for the fluctuations. For instance, neither the average energy of a gas nor its fluctuations change appreciably if the mass of one single molecule is increased a bit. This is simply because this particle is only one out of an Avogadro number. However, when condensation of fluctuations occurs, one can observe a giant response if the perturbed degree of freedom is exactly the one that contributes macroscopically to the fluctuation.

Singular probability distributions raise the question about the validity of the fluctuation relations (FRs). These relations have been extensively studied recently [40,41] because they reflect general symmetries of the deviations of certain quantities and are believed to contribute to a general understanding of non-equilibrium states. In particular, FRs connect the probability of observing events with a certain value *N* of the fluctuating variable, to the probability of the events associated to the opposite value −*N*. Among other open issues on the subject, one is represented by the case of singular fluctuations. Indeed, the singularity in *Nc* usually separates two regions where fluctuations have very different properties. For instance, on one side of *Nc* one can have a standard situation where all the degrees of freedom contribute, whereas on the other side fluctuations can condense and be determined by the contribution of a single degree. Clearly, if *N* is such that *N* and −*N* fall on different branches of *P*(*N*), namely on the two sides of *Nc*, the mechanism whereby an FR can be fulfilled must be highly non-trivial. In general, singular probability distributions may, or may not, exhibit the FR and a general understanding of this point is still not achieved.

This paper is a brief review devoted to the discussion of singular probability distributions where, without any presumption of neither completeness or mathematical rigor, we present examples of models where such non-analyticities show up, we highlight the mathematical mechanism producing condensation, and we discuss some relevant aspects related to the subject, such as those mentioned above. We do that in a physically oriented spirit, providing whenever possible an intuitive interpretation and a simple perspective. Non-differentiable probability distributions have been previously reviewed also in [42], where however the authors focus on different models and complementary aspects with respect to those addressed in this paper.

The paper is organized as follows. In Section 2 we recall some basic results of probability theory and introduce some notations. In Section 3 we present some models of statistical mechanics where non-differentiable probability distributions have been computed for different collective quantities. In Section 4 we illustrate in detail some phenomena related to the singular distribution function, mainly using the urn model as a paradigm, and discuss how similar behaviors arise in other systems. We also discuss the topic of the fluctuation relations. More specific features, such as giant response and observability, are then presented in Section 5, and, finally, some conclusions are drawn in Section 6.

#### **2. Probability Distributions: Generalities**

We consider a generic stochastic system, whose physical state is defined by the random variable *x* taking values on a suitable phase space. We will be mainly interested in the behavior of collective random variables, that are defined as the sum of a large number of microscopic random variables. For these quantities some general results can be derived [5]. As an example let us consider the sum *N* = ∑*<sup>M</sup> <sup>j</sup>*=<sup>1</sup> *xj* of a sequence of *M* random variables *xj*, with empirical mean

$$\boldsymbol{\rho} = \frac{\boldsymbol{N}}{\boldsymbol{M}} = \frac{1}{\boldsymbol{M}} \sum\_{j=1}^{M} \mathbf{x}\_{j}. \tag{1}$$

The quantities *xj* can represent a sequence of states of a system (for instance, the position of a particle along a trajectory) or an ensemble of variables describing its microscopic constituents (e.g., the energies of the single particles of a gas). In the case of independent identically distributed variables, with expectation *x* and finite variance *σ*, one has that the empirical mean tends to *x* for large *M*, namely

$$\lim\_{M \to \infty} p(\rho - \langle \mathbf{x} \rangle < \epsilon) \to 1,\tag{2}$$

where is a small quantity and hereafter *p*(*E*) (also *P*(*E*) or P(*E*)) is the probability of an event *E*. The above equation represents the Law of Large Numbers.

As a further step, one can describe the statistical behavior of the small fluctuations of *ρ* around the average *ρ*, *δρ* = *ρ* − *ρ*, introducing the quantity

$$z\_M = \frac{1}{\sigma\sqrt{M}}\sum\_{j=1}^{M} (x\_j - \langle \mathbf{x} \rangle),\tag{3}$$

which, for very large *M*, and for *δρ O*(*σ*/ <sup>√</sup>*M*), has the following distribution function

$$p(z\_M = z) \simeq \frac{1}{\sqrt{2\pi}} e^{-\frac{z^2}{2}}.\tag{4}$$

This result is the central limit theorem (CLT), that holds also in the case of weakly correlated variables.

More in general, fluctuations of arbitrary size of the quantity *ρ* can, under certain conditions, be characterized by the large deviation principle (LDP)

$$p(\rho = y) \sim e^{-MI(y)},\tag{5}$$

where *I*(*y*) is the so called rate function. When *p*(*ρ*) has a single absolute maximum (in *ρ*), the rate function is positive everywhere but for *y* = *ρ*, where it vanishes. It is easy to obtain the CLT Equation (4) from the LDP Equation (5) by expanding up to second order the function *I*(*y*) around *ρ*. However, as we will discuss in detail below, there are interesting cases where the LDP in the form Equation (5) is not satisfied.

A simple example where LDP holds and the rate function can be easily computed is obtained by considering {*xj*} as dichotomous variables taking the value +1 with probability *q* and −1 with probability 1 − *q*. Then, using the Stirling approximation, one obtains the explicit expression for the rate function:

$$I(y) = \frac{1+y}{2} \ln \frac{1+y}{2q} + \frac{1-y}{2} \ln \frac{1-y}{2(1-q)}.\tag{6}$$

Expanding Equation (6) around the mean *y* = 2*q* − 1 one has the CLT

$$I(y) \simeq \frac{(y - \langle y \rangle)^2}{2(1 - \langle y \rangle)}.\tag{7}$$

#### **3. Singular Probability Distributions: Examples**

As far as small deviations of a collective variable are considered, the associated probability distribution is usually regular, being a Gaussian when the hypotheses of the CLT are satisfied. Moving to the realm of large deviations, instead, can hold surprises as, for instance, the emergence of non-analyticities. Before deepening the meaning and the bearings of the singular behavior, in this section we first itemize some examples of systems where it has been observed. We will then study it in more detail in some specific models in the following sections.

#### *3.1. Gaussian Model*

The Gaussian model is a reference model of statistical mechanics. An order-parameter field *φ*(*x*) (which in the magnetic language can be thought of as a local magnetization at site *x*) is ruled by the following Hamiltonian

$$\mathcal{H}[\varphi] = \frac{1}{2} \int\_{V} d\vec{x} \left[ (\nabla \varphi)^2 + r\varphi^2(\vec{x}) \right],\tag{8}$$

where *r* > 0 is a parameter and *V* the volume. This simple model can be exactly solved and has a rather trivial phase-diagram without phase transitions.

Let us consider the collective variable

$$N[\varphi] = \int\_{V} d\vec{x} \,\varphi^2(\vec{x}) ,\tag{9}$$

namely the order parameter variance, and its density *ρ* = *N*/*V*. Its probability distribution was computed analytically in [34–36]. The (negative) rate function of this quantity, evaluated in equilibrium at a given temperature *T*, is plotted in Figure 1. The curve has a maximum in correspondence to the most probable value, where *I*(*ρ*) vanishes. Far from such maximum, in the large deviations regime, the rate function exhibits a singularity (marked with a green dot) at *ρ* = *ρc*. In this point the third derivative of the rate function has a discontinuity [34–36]. The existence of such a singularity is related to the fact that, as we will discuss later, fluctuations with *ρ* > *ρ<sup>c</sup>* have a different character with the respect to the ones in the region *ρ* < *ρ<sup>c</sup>* where the average, or typical, behavior of the system (i.e., the most probable value of *ρ*) is located.

**Figure 1.** The (negative) rate function *I*(*ρ*) of the variance *N* of the order parameter field in the Gaussian model in *d* = 3, with *r* = 1, in equilibrium at the temperature *T* = 0.2.

#### *3.2. Large-*N *Model*

Another reference model of statistical mechanics is the description of a magnetic system in terms of the Ginzburg–Landau Hamiltonian

$$\mathcal{H}[\varphi] = \frac{1}{2} \int\_{V} d\vec{x} \left[ (\nabla \varphi)^2 + r\varphi^2 (\vec{x}) + \frac{\mathcal{G}}{2\mathcal{N}} (\varphi^2)^2 \right],\tag{10}$$

where the N -components vectorial field *ϕ* has a meaning similar to that of the Gaussian model, and *r* < 0 and *g* > 0 are parameters. In the large-N limit (sometimes also denoted as the spherical limit) the model is exactly soluble. There is a phase transition at a finite critical temperature *Tc* separating a paramagnetic phase for *T* > *Tc* from a ferromagnetic one at *T* < *Tc*.

The probability distribution of the energy *N*(*t*, *tw*) = H[*ϕ*, *t*] − H[*ϕ*, *tw*] exchanged by the system in a time interval [*tw*, *t*] with a thermal bath was computed exactly in [37]. The (negative) rate function of the intensive quantity *ρ*(*t*, *tw*) = *N*(*t*, *tw*)/*V* is shown in Figure 2. This figure refers to the case of a system quenched from a very high temperature to another *T* < *Tc*. Also in this case there is a singularity corresponding to a certain value the quantity *ρ*(*t*, *tw*) = *ρ<sup>c</sup>* where the third derivative has a discontinuity, and this reflects a different mechanism of heat exchanges for *ρ* < *ρ<sup>c</sup>* and for *ρ* > *ρc*.

**Figure 2.** The (negative) rate function *I*(*ρ*) of the probability distribution *P*(*N*) of the energy *N* exchanged by the large-N model in *d* = 3, with *g* = −*r* = 1, with the environment after a quench to zero temperature.

#### *3.3. Urn Model*

Let us consider a set of integer variables *ni* ≥ 0 (*i* = 1, ... , *M*) equally distributed in such a way that the probability of having a certain value *n* of *ni* is

$$p(n) = \mathbb{Z}^{-1}(n+1)^{-k},\tag{11}$$

where *ζ* is a normalization constant and *k* a parameter. One can think of having *M* urns, each of them hosts a quantity *nm* of particles taken with probability Equation (11) from a reservoir. This setting is appropriate to describe a wealth of situations in many areas of science, from network dynamics to financial data. The probability distribution of the total number of particles

$$N = \sum\_{m=1}^{M} n\_m \tag{12}$$

was studied for large *M* in different contexts [14,17,21–23,43]. The (negative) rate function is shown in Figure 3. Also in this model it is found that, if *k* > 2, there is a singularity at *ρ* = *ρc*, that in this particular case coincides with the average value *ρ*. Notice that in this case, at variance with the previous examples, the rate function vanishes in the whole region *ρ* ≥ *ρc*. This is due to the fact that *P*(*ρ*) has a weaker dependence on *M* with respect to the exponential one of Equation (5), and hence the LDP is violated for *ρ* > *ρc*. We will comment later on that.

**Figure 3.** The rate function *I*(*ρ*) of the probability distribution *P*(*N*) of the total number of particles *N* in the urn model with *k* = 3.

#### *3.4. Stochastic Maxwell-Lorentz Particle Model*

The so-called stochastic Maxwell–Lorentz gas [44,45] consists of a probe particle of mass *m* whose velocity *v* changes due to the collisions with bath particles, of mass *M* at temperature *T*, and due to the acceleration produced by an external force field E. Collisions with the scatterers change instantaneously the probe's velocity from *v* to *v* and we assume the simple collision rule *v* = *V*, where *V* is the velocity of the scatterer, drawn from a Gaussian distribution:

$$P\_{\text{scatt}}(V) = \sqrt{\frac{\mathcal{M}}{2\pi T}} e^{-\frac{MV^2}{2T}}.\tag{13}$$

The scatterers play the role of a thermal bath in contact with the probe particle. This model is a particular case of a more general class of systems studied in [44,45]. During a time *τ* between two consecutive collisions, the probe performs a deterministic motion under the action of the field E. We assume that the duration of flight times *τ* is exponentially distributed *Pτ*(*τ*) = <sup>1</sup> *<sup>τ</sup><sup>c</sup>* exp(−*τ*/*τc*) and independent of the relative velocity of the particles. The system reaches a non-equilibrium stationary state characterized by a total entropy production Δ*s*tot, associated with the velocity *v*(*t*), defined as

$$\Delta s\_{\text{tot}}(t) = \ln \frac{P(\{\upsilon(s)\}\_0^t)}{P(\{\overline{\upsilon(s)}\}\_0^t)} \, \tag{14}$$

where *<sup>P</sup>*({*v*(*s*)}*<sup>t</sup>* <sup>0</sup>) and *<sup>P</sup>*({*v*(*s*)}*<sup>t</sup>* <sup>0</sup>) are, respectively, the pdf of a path {*v*(*s*)}*<sup>t</sup>* <sup>0</sup> spanning the time interval [0, *<sup>t</sup>*] and of the time-reversed path {*v*(*s*)}*<sup>t</sup>* <sup>0</sup> = {−*v*(*<sup>t</sup>* − *<sup>s</sup>*)}*<sup>t</sup>* <sup>0</sup> [46]. This fluctuating quantity takes contributions at any time and is therefore extensive in *t*. In this example it plays the role of the collective variable *N*, and *t* plays the role of the number *M* of elements contributing to it.

The rate function *I*(*ρ*) of the quantity *ρ* = Δ*s*tot/*t* was studied in [11] by means of numerical simulations for finite times and analytically in the limit *t* → ∞. This quantity is shown in Figure 4, where *<sup>ρ</sup><sup>c</sup>* = *<sup>m</sup>τc*E2/*θ*, with *<sup>θ</sup>* = *Tm*/*<sup>M</sup>* playing the role of an effective temperature [47]. Also in this case, as for the urn model, *I*(*ρ*) vanishes and *P*(Δ*stot*) does not satisfy a standard LDP for *ρ* > *ρc*. Indeed it can be shown that the far positive tail of *<sup>P</sup>*(Δ*stot*) scales exponentially with <sup>√</sup>*<sup>t</sup>* rather than with *t* [11], how it should be if the LDP (Equation (5)) holds. Recently, the nature of the singularities in *I*(*ρ*) and their physical meaning have been thoroughly discussed in a similar model in [48], where the observed non-analytical behaviors have been related to a first-order dynamical phase transition.

**Figure 4.** The rate function *I*(*ρ*) of the quantity *ρ* = Δ*s*tot/*t* for the Maxwell–Lorentz gas model [11], computed analytically in the limit *t* → ∞.

#### *3.5. Some Other Models*

We have discussed above some models where a singular probability distribution was found. All these cases can be grouped into two classes: the first contains the cases where the rate function is well defined, although it contains some non-analyticity point. The examples of Sections 3.1 and 3.2 behave in this way. The second class is the one represented by the urn model, where the probability distribution is still singular, but the the rate function is not defined in a certain region (that is to say it vanishes identically). The Maxwell–Lorentz gas is an example where the two behaviors are exhibited in different regions of the fluctuation spectrum.

Beyond the cases discussed before, other examples of singular behavior include the probability distribution of the work done by active particles [38], of the heat exchanged by harmonic oscillators during a quench with a thermal bath [39], of the magnetization in the spherical model [6,7], of the displacement of a Brownian walker with memory [10], of the work done in a quantum quench [12], and many others [4,25–33].

We also mention the case where the singularity appears as a "kink" at zero in the probability distribution, showing a linear regime for negative values. This behavior has been observed in the distribution of the entropy production and of other currents for a driven particle in periodic potentials [49–51], in a molecular motor model, described in [52], and in the experimental results reported in [53], where the large deviation function of the velocity of a granular rod was measured. In general, the presence of the kink can be related to different physical mechanisms [54], such as intermittency [55], detailed fluctuation theorem [56], and dynamical phase transitions [57].

#### **4. General Features of Singular Probability Distributions**

In this section we will discuss some general properties of singular probability distributions observed in the different models mentioned above, focusing on the common physical interpretation and on the underlying mathematical structure.

#### *4.1. Duality*

The singular behavior of the probability distribution seen in the examples of the previous section has an interpretation akin to the occurrence of phase transitions in ordinary critical phenomena. In order to discuss this point we can refer to the Gaussian model as a paradigm. The partition function is

*Z* = *δϕ* P(*ϕ*), (15)

where P is the probability of microscopic configurations as specified by the field *ϕ*. For instance, in a canonical setting it is P(*ϕ*) = *<sup>Z</sup>*−<sup>1</sup> exp[−*β*H(*ϕ*)], where *<sup>β</sup>* is the inverse temperature *<sup>β</sup>* = 1/(*kBT*); in this case *Z* depends on *T* and *V*, the volume. On the other hand the probability of the collective variable *N* of Equation (9) can be written as

$$P(N) = \int \delta\varphi \,\mathcal{P}(\varphi) \,\delta\left(\int\_V d\vec{x} \varphi^2(\vec{x}) - N\right). \tag{16}$$

In view of Equation (15), one can recognize Equation (16) as a partition function as well. However this is not the partition function of the original model that is, in this example, the Gaussian one. Instead, *P*(*N*) in Equation (16) can be interpreted as the partition function of a dual system that can be obtained from the original one upon removing all the configurations such that the argument of the delta function in Equation (16) does not vanish. In other words, this is the model one arrives at upon constraining configurations in a certain way. In this case the requirement is that the variance of *ϕ* must equal a given value *N*. Such a system, a Gaussian model with a constraint on the variance, is the spherical model of Berlin and Kac [58].

The equilibrium properties of the spherical model are exactly known. For fixed *N*, there is a phase-transition at a critical temperature *Tc*, from a disordered phase for *T* ≥ *Tc* to an ordered one below *Tc*. Equivalently, still in the Berlin–Kac model, if one keeps *T* fixed, the transition occurs changing the variance *N*[*ϕ*] defined in Equation (9) upon crossing a critical value *Nc*. The ordered phase is found for *N* > *Nc*, in this case. The presence of such a phase transition crossing *Nc* determines a singularity of the partition function *P*(*N*) of the spherical model (Equation (16)) at *N* = *Nc*. However the same quantity *P*(*N*) is also the probability distribution of the quantity *N*[*ϕ*] in the context originally considered, the Gaussian model. This explains what one observes in Figure 1. *Nc* is the value of *N* marked by a dot in this figure, where the singular behavior shows up.

This dual interpretation of *P*(*N*), as a probability distribution of a collective variable in the original model, or as a partition function in a dual model, may help to understand why singularities are manifested in the probability distributions. Indeed, if one asks the question: why a simple model without phase transitions, such as the Gaussian model, exhibits a non trivial singularity in the probability distribution *P*(*N*), the answer can be that, although the original model is quite simple, the dual one is far from being trivial, with a phase-transition induced by the presence of the constraint. This generates anomalous behavior in the fluctuation spectrum of the original model.

We have discussed the fact that imposing a constraint to the Gaussian model we change the system into a dual one that is interested by a phase transition, since this is the spherical model. Is this an isolated example or has this feature some generality? The answer is that it occurs quite often. Besides the above mentioned spherical model, another well known example where the same mechanism is at work is the perfect boson gas. There is no phase transition in a gas with a non conserved number *N* of bosons, as in the case of photons, but if the number *N* of particles is fixed Bose–Einstein condensation happens. The partition function of the conserved bosons, for a given volume and temperature, has a singularity at a certain value of the boson number *N* = *Nc* (or density). This singularity corresponds to the critical number of particles below which the condensed phase develops. According to the duality principle discussed above, this implies that the probability distribution of the number of bosons in a system of, say, photons, where this number is allowed to fluctuate, will be singular at the same value *Nc* of the random variable *N* [33]. The very urn model is another instructive example. One can consider a model, dual to the one discussed in Section 3.3, where the total number of balls is conserved [21]. Marbles can only be exchanged among boxes and their density *ρ* is an external control parameter. This model is known to be interested, for *k* > 2, by a phase transition crossing *ρ* = *ρc*. Notice that, since *ρ* is a control parameter, having *ρ* > *ρ<sup>c</sup>* in this dual model is not a rare event (as in the model introduced in Section 3.3). A similar situation is found in related models such as the zero range process [18,21,28].

#### *4.2. Condensation*

In order to see how singularities may come about in another perspective we will discuss the phenomenon in the framework of the urn models, where the physical meaning is probably more transparent in term of a condensation mechanism. Something similar occurs also in the other models considered in Section 3, regardless of the fact that the rate function is well defined or not.

Let us consider the conditional probability *π*(*n*, *N*, *M*) that one of the *M* a priori equivalent urns contains *n* particles, given that there are *N* particles in the whole system. This quantity can be evaluated exactly and is shown in Figure 5 (normalized by its value in *n* = 0 to better compare curves with different *N* in a single figure). Let us discuss its properties. First of all *π* vanishes for *n* > *N*, since it is impossible that an urn contains more particles than the whole system. Secondly, for small *n* one has *π*(*n*, *N*, *M*) ∝ *p*(*n*) (dotted green line in Figure 5). This means that, as far as very few particles are stored in the tagged urn, the condition on the total number of balls is irrelevant.

**Figure 5.** The function *π*(*n*, *N*, *M*) is plotted, for *k* = 3 and *M* = 100, against *n* + 1 for two values of *N*: *N* = 35, corresponding to a case without condensations, and *N* = 300, corresponding to a condensed situation. The dotted green curve is the power-law *x*−*k*.

More interestingly, at large *n*, *n N*, different behaviors are observed in the region of relatively small *N*, and in the one with relatively large *N*, exemplified by *N* = 35 and *N* = 300, respectively, in Figure 5. In the former case *π* is exponentially damped at large *n*, meaning that accommodating many particles in a single urn is probabilistically very unfavorable. In the latter case there is a peak at a value of *n* or of order *N*. This means that a significant fraction of the total number of particles is located in a single urn. This is the condensation phenomenon. (We will see in the next section that in this particular model occurs when *k* > 2 and for sufficiently large densities). The essence of a condensation phenomenon is that a given quantity is not fairly distributed among many degrees of freedom, but is concentrated in a single one. This is particularly clear in the urn model, where one particular urn contains a macroscopic fraction of balls.

One easily realizes that something similar occurs in the other example models discussed in Section 3. For instance, in the model of Section 3.1, writing *N*[*ϕ*] in terms of the Fourier components *ϕ<sup>k</sup>* of the field *ϕ* as

$$N[\varphi] = \frac{1}{V} \sum\_{\vec{k}} |\varphi\_{\vec{k}}|^2 \,. \tag{17}$$

one can show that while for *N* ≤ *Nc* (or equivalently *ρ* ≤ *ρc*) all the Fourier components add up to realize the sum in Equation (17) in a comparable way, for *N* > *Nc* the term with *k* = 0 alone provides the most important contribution to the sum. A similar mechanism, with the dominance of the *k* = 0 term, is also at work in the example of Section 3.2. In the Maxwell–Lorentz particle model (Section 3.4) one has that normal entropy fluctuations are formed by the addition of many contributions associated to many short flights of the probe particle. However above the critical threshold *ρ<sup>c</sup>* they are associated to a single event which is responsible for a macroscopic contribution to the entropy production. This event is a long flight of the probe particle with no collisions with the scatterers. For more details and a very accurate analytical description of these kinds of behaviors in a similar system re-framed in the context of active particles, see the recent work [48].

#### *4.3. Mathematical Mechanism*

In the previous section we have discussed the phenomenon of condensation on physical grounds. In this section we show the underlying mathematical mechanism. We will give a description as simple as possible, without presumption of mathematical rigor, in the framework of the urn model.

The probability distribution of the total number of particles *N* reads

$$P(N,M) = \sum\_{n\_1, n\_2, \dots, n\_M} p(n\_1)p(n\_2)\dots p(n\_M) \,\delta\_{\sum\_{m=1}^M n\_m N} \,\,\,\,\tag{18}$$

where *δa*,*<sup>b</sup>* is the Kronecker function and in the leftmost sum the variables *n*1, *n*2, ... , *nM* run from 0 to ∞. Using the representation

$$\delta\_{a,b} = \frac{1}{2\pi i} \oint dz \, z^{-(b-a+1)}\tag{19}$$

of the *δ* function one arrives at

$$P(N,M) = \frac{1}{2\pi i} \oint dz \, e^{M[\ln Q(z) - \rho \ln z]},\tag{20}$$

where

$$Q(z) = \sum\_{n=0}^{\infty} p(n)z^n,\tag{21}$$

and we have confused *<sup>N</sup>*+<sup>1</sup> *<sup>M</sup>* with *<sup>ρ</sup>* <sup>=</sup> *<sup>N</sup>*/*<sup>M</sup>* for large *<sup>M</sup>*. Still for large *<sup>M</sup>*, the integral in Equation (20) can be evaluated by the steepest descent method as

$$P(N,M) \simeq \mathfrak{e}^{-MR(\rho)},\tag{22}$$

where

$$R(\rho) = -\ln Q[z^\*(\rho)] + \rho \ln z^\*(\rho),\tag{23}$$

with *z*∗ the value of *z* for which the argument in the exponential of Equation (20) has its maximum value. This in turn is given by the following implicit saddle-point equation

$$
\Theta(z^\*) = \rho,\tag{24}
$$

where

$$
\Theta(z^\*) = z^\* \frac{Q'(z^\*)}{Q(z^\*)}.\tag{25}
$$

Let us study this equation. Clearly, it must be *z* ≤ 1 in order for the sums hidden in *Q* and *Q* to converge. It can also be easily seen that Θ(0) = 0 and that this function increases with *z* up to

$$\Theta(1) = \begin{cases} \infty, & k \le 2 \\ \Theta\_{M\nu} & k > 2, \end{cases} \tag{26}$$

where Θ*<sup>M</sup>* is a finite positive number. The function Θ(*z*) is shown in Figure 6, for two values of the parameter *k*. As it is clear from this figure, for *k* > 2 the saddle point Equation (24) admits a solution only for 0 ≤ *ρ* ≤ *ρ<sup>c</sup>* = Θ(1). It is trivial to show that *ρ<sup>c</sup>* ≡ *ρ* = ∑*<sup>n</sup> np*(*n*). However nothing prevents fluctuations with *ρ* > *ρ* to occur. How can we recover the model solution for *ρ* > *ρ*? We know that for such high densities urns are no longer equivalent: there is one—say the first—which hosts an extensive number of particles and condensation occurs. In a physically oriented approach, we can take into account this fact by writing, in place of Equation (18), the following

$$P(N,M) = M \sum\_{n\_1=0}^{\infty} p(n\_1) \sum\_{n\_2, n\_3, \dots, n\_M} p(n\_2) p(n\_3) \dots p(n\_M) \,\delta\_{\sum\_{m=2}^{M} n\_m N - n\_1}.\tag{27}$$

The factor *M* in front of the r.h.s. stems from the fact that there are *M* ways to chose the urn (denoted as 1) to be singled out. Repeating the mathematical manipulations as in Equations (18) and (20), but only on the sum ∑*n*2,*n*3,...,*nM* . . . , one arrives at

$$P(N,M) = \frac{M}{2\pi i} \sum\_{n\_1=0}^{\infty} p(n\_1) \oint dz \, e^{M[\ln Q(z) - (\rho - \frac{n\_1}{M})\ln z]}.\tag{28}$$

Evaluating the integral with the steepest descent method, the saddle point equation is now

$$
\Theta(z^\*) = \rho - \frac{n\_1}{M}.\tag{29}
$$

Notice that in a normal situation, where condensation does not occur, in the thermodynamic limit where *M* → ∞ with fixed *ρ*, the typical number of particles in a single urn does not depend on the number of urns. Therefore the last term on the r.h.s. of Equation (29) is negligible and one goes back to the previous saddle point Equation (24). However, when condensation occurs (i.e., with *k* > 2 and *ρ* > *ρ*) the only possibility to close the model equations is to have the last term in Equation (29) finite. In conclusion one has

$$\begin{cases} \; z^\* < 1, & \frac{n\_1}{M} \simeq 0 \\\; z^\* = 1, & \frac{n\_1}{M} = \rho - \langle \rho \rangle \qquad \text{condensation}. \end{cases} \tag{30}$$

**Figure 6.** The function Θ(*z*) is shown for *k* = 1.5 and *k* = 2.5.

Clearly we are in the presence of a phase-transition resembling the ferro-paramagnetic or the gas–liquid transitions. There are two phases with qualitatively different behaviors. However, at variance with usual phase transitions, here the parameter producing the transition is not an external one that can be varied at will, but the value of the spontaneously fluctuating variable *N*. Another difference with usual phase transitions is the fact that here there is no interaction among urns. Despite that, urns are not completely independent due to the constraint over the number of particles

represented by the Kronecker function in Equations (18) and (27). This constraint can be regarded as an effective interaction determining the transition (it can be easily seen, in fact, that without such conservation there is no transition).

Notice that it is *n*1/*M* = 0 in the normal phase and *n*1/*M* = 0 in the condensed phase, therefore this quantity represents the order parameter of the transition. Despite the fact that a priori the system (i.e., the Hamiltonian) is invariant under a permutation of the urns, namely all boxes are equal, this property is not shared by the physical realization of the actual state of the system when condensation occurs, since one urn behaves very differently from the others. We are in the presence of spontaneous symmetry breaking.

As a final remark, let us note that the phenomenon of condensation in the sum of many identically distributed variables is not specific to an algebraic decay of *p*(*n*), or to the discrete value of the variable *n*. Indeed it is found [31] that it occurs provided that ∑*<sup>n</sup> np*(*n*) < ∞. Condensation in the presence of a stretched exponential *p*(*n*), for instance, has been discussed in [59,60]. Finally, we mention the fact that in the context of Lévy flights the phenomenon of condensation is usually referred to as the big jump principle [61].

#### *4.4. Fluctuation Relation*

The Fluctuation Relation is one of the few general results of non-equilibrium statistical mechanics, expressing an asymmetry property of the fluctuations of some extensive (in time or in number of degrees of freedom) quantities *N* [40]. The FR reads

$$\frac{P(N/M=\rho)}{P(N/M=-\rho)} = e^{cM\rho + o(M)},\tag{31}$$

where *c* is a constant, and *o*(*M*) stands for sub-linear corrections in *M*. Usually, the exponential form of the FR is related to two properties of *P*(*N*/*M*): (i) it satisfies a LDP Equation (5), and (ii) the rate function *I*(*ρ*) has the symmetry:

$$I(-\rho) - I(\rho) = \mathfrak{c}\rho.\tag{32}$$

These two conditions, with *I*(*ρ*) different from 0 and ∞, are known to be sufficient for *N*/*M* to satisfy a FR (see, e.g., [4] and references therein).

It is interesting to consider the validity of an FR in the case of probability distributions with singularities. First, let us note that, when the singularity appears in zero, as in the case of the "kink" mentioned in Section 3.5, then the validity of an FR is clearly not affected by the singularity. More in general, the FR can also be satisfied by a pdf for which a standard (namely, with a leading exponential scaling in *M*) LDP does not hold. This can be observed for instance in the driven Maxwell–Lorentz gas described in Section 3.4. In this model it has been shown [11] that the entropy production calculated over a time *t* satisfies an FR, even though the far positive tail of its pdf scales exponentially with <sup>√</sup>*<sup>t</sup>* rather than *<sup>t</sup>*. In this case the validity of the FR can be exploited to extract some information on the behavior of the probability distribution in the regions where the stretched-exponential scaling takes place.

The FR Equation (32) in the presence of a singular rate function has been also observed [37], besides the already mentioned Maxwell-Lorentz case , in some large time limit for the exchanged heat, in the large-N model of Section 3.2. More recently, it has been shown [39] that the rate function of the heat exchanged by a set of uncoupled Brownian oscillators with the thermostat during a non-stationary relaxation process does not satisfy an FR in the form Equation (31). Although, even in this case, the rate function shows a singular behavior in the limit of a large number of degrees of freedom, the lack of a standard FR is not necessarily related to the presence of the singularity.

#### **5. Some Peculiarities of Singular Distributions**

#### *5.1. Giant Response*

Generally, the behavior of a collective quantity such as the empirical mean Equation (1) is not substantially altered if, for large *M*, the properties of only one out of *M* variables is slightly modified. For instance, one does not expect to observe any significant change in the thermodynamic properties of a gas of identical molecules if one is replaced with another of a different substance. This is because the collective properties are determined by the synergic contribution of a huge number *M* of constituents, and hence the features of a single molecule are negligible. This is true not only for the typical properties but also for the fluctuation distribution. However, the situation can be dramatically different when singular probability distributions enter the game.

Let us show this with the prototypical example of the urn model. We consider a slightly modified version of the model defined in Section 3.3, where a single variable, say *n*, is distributed as in Equation (11) but with an exponent *k* that may be different from the one, *k*, of all the remaining ones. Let us now look at Equations (28)–(30). In a situation where condensation does not occur, as we remarked earlier, the effect of a single variable is negligible, the first line of Equation (30) applies and hence *nm <sup>M</sup>* 0, for any *m*. On the other hand, in the presence of condensation, the second line of Equation (30) holds. In the case of equally distributed variables condensation occurs with equal probability in any of the urns. However, if the -th variable behaves differently, one has to understand if the condensing variable could be the -th, or any of the remaining ones. Both the cases can occur, depending on the values of the exponents *k* and *k*.

For *k* > *k* > 2 (the latter inequivalence being needed for condensation) it is *p*(*n* = *n*) *p*(*nm* = *n*) for large *n* (with = *m*). Hence the condensation phenomenon, which occurs by letting a huge amount of particles occupy a single urn, is unfavoured in the -th urn. The situation in this case is analogous to the one discussed before with equally distributed variables, i.e., with *k* = *k*. However for *k* > *k* > 2 the opposite occurs, the condensing variable is the -th. Hence Equation (28) applies with *n*<sup>1</sup> replaced by *n*. One sees from Equation (28) that, when condensation occurs, *P*(*N*, *M*) is proportional to *p*(*n*). Since *k* = *k*, *P*(*N*, *M*) turns out to be different from the one found for equally distributed variables. Hence, in this case, an even small change of the properties of a single variable can trigger the form of the probability distribution of the collective variable *N*, a fact that is sometimes referred to as giant response.

This is illustrated in Figure 7. Here *P*(*N*, *M*) is compared for three different choices of the exponents *k*, *k*. The continuous blue curve with asterisks refers to the case (i) with identically distributed variables with *k* = *k* = 3. Similarly, the dot-dashed green curve with squares corresponds to the situation with (ii) *k* = *k* = 6. Instead, the dashed-magenta curve with circles corresponds to non-identically distributed variables with (iii) *k* = 3 and *k* = 6. Notice that in the region to the left of the maximum, where condensation does not occur (because in this region *ρ* < *ρ*), the curves of the cases (ii) and (iii) coincide. This nicely shows that in the absence of condensation the shift of the properties of a single variable does not influence the collective behavior of the system. For *ρ* > *ρ*, on the other hand, the form of *P* drastically changes in going from (ii) to (iii), namely by perturbing the properties of one single variable. Even more impressive, the slope of the curve for case (iii) is the same as that of case (i), showing that this feature is dictated by the sole properties of the variable, *n*, which in case (iii) behaves as in (i).

**Figure 7.** *P* is plotted for *M* = 333 and the three different choices (see text) (i) *k* ≡ *k* = 3, continuous blue with asterisks, (ii) *k* ≡ *k* = 6, dot-dashed green with squares and (iii) *k* = 3 , *k* = 6, dashed magenta with a circles.

#### *5.2. Development of a Singular Fluctuation*

We have seen in Section 4.1 that a singularity in the probability distribution can be interpreted as a phase transition occurring at a critical value of *ρ*, playing the role of a control parameter. The analogy can be pushed a step further. When a system is prepared in a certain equilibrium state and then a control parameter is changed as to make it cross a phase transition, the ensuing dynamics can be slow and characterized by a dynamical scaling symmetry associated most of the times with an ever growing length scale [62–65]. Typical examples are magnets and binary systems quenched across the critical temperature, and glassy systems.

Building on the analogy above, one might expect something similar to happen if one prepares a system with a singular *P*(*N*) in a state such that the fluctuating collective variable *N* takes a definite value *N*<sup>0</sup> on one side (say the left) of the critical value *Nc* where the singularity takes place. If the system is then left to evolve freely, all possible fluctuations will take place, including those on the other side (say the right) of the singularity. Due to the duality principle, this process should occur in a way akin to the kinetics of a system brought across a phase-transition. Hence slow evolution and dynamical scaling should be observed. This has actually been shown to be the case, as we discuss below.

Upon supplementing the urn model of Section 3.3 with a kinetic rule allowing the system to exchange single particles with an external reservoir in such a way that the stationary occupation probability of any urn is given by Equation (11), one can solve exactly [43] the evolution of a system whose initial state is such that condensation is not present. In the following we will discuss the case in which the initial value of the density is *ρ* = *ρ*. Starting from this configuration, corresponding to an initial form *P*(*N*, *M*, *t* = 0) of the probability distribution of the collective variable, the system will evolve as to produce all the allowed fluctuations. Hence *P*(*N*, *M*, *t*) becomes time-dependent. Clearly, for long times it is expected to approach the stationary value *P*(*N*), with the singular behavior already discussed. This curve is plotted in Figure 8, with a dotted green line.

In this figure one sees that the time evolution of the probability *P*(*N*, *M*, *t*) towards this asymptotic form is much different on the two sides of the singularity. For *N* < *N*, in the normal region without condensation, the evolution is fast and the asymptotic form of the probability is attained at relatively short times. Indeed, already the red curve for *<sup>t</sup>* = 1.2×106 is indistinguishable from the stationary form and increasing time does not change anything. Conversely, the evolution is slow in the condensing region for *N* > *Nc*. Here one sees that, at any time, the asymptotic form is only attained up to a value *N* = *ν*(*t*), beyond which *P*(*N*, *M*, *T*) drops much faster than what expected asymptotically. It can be shown that *ν*(*t*) grows indefinitely in an algebraic way, much in the same way as a characteristic

growing length does in systems quenched across a phase transition. In addition, a dynamical scaling symmetry can be shown to be at work also in this case. The origin of this slow kinetics is clearly due to the difficulty to condense a huge amount of particles in a single urn by exchanging single particles with the reservoir.

**Figure 8.** The probability *P*(*N*, *M*, *t*) with *k* = 3 is plotted against *N* with double logarithmic scales for different times (see key), exponentially spaced. The dotted green line is the asymptotic form.

#### *5.3. Observability*

In the previous sections we discussed some peculiar properties of singular distribution functions. A natural question is if such features can be observed in practical situations. Indeed, the non-analycities of the probability distributions are observed in the regime of large deviations, namely outside the range of small fluctuations which are generally described by the central limit theorem and are more likely to be observed.

To make more clear this point let us make reference to the Gaussian model and, specifically, to Figure 1. In this case, in order to detect singular deviations, *ρ* = *ρ<sup>c</sup>* must be exceeded. Namely, the system has to move quite far from the most likely observed value—the maximum of the distribution. If the LDP Equation (5) holds (it does so in this model) the possibility to observe such a large fluctuation is extremely small already for moderately large values of the number of constituents *M* (or volume *V*), due to the exponential damping in *M* expressed by Equation (5). But the situation is different if the LDP is violated. This occurs, for instance, in the urn model or in the Maxwell–Lorentz gas, in the fluctuation range where the rate function vanishes. In the former model one can easily check from Equations (28)–(30) that the LDP is obeyed in the non-condensing regime but it is violated when condensation occurs. In fact, it is trivial to see that with *z*∗ = 1 the saddle point evaluation of the integral in Equation (28) gives an exponential with an argument that is identically vanishing. As a consequence fluctuations away from the average are no longer damped exponentially in *M*, but only as *M*1−*<sup>k</sup>* (keeping *ρ* fixed). This is why the rate function of the model vanishes in the whole sector *ρ* > *ρ<sup>c</sup>* where condensation occurs (see Figure 3), despite the fact that *P*(*N*, *M*) decays for *ρ* > *ρ*, as it can be seen in Figure 7. Due to this much softer decay, there is a better chance to observe singular fluctuations in this model than in others, e.g., the Gaussian model, where the LDP holds. A similar situation, with LDP violations, is observed also in the Maxwell–Lorentz particle model (for *ρ* = Δ*stot*/*t* > *ρc*) [11] and in Bose–Einstein condensates [33].

#### **6. Summary and Conclusions**

In this paper we have shortly reviewed the issue of probability distributions characterized by non-analyticities. Naively, this feature could be considered as a rare manifestation of curious mathematical pathologies occurring in scholarly model with uncertain relations to the physical world. In reality, singular probability distributions have been shown analytically to occur in very simple and fundamental models of statistical mechanics, such as the Gaussian one, and not only in weird non-equilibrium states but also in equilibrium. Furthermore, they have been detected in numerical simulations and, most importantly, also in real experiments. This widespread occurrence points towards an underlying general mechanism for the development of singularities in the fluctuation probability. This paper has been conceived in order to highlight and discuss, at a simple and physically oriented level, at least some of such general features.

In the first part of the paper, after recalling basic and general concepts of probability theory we have reviewed some models where singular fluctuation spectra have been observed. These range from the aforementioned Gaussian model to the spherical limit of a ferromagnet, from the so-called urn model to a description of the Maxwell–Lorentz gas. In all these cases the deviations of certain collective observables are described by non-analytical probability distributions, which, in the case when LDP holds, are characterized by the presence of exponential branches.

The non-analytical behavior has been interpreted as due to the same mechanism whereby singularities develop in the thermodynamic functions of systems experiencing phase transitions. Indeed we have discussed the fact that a singular fluctuation distribution function can be mapped onto a thermodynamic potential of a dual model with a critical point. The singularity appears similarly to what one observes in thermodynamic functions when a condensation transition is present. When such feature occurs at the level of fluctuations, at variance with the usual examples of condensation, one speaks of condensation of fluctuations.

Singularities of the probability distributions can have a scarce practical relevance if they occur in regions where fluctuations have a negligible chance to be observed. However, in some of the cases considered in this paper the non-analytical behavior is associated to the breakdown of the large deviation principle. As a result, large fluctuations of macrovariables have a better chance to be observed even in systems with a relatively large number of degrees of freedom. In this case the presence of singularities not only can be observed, but its effects can be appreciated. Perhaps, one of the most intriguing one is the so called giant susceptibility, whereby slightly tuning the properties of even one single component, say a molecule of a gas, can have catastrophic consequences on the behavior of the whole system.

Non-analyticity points in the probability distributions also influence the way in which rare fluctuations are developed out of typical state where they are absent. Indeed, it has been shown that large fluctuations in the region where condensation occurs are formed by means of a complex slow dynamics which resembles, once again a manifestation of a dual behavior, that of systems brought across a phase transition. The knowledge of the dynamical path leading to a rare fluctuation may have important consequences in those cases when such deviations lead to catastrophic events, as in the case of extinctions or bankruptcies.

Among the several perspectives of future studies in this context, we mention the possibility to explore the role of correlated noise on the large deviations, for instance in models of active particles where some analytical results can be obtained [66]; the meaning of singularities, which are related to non-equilibrium phase transitions, within the general framework of the macroscopic fluctuation theory [67]; the relation between the presence of singularities and the validity of the fluctuation relation for entropy production or related quantities in more general cases; the role of correlations among random variables in the anomalous large deviations, as observed for instance in conditioned random walks [68] and Brownian motion [69]; the effect of inhomogeneous rates in bulk-driven exclusion processes [70].

**Author Contributions:** F.C. and A.S. equally contributed to conceptualization, investigation and writing of the paper.

**Funding:** F.C. acknowledges funding from grant PRIN 2015K7KK8L covering also the publication costs of the present open access article. A.S. acknowledges support from "Programma VALERE" of University of Campania "L. Vanvitelli".

**Acknowledgments:** F.C. and A.S. acknowledge A. Crisanti, L.F. Cugliandolo, G. Gonnella, G. Gradenigo, A. Piscitelli, A. Puglisi, H. Touchette, A. Vulpiani and M. Zannetti for a long collaboration on some of the issues here discussed.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
