How has Parisi's nobel-prize winning work been applied to all kinds of complex systems? As discussed briefly in this APS Physics editorial, Nobel Prize: Complexity, from Atoms to Atmospheres, the most important works of the recent Nobel prize winning physicist concern the study of the statistical physics of spin glasses. Taking a look at his original papers on the subject (linked to from the editorial), it seems to me that his revolutionary core insight was that order parameters for such systems must be functions instead of scalar parameters.
Now, it has been often said, that his results have found widespread application in all kinds of systems, up to 'flocks of birds'.
Is there a model-free explanation of the technique, making it part of the standard 'toolkit' for complex systems? And what kind of problems does it address, specifically?
 A: The Forbes article John called our attention to summarizes it in the following way:

When you have a large number of things interacting with each other, and an element of disorder, Parisi’s replica symmetry breaking technique can often be used to predict aspects of the collective behavior.

This is from the article's closing paragraph. The rest of it explains "disorder" (such as that of a glass, in contrast to the regular pattern of a crystal), calling attention to the fact that its configurations, though random, are weirdly persistent; and gives an idea of the "replica trick" (which roughly is to do the calculations not for one given configuration, but for ensembles of copies of the system with the same basic structure but in different states); and, last, tells that the key contribution from Parisi is to clarify the meaning and consequences of this trick — such as that configurations of similar energy might be connected only by an energetically expensive path in configuration space (which explains why a given random configuration can be persistent).
And then you see why this can find wide application — interacting systems with some degree of disorder are ubiquitous and a number of them can benefit from the replica trick approach, including: "structural glasses, granular systems, protein configurations, financial markets, machine-learning algorithms", and "neural and communications networks, granular materials, and flocks of birds."
A: Here is the model-free explanation: statistical mechanics concerns systems of many (N>>1) interacting elements. Observables can be computed from the partition function $Z$, possibly with added fields. In principle this requires some microscopic probability distribution over dynamical elements, although maximum entropy approaches are also widely used. If we have universality, then not all details of this distribution will be important.
More precisely, observables are obtained by derivatives of $\log Z$. In the thermodynamic limit, the finite observables are of the form
$\mathcal{O} = \lim_{N \to \infty} \frac{1}{N}\frac{\partial \log Z}{\partial h}$
where $h$ is some field-like parameter.
In the simplest class of systems, those we call equilibrium, all but a finite number of degrees of freedom are dynamical: we average over them to obtain expected values of observables. In these systems we compute $Z$, possibly under some approximation, and compute observables as above.
However, in the real world many systems have a large number of degrees of freedom that are not dynamical on the timescales of interest. For example, a spin glass has impurities that don't move; in a glass, the molecules are not free to move arbitrarily, but maintain their neighbours (roughly speaking); in an economy, we may be interested in dynamics over a background of fixed technology and resources. Such systems are called disordered, where the `disorder' is the set of fixed DOF. Call this $J$.
To treat disordered systems, ideally we would like to compute $\mathcal{O}$ at a fixed value of disorder, i.e. $\mathcal{O}[J]$. This is generally impossible because the disorder creates too many parameters to keep track of. So we assign some probability distribution to the disorder $J$, and then we would be happy to compute
$\overline{\mathcal{O}[J]}$
where $\overline{ \cdot }$ denotes the average over the distribution of disorder. Now the problem is that to average $\log Z[J]$, we seem to need to compute $Z[J]$, which is the impossible task. The replica method is based on the trivial identity
$\log z = \left.\frac{\partial z^n}{\partial n}\right|_{n=0}$
which allows us to replace computation of $\log Z[J]$ by $Z[J]^n$. If we can compute this for general $n$, then we can compute $\overline{\mathcal{O}[J]}$ as desired.
Now if $n$ is an integer, then $Z[J]^n$ is the partition function for a system of $n$ replicas of the original system, where the disorder is the same in all the replicas. Generally one considers models where it is easy to average $Z[J]^n$ over the disorder, so we are back to a pure system. But, in this new pure system, the replicas are coupled. These new interactions have the following physical interpretation: in our original system, we are exploring dynamical configurations over a landscape of fixed disorder. In general the disorder restricts the configurations that are possible (or that have a non-negligible probability). (This is what is called `frustration'). Thus, for a given disorder, the different dynamical configurations are correlated. In the replica method, once we have averaged over the disorder, what remains is this inter-replica coupling.
What happens typically is that if frustration is weak, the inter-replica coupling is not important, $\overline{Z^n} \approx \overline{Z}^n$, and the replicas drop out of the problem. In this regime it is equivalent to average $Z$ over the disorder, so physically the disorder may as well have been an extra set of dynamical variables.
Instead if frustration is strong, the inter-replica coupling is important, and the disorder-averaged partition function does not factorize over replicas. In mean-field models, the latter can be reduced to depend only on an $n\times n$ matrix $Q$, called the overlap matrix.
Originally there was a permutation symmetry over the replicas. But, like in any statistical mechanical problem, this symmetry can be broken. Giorgio Parisi understood what is the correct pattern of replica symmetry breaking that applies to the overlap matrix.
It turns out that there is a family of matrices, called hierarchical, in terms of which one can take the $n \to 0$ limit required above. It is important that the matrix size $n$ is under analytical control. Amazingly, in the limit $n \to 0$ there is a new phase of matter, the Parisi-Gardner phase${}^*$, in which the energy landscape is fractal: basins have sub-basins have sub-sub-basins, and so on, to eternity. This is a new phase of matter in which disorder has a very strong effect on the physics. For examples, response functions are generally jerky, leading to avalanches, etc.
Giorgio and collaborators worked out many of the physical consequences of this phase, and applied these constructions to many problems, among them spin glasses, structural glasses, neural networks, random manifolds, and constraint satisfaction problems. Others have applied it to a variety of problems (economies, ecosystems, languages, quantum gravity).
${}^*$ Elisabeth Gardner showed in 1985 that one can have a phase transition from a simple glass to the more exotic Parisi-Gardner phase.
