# Monte Carlo and nonlinearities

Jérémi Dauchet1 Jean-Jacques Bezian2 Stéphane Blanco3 Cyril Caliot4 Julien
Charon1 Christophe Coustet5 Mouna El Hafi2 Vincent Eymet5 Olivier Farges6
Vincent Forest5 Richard Fournier3 Mathieu Galtier7 Jacques Gautrais8 Anaïs
Khuong8 Lionel Pelissier9 Benjamin Piaud5 Maxime Roger7 Guillaume Terrée2
Sebastian Weitz1

###### Abstract

The Monte Carlo method is widely used to numerically predict systems
behaviour. However, its powerful incremental design assumes a strong premise
which has severely limited application so far: the estimation process must
combine linearly over dimensions. Here we show that this premise can be
alleviated by projecting nonlinearities on a polynomial basis and increasing
the configuration-space dimension. Considering phytoplankton growth in light-
limited environments, radiative transfer in planetary atmospheres,
electromagnetic scattering by particles and concentrated-solar-power-plant
productions, we prove the real world usability of this advance on four test-
cases that were so far regarded as impracticable by Monte Carlo approaches. We
also illustrate an outstanding feature of our method when applied to sharp
problems with interacting particles: handling rare events is now
straightforward. Overall, our extension preserves the features that made the
method popular: addressing nonlinearities does not compromise on model-
refinement nor system complexity and convergence rates remain independent of
dimension.

  1. 1.

Université Clermont Auvergne, Sigma-Clermont, Institut Pascal, BP 10448,
F-63000 Clermont-Ferrand, France.

  2. 2.

Université Fédérale de Toulouse Midi-Pyrénées, Mines Albi, UMR CNRS 5302,
Centre RAPSODEE, Campus Jarlard, F-81013 Albi CT Cedex.

  3. 3.

LAPLACE, Université de Toulouse, CNRS, INPT, UPS, France.

  4. 4.

Processes, Materials and Solar Energy Laboratory, PROMES, CNRS, 7 rue du Four
Solaire, 66120 Font-Romeu-Odeillo-Via, France.

  5. 5.

Méso-Star SAS, www.meso-star.com.

  6. 6.

Total Energies Nouvelles, R&D - Concentrated Solar Technologies, Tour
Michelet, 24 cours Michelet, 92069 Paris La Défense Cedex, France.

  7. 7.

Univ Lyon, CNRS, INSA-Lyon, Université Claude Bernard Lyon 1, CETHIL UMR5008,
F-69621, Villeurbanne, France.

  8. 8.

Research Center on Animal Cognition (CRCA), Center for Integrative Biology
(CBI), Toulouse University, CNRS, UPS, France.

  9. 9.

UMR MA 122 EFTS (Education, Formation, Travail, Savoirs), Université Toulouse
2, Toulouse.

The standard Monte Carlo (MC) method is a method to predict a physical
observable by numerically estimating a statistical expectation over a multi-
dimensional configuration space[1]. The reason why this method is so popular
in all fields of scientific research is its intuitive nature. In the most
current practice, simulation tools are designed in direct relation with the
pictures available in each discipline and later refinements are gradual and
straightforward. Model refinements merely extend sampling to new suitable
dimensions. The method is nonetheless mathematically rigorous: specialists
specify observables that are implicitly translated into integral quantities
that are estimated using random sampling in each direction of the
configuration space. This choice of a statistical viewpoint is highly powerful
because the algorithm can be designed directly from the description of the
system, whether it is deterministic or not, with no reworking nor
approximation.

Let us illustrate how MC is used in engineering with a typical example: the
optimal design of a concentrated solar plant[2] (see Fig.1-a). The power
collected by the central-receiver results from all the rays of sunlight that
reach it after reflection by heliostats, so it depends on the complex geometry
of the heliostats. Moreover, the heliostats being adaptively oriented to
follow the sun position, they can mask one another at some times of the day.
To estimate by MC the received power at one time of interest, i.e. for a given
geometry of the heliostats: pick an optical path among the ones that join the
sun to the central-receiver via an heliostat, check whether this path is
shadowed or blocked by another heliostat, and retain a Monte Carlo weight
equal to $0$ or $1$ depending on the transmission success. Let $\bf X$ be the
random variable denoting the transmission success. The collected fraction of
the available sun power is then the expectation ${\cal E}_{\bf X}(\bf X)$ of
$\bf X$ and can be evaluated with no bias as the average of such weights over
a large number of sampled paths.

This approach robustly complies with expanded descriptions of the physical
observable to be addressed. For instance, the fraction of the available sun
power collected on average over the entire lifetime of the plant (operating
typically 30 years) can be predicted as the expectation over time of ${\cal
E}_{\bf X}(\bf X)$ varying in time. Denoting ${\cal E}_{\bf X|Y}(\bf X|Y)$ the
collected fraction at random time $\bf Y$ within the 30 years, the time-
averaged fraction is given by ${\cal E}_{\bf Y}({\cal E}_{\bf X|Y}({\bf
X|Y}))={\cal E}_{{\bf Y},{\bf X|Y}}({\bf X|Y})$. The basic algorithm above can
then be encapsulated within time sampling: first pick a date for $\bf Y$, then
pick a path at that date for $\bf X|Y$. In the end, estimate ${\cal E}_{{\bf
Y},{\bf X|Y}}({\bf X|Y})$ by computing the average of the transmission success
over all combined pairs (date, path). Meanwhile, one can as well account for
sun power fluctuations by estimating the atmospheric transmission at each
picked date. The choice of the statistical viewpoint thus allows to
incorporate into one single statistical question as many ingredients as it
takes: the geometrical complexity of the heliostats[3], the daily course of
the sun, the seasonal scale as well as the hour scale weather fluctuations[6].
Remarkably, the latter question is nearly as simple to address as the
estimation of the power collected at one single date: the algorithmic design
can map the full conceptual description, yet computational costs are poorly
affected. Contrastingly, deterministic approaches would translate into
unpractical computation times or require simplified and approximate
descriptions, so MC has become the only practical solution in many such
engineering contexts. Becoming such a standard practice, MC has prompted
numerous theoretical developments[6, 7, 5, 8].

Figure 1 — Complex systems with non linear outputs: four real-world examples.

#### Fig. 1 Complex systems with nonlinear outputs: four real world examples.

a, Solar-driven high-temperature thermal reduction of zinc oxide, as the first
part of a two step water splitting cycle. Photons emitted from the sun are
reflected on heliostats and concentrated at the entrance of the chemical
reactor in which $ZnO$ dissociation is carried out. Depending on their
transmission success $\bf X$, the solar power ${\cal E}_{\bf X|\bf
Y}\left({\bf X|\bf Y}\right)$ absorbed by the receiver at a random instant
${\bf Y}$ of lifetime determines the nonlinear chemical conversion rate of the
reaction $ZnO\rightarrow Zn+\frac{1}{2}O$. Here we address the estimation of
the annual solar-plant’s conversion rate $C(\bf p)$ at different earth
locations $\bf p$, by averaging the instantaneous conversion rates over the
statistics of sun position and incident Direct Normal Irradiance (DNI), which
fluctuates with time and weather (see also EDF1 and SI1).

b, Wave scattering by a complex-shaped and optically-soft scatterer
(cyanobacterium Arthrospira). An incident plane wave with propagation
direction ${\bf e}_{i}$ is scattered by the helical cyanobacterium. The
bacterium has low relative refractive index and is much larger than wavelength
(optically-soft particle). The complex scattering-amplitude ${\cal E}_{\bf
X|\bf Y}\left({\bf X|\bf Y}\right)$ in the forward directions is the sum of
secondary waves contributions ${\bf X|\bf Y}$ (interference) originating from
the scatterer projected surface. This surface depends on the scatterer
orientation ${\bf Y}$. Here we address the estimation of $S({\bf e}_{s})$ the
single-scattering differential cross-section in direction ${\bf e}_{s}$ for a
suspension of particles, assuming independent scattering, by averaging the
squared modulus of ${\cal E}_{\bf X|\bf Y}\left({\bf X|\bf Y}\right)$ over the
statistics of orientations ${\bf Y}$ (see also EDF2 and SI2).

c, Phytoplankton growth in light-limited environments. Phytoplankton is put to
grow in a continuous stirred tank photobioreactor internally illuminated by
optical fibres ${\cal F}$ immersed in the culture. The local rate of photon
absorption ${\cal E}_{\bf X|\bf Y}\left({\bf X|\bf Y}\right)$ at location $\bf
Y$ is the average of the contributions ${\bf X}|{\bf Y}$ of every optical path
from fibres to $\bf Y$ through the scattering and absorbing suspension. ${\cal
E}_{\bf X|\bf Y}\left({\bf X|\bf Y}\right)$ determines the nonlinear
photosynthetic growth-rate at location ${\bf Y}$. Here, we address the Monte
Carlo estimation of $R(C)$ the global growth-rate in the whole culture volume
as a function of biomass concentration $C$, by averaging the local rate over
locations in the volume (see also EDF3 and SI3).

d, Atmospheric radiative transfer : top-of-atmosphere (TOA) specific intensity
(from earth toward the outer space). Photons emitted by the atmosphere at
infrared frequencies are due to random emission-transitions ${\bf Y}$, from an
upper to a lower energy-state, of mainly $\mathrm{CO_{2}}$ and
$\mathrm{H_{2}O}$ molecules of concentration $C_{\bf Y}$ at altitude $H_{\bf
Y}$. The corresponding source ${\cal S}_{\bf Y}$ depends on the thermodynamic
state of the atmosphere, mainly temperature $T_{\bf Y}$ (defining the energy-
state population) and pressure $P_{\bf Y}$ (defining most of the line-width,
i.e. the uncertainty of the emission-frequency $\nu_{\bf Y}$). This source is
then exponentially extinct by atmospheric absorption, i.e. by all random
absorption-transitions ${\bf X}|{\bf Y}$, from a lower to an upper energy-
state, occurring at altitude $H_{\bf X|Y}$ between $H_{\bf Y}$ and the top of
atmosphere (see also EDF4 and SI4).

Still, MC cannot handle _any_ question to date. As a matter of fact, it has
been identified early on that “the extension of Monte Carlo methods to
nonlinear processes may be impossible”[9] and it is a prevalent opinion
nowadays that “Monte Carlo methods are not generally effective for nonlinear
problems mainly because expectations are linear in character”[10] so that “a
nonlinear problem must usually be linearized in order to use Monte Carlo
technique”[11]. We are aware of only one attempt so far to bypass this
failing: the proposition by Vajargah[12] or Dimov[1] of using branching
processes to solve Fredholm-type integral equations with polynomial
nonlinearity.

Unluckily, most real world problems are nonlinear. And indeed, if the question
were now to evaluate the final return on investment of the solar plant, namely
how much electrical power it would deliver over its lifetime, standard MC
would fail because the instantaneous conversion efficiency from collected
solar-power to electrical power is not linear. Let us consider, as a toy
example, a basic nonlinear case where the electrical power would be
proportional to the square of the instantaneous collected solar power ${\cal
E}_{\bf X|Y}(\bf X|Y)$ at date $\bf Y$. In Monte-Carlo terms, the question
would be then to estimate ${\cal E}_{\bf Y}({\cal E}_{\bf X|Y}({\bf
X|Y})^{2})$ over the plant lifetime. In this case, the optical and temporal
expectations cannot be combined any more because it would be wrong to first
estimate, as above, the total solar power collected over lifetime, and then
apply the conversion efficiency in the end (basically, ${\cal E}_{\bf Y}({\cal
E}_{\bf X|Y}({\bf X|Y})^{2})\neq{\cal E}_{\bf Y}({\cal E}_{\bf X|Y}({\bf
X|Y}))^{2}$ the same way as $a^{2}+b^{2}\neq(a+b)^{2}$). Instead, one would
have to sample dates (say $M$ dates, millions over 30 years), estimate the
solar power collected at each date by averaging over numerous optical paths
transmission successes (say $N$ paths, millions for each date), apply the
nonlinear conversion to the result at that date, and then average over all
dates. Doing so, MC would now require $M\times N$ samples, and, even worse,
further levels of complexity (each adding a nonlinearity to the problem) would
similarly multiply the computation time. Moreover, the result would be biased
due to the finite sampling sizes of the innermost dimensions. In short, MC
distinctive features are dismissed and exact life-time integration looks
impossible.

Bearing in mind our earlier theoretical works about MC integral
formulations[2], we found a way to bypass this obstacle for a large class of
nonlinear problems, based on the very statistical nature of MC. In the case of
our toy example, we use the fact that:

${\cal E}_{\bf Y}({\cal E}_{\bf X|Y}({\bf X|Y})^{2})={\cal E}_{{\bf
Y},{\bf(X_{1},\,X_{2})|Y}}({\bf X_{1}\,X_{2}|Y})$ (1)

where ${\bf X}_{1}$ and ${\bf X}_{2}$ are two independent variables
identically distributed as ${\bf X}$ (see Methods). Translated in sampling
algorithm, the solution is then to sample optical paths in pairs
${\bf(X_{1},X_{2})|Y}$ (instead of millions) at each sampled date and now
retaining the pair product ${\bf X_{1}X_{2}|Y}$ of their transmission
successes. The optical and temporal statistics can then actually be sampled
together and yield the unbiased result with no combinatorial explosion. This
reformulation can be generalised to any nonlinearity of polynomial shape.
First, monomials of any degree can indeed be estimated using the same sampling
property as the one used above for $n=2$:

${\cal E}_{\bf Y}({\cal E}_{\bf X|Y}({\bf X|Y})^{n})={\cal E}_{{\bf
Y},{\bf(X_{1},\,X_{2},...,\,X_{n})|Y}}({\bf X_{1}\,X_{2}...X_{n}|Y})$ (2)

where the ${\bf X}_{i}$ are $n$ independent random variables identically
distributed as ${\bf X}$. For any monomial of degree $n$, the expectation can
then be computed by sampling series of $n$ independent realisations of $\bf
X|Y$, and averaging the series products. The linear case, solved by standard
MC, corresponds to $n=1$. Secondly, since polynomials are just linear
combinations of monomials, the expectation for any polynomial function of
${\cal E}_{\bf X|Y}({\bf X|Y})$ of degree $n$ can be translated into a Monte
Carlo algorithm, sampling first a degree in the polynomial, and then sampling
as many independent realisations of ${\bf X|Y}$ as this random degree (see
Methods). For a polynomial function of degree $n$, the corresponding Non-
Linear Monte Carlo (NLMC) algorithm is then:

  * •

pick a sample $y$ of ${\bf Y}$,

  * •

pick a monomial degree value $d\leq n$,

  * •

draw $d$ independent samples of ${\bf X|Y=y}$ and retain their product,

repeat this sampling procedure and compute the estimate as the average of
retained products.

Moreover, if polynomial forms of any dimension are now solvable with no
approximation, so is the projection of any nonlinear function on a polynomial
basis of any dimension, even of infinite dimension if needed (full details
using the Taylor expansion are given in Methods). As a result, any hierarchy
of nested statistical processes that combine nonlinearly can now be, in
theory, exactly addressed within the Monte Carlo framework. The deep rationale
of our reformulation of nonlinear expectations can be grounded on the
fundamentals of functional statistics: nonlinear processes are formally
equivalent to infinite-dimension processes and MC algorithms are well known to
smartly address, by design, expectations over domains of infinite dimension.

To our knowledge, this analysis has never been made. It has yet major
practical consequences for real world problems, provided the polynomial
sampling, which is the price to be paid for tackling nonlinearities exactly,
remains tractable. For instance, let’s go back to our illustrative solar-
power-plant example, and use now the actual expression for the conversion rate
and its Taylor expansion: for each date, once a sun-position and climate
conditions have been fixed, one would have to pick a random number of
independent optical paths (instead of one optical path in the linear case),
keep the product of their transmission success, and make the average of many
such products in the end. Doing so, it becomes possible to integrate hourly
solar input fluctuations over 30 years in the full geometry of a kilometre-
wide heliostats field in order to optimise the nonlinear solar-to-electric
conversion over the plant lifetime (Fig. 1A). The same line of thought can be
used to predict wave scattering by a tiny scale complex-shaped scatterer[7]
such as an helicoidal cyanobacterium (Fig. 1B). The biomass production example
(Fig. 1C), where incoming light favours the photosynthetic growth that in turn
blocks the incoming light, illustrates how our method also handles nonlinear
feedbacks[14]. Finally, with the estimation of Earth radiative cooling (Fig.
1D), we reproduce quite classical results, yet with a pure statistical
approach: by sampling directly the state transitions of greenhouse gases, we
avoid costly deterministic computations that the standard linear Monte Carlo
approach requires in order to by-pass the nonlinearity of Beer Extinction
Law[15]. In each of the four cases, it appeared that the additional
computations were well affordable using only ordinary computing power (the
complete physical descriptions of the four problems, the nonlinearities
involved and their translation in NLMC can be found in their respective
Extended Data Figures and Supplemental Information, Solar Plant: EDF1, SI1 ;
Complex-shaped Scatterer: EDF2, SI2 ; Biomass production: EDF3, SI3 ; Earth
radiative cooling: EDF4, SI4).

From a theoretical point of view, in these four cases, the model is directly
enunciated in statistical terms, defining two random variables ${\bf X}$ and
${\bf Y}$ from the start. More broadly, standard MC practice can start as well
from a deterministic description (see Methods), most commonly from a linear
partial differential equation (PDE). The formal equivalence between the
solution of a linear PDE and the expectation of a random variable has been
established for long[16]. And indeed PDE-to-MC translations are essential to
nanoscale mechanics (Quantum Monte Carlo[17]) or nuclear sciences. NLMC allows
such translations for nonlinear PDE.

As a touchstone, we address a prominent example of nonlinear PDE in
statistical physics, the Boltzmann equation which governs the spatiotemporal
density of interacting particles in a dilute gas (full details in SI5). The
corresponding physics is easy to visualise: a particle just follows its
ballistic flight until it collides with another particle. The collisions are
considered as instantaneous and only modify the two particles velocities. The
equation for the evolution of particle density in phase-space (position,
velocity) is nonlinear because the collision rate depends on the density
itself. In order to project this nonlinearity on the proper polynomial basis
of infinite dimension, this PDE is first translated into its Fredholm integral
counterpart (a step reminiscent of the aforementioned Dimov’s proposition[1]).
This Fredholm integral expresses the density in the phase space at some
location for some velocity at some time, as if putting a probe in space-time.
It is estimated by Monte Carlo, tracking the dynamics backward in time up to
the initial condition (or boundary conditions). Importantly, such a probe
estimation does not require the exhaustive resolution of the whole field at
previous times: like in standard backward MC for linear transport, the
information about previous states of the field is reconstructed along each
path only where and when it is required[20]. The contrast with linear MC is
that nonlinearity due to collisions translate into branching paths.

This extension deals very efficiently with extremely rare events because it
preserves an essential feature of MC: avoiding time / space / velocity
discretisation[19, 20, 21], very low densities can be estimated with no bias,
and the only source of uncertainty is the finite number of sampled events
(i.e. the confidence interval around the estimated density). As a test bed, we
consider a case for which analytical solutions have been published: the early
Krook’s analysis of the distribution of speeds in far out-of-equilibrium
conditions[22, 23]. Krook’s analysis was outstanding because it provided an
analytical solution to a problem which looked impossible to solve numerically:
events with the greatest consequences, namely particles with the highest
energies (i.e. high speed particles, of tremendous importance in nuclear
chemistry) lie far in the tail of speed distribution and have very low
probability of occurrence (rare events). Using our NLMC design, the fractions
of particles which have a kinetic energy higher than $10^{6}$ times the
average value, and which correspond to a fraction as low as $10^{-11}$ of the
total, can now be quantified as accurately as desired, and fit perfectly the
analytical solution (fig. 2a).

Being validated on Krook’s case, this extension opens the way to solve systems
for which no analytical solutions are available. As an example, we consider
now a fully spatialised system in which the particles are confined by an
outside harmonic potential, leading to a so-called _breathing mode_ of the gas
density. Such a solution to the Boltzmann equation had been early identified
by Boltzmann himself[24], but has been recently revisited and generalised in
the context of shortcut to adiabaticity techniques for classical gases[25].
Exact solutions are available only under the constraint that gas is at local
equilibrium, in which case density displays a permanent oscillation. Here
again, these analytical solutions are exactly recovered. Yet, NLMC allows to
go beyond this constraint and to explore the gas behaviour when the local
equilibrium constraint is alleviated: starting from a state far from local
equilibrium, it is now possible to estimate how fast the velocities
redistribution induced by collisions actually dampens the oscillation (fig.
2b).

Figure 2 — Non linear Monte Carlo for gas kinetics.

#### Fig. 2 Non-Linear Monte Carlo for gas kinetics.

a, Relaxation of speed distribution to equilibrium (a1: tail distribution of
particles speeds (fraction of particles faster than $\bf S$), a2: probability
density function of velocity, a3: evolution of the fraction of particles
faster than 6 speed units). In an homogenous gas, collisions between particles
redistribute velocities so that the speed distribution tends to equilibrium
(Maxwellian distribution). Starting from a distribution far from equilibrium
(black curves), we compute the relaxation to the Maxwellian distribution (blue
curves) by estimating the tail distribution at different times (e.g. red
curves correspond to the system state at 1 unit time). The continuous lines
correspond to the analytical solutions and each point corresponds to an
independent NLMC estimation. The fraction of particles faster than 6 speed
units (a3) illustrates how NLMC well accounts for the 1000-fold increase of
the rarest high-speed particles, with no space nor time discretisation.
Remarkably, rare events are estimated with the same relative uncertainty as
frequent events ($10^{4}$ samples for each estimate, confidence intervals of
all estimates are contained within the point thickness).

b, Dampening of breathing mode. A dilute gas confined by an outside harmonic
potential display a _breathing mode_. We estimate the density at probe
position $(1.75,0,0)$ and velocity $(0.35,0,0)$ (in adimensional units, see
SI5) at different times (adim. unit). Starting from a distribution complying
with local equilibrium, the density displays a perpetual oscillation at twice
the trap frequency (blue curve, analytical solution), independent of the
collision rate (or, equivalently, of the elastic cross section). Starting from
a distribution far from local equilibrium (the same initial distribution as in
a), the density still pulsates in absence of collisions (null cross section,
black curve: analytical solution, black points: probe estimates). Introducing
collisions (raising cross section), the velocities redistribution induced by
collisions dampens the oscillation (red points: probe estimates, no analytical
solution). The inset illustrates that probes can be bunched to zoom in on a
period of special interest (e.g. estimating the peak values at each cycle).
Each point is estimated independently using $10^{7}$ samples, confidence
intervals of all estimates are contained within the point thickness in the
main figure, confidence intervals are represented by the gray background area
in the inset).

## Conclusions

From now on, the Monte Carlo Method is no more restricted to linear problems.
The five examples exposed above have been worked out by teams gathering
specialists of Monte Carlo method and specialists of the physical problem
under consideration. By their full exposition, we offer the readers all the
details to implement their own applications. As a guideline, the first step is
to formulate the physical observable under its expectation form, including the
non linearities and integrating all levels of complexity. The second step is
to reformulate this expectation as a formulation compliant with standard Monte
Carlo Method, according to the kind of non linearities. For polynomial non
linearities, use i.i.d. series products. For other differentiable forms, use
Taylor expansion around an upper-bound of the innermost random variable in
order to regain a polynomial form. Using this MC-compliant formulation, every
advanced MC technique can then be applied: parallel implementation, complex
geometry, null-collisions, zero-variance, control variate, importance
sampling, sensitivities analysis… As illustrated by the variety of our seminal
examples, this guideline covers a large set of non linear academic and real
world problems.

## References

  * [1] Dimov, I. T. & McKee, S. Monte Carlo Methods for Applied Scientists. (World Scientific Publishing, 2008).
  * [2] Delatorre, J. et al. Monte Carlo advances and concentrated solar applications. Sol. Energy 103, 653-681 (2014).
  * [3] Siala, F. M. F. & Elayeb, M. E. Mathematical formulation of a graphical method for a no-blocking heliostat field layout. Renew. energy 23, 77-92 (2001).
  * [4] Farges, O. et al. Life-time integration using Monte Carlo Methods when optimizing the design of concentrated solar power plants. Sol. Energy 113, 57-62 (2015).
  * [5] Assaraf, R. & Caffarel, M. Zero-variance principle for Monte Carlo algorithms. Phys. Rev. Lett. 83, 4682 (1999).
  * [6] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087-1092 (1953).
  * [7] Hammersley, J. M. & Handscomb, D. C. Monte Carlo Methods. (Springer Netherlands, 1964).
  * [8] Roger, M., Blanco, S., El Hafi, M. & Fournier, R. Monte Carlo Estimates of Domain-Deformation Sensitivities. Phys. Rev. Lett. 95, 180601 (2005).
  * [9] Curtiss, J. H. Monte Carlo Methods for the Iteration of Linear Operators. J. Math. Phys. 32, 209-232 (1953).
  * [10] Kalos, M.H. & Whitlock, P.A. Monte Carlo Methods, second ed. (Wiley-VCH, Weinheim 2008).
  * [11] Chatterjee, K., Roadcap, J. R. and Singh, S. A new Green’s function Monte Carlo algorithm for the solution of the two-dimensional nonlinear Poisson-Boltzmann equation: Application to the modeling of the communication breakdown problem in space vehicles during re-entry. J. Comput. Phys. 276, 479-485 (2014).
  * [12] Vajargah, B. F. & Moradi, M. Monte Carlo algorithms for solving Fredholm integral equations and Fredholm differential integral equations. Appl. Math. Sci. 1, 463-470 (2007).
  * [13] Charon, J. et al. Monte Carlo implementation of Schiff’s approximation for estimating radiative properties of homogeneous, simple-shaped and optically soft particles: Application to photosynthetic micro-organisms. J. Quant. Spectrosc. Radiat. Transf. 172, 3-23 (2016).
  * [14] Cornet, J. F. Calculation of optimal design and ideal productivities of volumetrically lightened photobioreactors using the constructal approach. Chem. Eng. Sci. 65, 985-998 (2010).
  * [15] Galtier, M. et al. Radiative transfer and spectroscopic databases: A line-sampling Monte Carlo approach. J. Quant. Spectrosc. Radiat. Transf. 172, 83-97 (2016).
  * [16] Kac, M. On some connections between probability theory and differential and integral equations. Proc. Second Berkeley Symp. Math. Statistics Probab. 189 (1951).
  * [17] Corney, J. F. & Drummond, P. D. Gaussian quantum Monte Carlo methods for fermions and bosons. Phys. Rev. Lett. 93, 2-5 (2004).
  * [18] Galtier, M. et al. Integral formulation of null-collision Monte Carlo algorithms. J. Quant. Spectrosc. Radiat. Transf. 125, 57-68 (2013).
  * [19] Wagner, W. Stochastic particle methods and approximation of the Boltzmann equation. Math. Comput. Simul. 38, 211-216 (1995).
  * [20] Rjasanow, S. A Stochastic Weighted Particle Method for the Boltzmann Equation. J. Comput. Phys. 124, 243-253 (1996).
  * [21] Rjasanow, S. & Wagner, W. Simulation of rare events by the stochastic weighted particle method for the Boltzmann equation. Math. Comput. Model. 33, 907-926 (2001).
  * [22] Krook, M. & Wu, T. T. Formation of Maxwellian Tails. Phys. Rev. Lett. 36, 1107-1109 (1976).
  * [23] Krook, M. & Wu, T. T. Exact solutions of the Boltzmann equation. Phys. Fluids 20, 1589-1595 (1977).
  * [24] L. Boltzmann, in Wissenschaftliche Abhandlungen, edited by F. Hasenorl (J.A. Barth, Leipzig, 1909), Vol. II, p. 83.
  * [25] Guéry-Odelin, D., Muga, J. G., Ruiz-Montero, M. J. & Trizac, E. Nonequilibrium Solutions of the Boltzmann Equation under the Action of an External Force. Phys. Rev. Lett. 112, 180602 (2014).
  * [26] Weitz, S. et al. Monte Carlo efficiency improvement by multiple sampling of conditioned integration variables. J. Comput. Phys. 326, 30-34 (2016).

Competing Interests

The authors declare that they have no competing financial interests.

Correspondence

Correspondence and requests for materials should be addressed to Jérémi
Dauchet (email: jeremi.dauchet@sigma-clermont.fr).

Authors contributions

All authors contributed extensively to the theoretical developments presented
in this paper. Each author contributed to the practical applications according
to its scientific expertise: J.D., M.E.H., V.E., R.F. and M.G. in Atmospheric
sciences, J.D., S.B., C.Ca., M.E.H., V.E., O.F., R.F., M.G. and M.R. in
Radiative Transfert, J.D., S.B., R.F., J.G., A.K. and S.W. in Complex Systems
in Biology, J.D., J.J.B., S.B., C.Ca., M.E.H., V.E., O.F. and R.F. in Solar
Energy, J.D., S.B., J.C., M.E.H. and R.F. in Electromagnetism and quantum
mechanics, J.D., S.B., M.E.H., R.F., J.G., B.P. and G.T. in Statistical
Physics. C.Co., V.E., V.F. and B.P. (www.meso-star.com) performed the
numerical implementations. J.D., S.B., R.F. and J.G. wrote the paper.

## Methods

#### Basics of Monte Carlo Methods.

Let us estimate $E=1+4$ by repeatedly tossing a (fair) coin. The tossing
process is described by a random variable $R\in\\{0,1\\}$ that takes the value
$0$ for head (probability $P_{R}(0)=\frac{1}{2}$) and $1$ for tail
(probability $P_{R}(1)=\frac{1}{2}$).

Now, to estimate any process (e.g. a process output: $E=1+4$), we can affect
arbitrary weights $w(R)$ to values $\\{0,1\\}$ in order to write $E$ as an
expectation of the weighted process, following:

$E=1+4=P_{R}(0)w(0)+P_{R}(1)w(1)={\cal E}_{R}\left(w(R)\right)$ (3)

with $w(0)=\frac{1}{P_{R}(0)}=2$ and $w(1)=\frac{4}{P_{R}(1)}=8$ and where
${\cal E}_{R}$ denotes expectation with respect to $R$. Using the results
$r_{1}...r_{N}$ of $N$ successive tosses (independent realisations of $R$), we
can then estimate $E={\cal E}_{R}\left(w(R)\right)$ from the weighted average
of the tosses’ results $\frac{1}{N}\sum_{i=1}^{N}w(r_{i})$ since $E=5$ is
indeed the average of Monte Carlo weights that take the values $2$ and $8$
with equal probabilities.

Such an approach is at the basement of Monte Carlo techniques: define the
weights according to the problem to be solved, sample the process repeatedly,
and take the average. Depending on the physical description of the value to be
estimated, this approach still holds for an infinite number of terms and can
also be extended to integrals formulation using continuous random variables:

${\cal E}_{\bf Y}\left(w({\bf Y})\right)=\int_{{\cal D}_{\bf Y}}d{\bf
y}\,p_{\bf Y}({\bf y})\,w({\bf y})$ (4)

which can be estimated by $\frac{1}{N}\sum_{i=1}^{N}w({\bf y}_{i})$, where the
$y_{i}$ are $N$ realisations of the random variable $Y$ with probability
density function $p_{Y}$ and domain of definition ${\cal D}_{\bf Y}$.

#### Basics of Non-Linear Monte Carlo Methods.

In order to estimate

$E={\cal E}_{\bf Y}({\cal E}_{\bf X|Y}({\bf X|Y})^{2})=\int_{{\cal D}_{\bf
Y}}d{\bf y}\,p_{\bf Y}(y)\left(\int_{{\cal D}_{\bf X|Y}}dx\,p_{\bf
X|Y}(x|y)\;x\right)^{2}$ (5)

we introduce two independent variables ${\bf X}_{1}$ and ${\bf X}_{2}$
identically distributed as ${\bf X}$ (still conditioned by the same ${\bf
Y}$):

$\begin{array}[]{lcl}E&=&{\cal E}_{\bf Y}\left({\cal E}_{\bf X_{1}|Y}({\bf
X_{1}|Y})\,{\cal E}_{\bf X_{2}|Y}({\bf X_{2}|Y})\right)\\\
&=&\displaystyle\int_{{\cal D}_{\bf Y}}dy\,p_{\bf Y}(y)\left(\int_{{\cal
D}_{\bf X|Y}}dx_{1}\,p_{\bf X|Y}(x_{1}|y)\;x_{1}\right)\left(\int_{{\cal
D}_{\bf X|Y}}dx_{2}\,p_{\bf X|Y}(x_{2}|y)\;x_{2}\right)\end{array}$ (6)

Since ${\bf X}_{1}$ and ${\bf X}_{2}$ are independent, and conditionally
independent given $\bf Y$:

$\begin{array}[]{lcl}E&=&\displaystyle\int_{{\cal D}_{\bf Y}}dy\,p_{\bf
Y}(y)\left(\iint_{{\cal D}^{2}_{\bf X|Y}}dx_{1}\,p_{\bf
X|Y}(x_{1}|y)\;dx_{2}\,p_{\bf X|Y}(x_{2}|y)\;x_{1}x_{2}\right)\\\ &=&{\cal
E}_{\bf Y}\left({\cal E}_{\bf(X_{1},\,X_{2})|Y}({\bf
X_{1}X_{2}|Y})\right)\end{array}$ (7)

Hence

$\begin{array}[]{lcl}E&=&\displaystyle\iiint_{{\cal D}_{\bf Y}\times{\cal
D}^{2}_{\bf X|Y}}dy\,p_{\bf Y}(y)\;dx_{1}\,p_{\bf
X|Y}(x_{1}|y)\;dx_{2}\,p_{\bf X|Y}(x_{2}|y)\;x_{1}x_{2}\\\ &=&{\cal E}_{{\bf
Y},{\bf(X_{1},\,X_{2})|Y}}\left({\bf X_{1}X_{2}|Y}\right)\end{array}$ (8)

The same demonstration can be drawn to establish that:

${\cal E}_{\bf Y}({\cal E}_{\bf X|Y}({\bf X|Y})^{n})={\cal E}_{{\bf
Y},{\bf(X_{1},\,X_{2},...,\,X_{n})|Y}}({\bf X_{1}~{}X_{2}...X_{n}|Y})$ (9)

Let us now assume that the weights associated to the random variable ${\bf Y}$
are described by a nonlinear function ${\bf f}({\bf Z_{\bf Y}})$ of the
conditional expectation ${\bf Z_{\bf Y}}={\cal E}_{\bf X|Y}\left({\bf X|\bf
Y}\right)$. The problem then becomes to compute:

$E={\cal E_{\bf Y}}{\bf\left(f(Z_{\bf Y})\right)}={\cal E}_{\bf Y}\left({\bf
f}\left({\cal E}_{\bf X|Y}\left({\bf X|\bf Y}\right)\right)\right)$ (10)

Such nonlinearity can be treated with no approximation using a projection on
infinite basis. In all the examples presented in this article, we have used
Taylor polynomials basis, which means that ${\bf f}(x)$ is expanded around
$x_{0}$

${\bf f}(x)=\sum_{n=0}^{+\infty}\frac{\partial^{n}{\bf
f}(x_{0})}{n!}\left(x-x_{0}\right)^{n}$ (11)

We note that both $x_{0}$ and ${\bf f}$ can be conditioned by $\bf Y$. Now,
following the same line as explained above for the Basics of Monte Carlo
Methods, we regard the sum in the expansion of ${\bf f}$ as an expectation,
writing:

${\bf f}(x)={\cal E}_{N}\left(\frac{\partial^{N}{\bf
f}(x_{0})}{P_{N}(N)N!}\left(x-x_{0}\right)^{N}\right)$ (12)

where the random variable $N$ (of probability law $P_{N}$) is the degree of
one monomial in the Taylor polynomial. This step only requires to define an
infinite set of probabilities (instead of two in Eq. 3), with
$\sum_{n=0}^{+\infty}P_{N}(n)=1$. For the choice of $x_{0}$, we have used an
upper bound of $\bf X|Y$, so that $P_{N}$ is implicitly derived using
Bernoulli random variables (see details in SI 1-3).

Equation 10 can then be rewritten as:

$E={\cal E}_{\bf Y}\left({\bf f}\left({\cal E}_{\bf X|Y}\left({\bf X|\bf
Y}\right)\right)\right)={\cal E}_{Y,N}\left(\frac{\partial^{N}{\bf
f}(x_{0})}{P_{N}(N)N!}\left({\cal E}_{\bf X|Y}({\bf
X|Y})-x_{0}\right)^{N}\right)$ (13)

Defining independent and identically distributed random variables ${\bf
X}_{q}$, with same distribution as ${\bf X}$, the innermost term rewrites

$E={\cal E}_{Y,N}\left(\frac{\partial^{N}{\bf
f}(x_{0})}{P_{N}(N)N!}\prod_{q=1}^{N}\left({\cal E}_{\bf X_{q}|Y}({\bf
X}_{q}|{\bf Y})-x_{0}\right)\right)$ (14)

or, equivalently:

$E={\cal E}_{Y,N}\left(\frac{\partial^{N}{\bf
f}(x_{0})}{P_{N}(N)N!}\prod_{q=1}^{N}{\cal E}_{\bf X_{q}|Y}({\bf X}_{q}|{\bf
Y}-x_{0})\right)$ (15)

Since the variables ${\bf X}_{q}|{\bf Y}$ are independent in the innermost
term, we have:

$\prod_{q=1}^{N}{\cal E}_{\bf X_{q}|Y}({\bf X}_{q}|{\bf Y}-x_{0})={\cal
E}_{\bf(X_{1}\,,X_{2},...,X_{N})|Y}\left(\prod_{q=1}^{N}\left({\bf X}_{q}|{\bf
Y}-x_{0}\right)\right)$ (16)

so that:

$E={\cal E}_{\bf Y,N}\left(\frac{\partial^{N}{\bf
f}(x_{0})}{P_{N}(N)N!}\,{\cal
E}_{\bf(X_{1}\,,X_{2},...,X_{N})|Y}\left(\prod_{q=1}^{N}\left({\bf X}_{q}|{\bf
Y}-x_{0}\right)\right)\right)$ (17)

and we finally have:

$E={\cal E}_{\bf Y,N,(X_{1}\,,X_{2},...,X_{N})|Y}\left(\frac{\partial^{N}{\bf
f}(x_{0})}{P_{N}(N)N!}\prod_{q=1}^{N}\left({\bf X}_{q}|{\bf
Y}-x_{0}\right)\right)$ (18)

which can be read as:

$E={\cal E}_{\bf Y,N,(X_{1}\,,X_{2},...,X_{N})|Y}\left(w({\bf Y},{\bf N},{\bf
X_{1}},{\bf X_{2}},...,{\bf X_{N}})\right)$ (19)

with

$w({\bf Y},{\bf N},{\bf X_{1}},{\bf X_{2}},...,{\bf
X_{N}})=\frac{\partial^{N}{\bf f}(x_{0})}{P_{N}(N)N!}\prod_{q=1}^{N}\left({\bf
X}_{q}|{\bf Y}-x_{0}\right)$ (20)

The translation into a Monte Carlo algorithm then follows:

  * •

sample a realisation ${\bf y}$ of ${\bf Y}$ (and set $x_{0}$ and ${\bf f}$
accordingly if they depend on ${\bf y}$)

  * •

sample a realisation $n$ of $N$

  * •

sample $n$ independent realisations ${\bf x}_{q=1,...,n}$ of ${\bf X}$
conditioned by $\bf y$

  * •

keep

$\hat{w}=w({\bf y},n,{\bf x}_{1},...,{\bf x}_{n})=\frac{\partial^{n}{\bf
f}(x_{0})}{P_{N}(n)n!}\prod_{q=1}^{n}\left({\bf x}_{q}-x_{0}\right)$

and estimate $E$ as the average of many realisations $\hat{w}$.

## Supplemental Information

##  1 SI1 — Solar thermochemical reduction of zinc oxide averaged over the
year

###  1.1 The problem

Solar-driven high temperature thermochemical cycles processes, commonly based
on metal oxides reduction [1, 2], are an alternative to fossil fuel-based
method for $H_{2}$ generation. Here we focus on thermal reduction zinc oxide,
as the first part of a two step water splitting cycle. Photons emitted from
the sun are reflected on heliostats and concentrated at the entrance of the
chemical reactor in which $\mathrm{ZnO}$ dissociation is carried out. Solar
power received by the reactor at a given instant determines the chemical
conversion rate of the reaction $\mathrm{ZnO\rightarrow Zn+\frac{1}{2}O}$. The
present Monte Carlo estimation of the solar-plant’s annual conversion rate $C$
associates the thermochemical knowledge of the non-linear kinetics of zinc
oxide dissociation [3, 4] to the description of radiative transfer on a
multiple-reflection solar receiver [5]. The nonlinearity lies in the
instantaneous coupling between photon transport and zinc-oxide reduction.

The random variable ${\bf X}$ is the contribution of an optical path from the
sun to the entrance of the chemical reactor. Its expectation ${\cal E}_{{\bf
X}|{\bf Y}}({\bf X}|{\bf Y})$ is the instantaneous thermal-power fraction
collected at a given moment $\bf Y$ of the year[6]. Solar power is used to
reduce the zinc oxide with a nonlinear conversion rate $\mathrm{f}\,({\cal
E}_{{\bf X}|{\bf Y}}({\bf X}|{\bf Y}))$ (see[3, 4]).

###  1.2 Non-Linear Monte Carlo formulation

The annual conversion rate $C$ is reformulated as its Taylor expansion around
$x_{0}$. First, independent and identically distributed (i.i.d) optical-path
random variables are introduced, each defining a random contribution ${\bf
X}_{i}|{\bf Y}$ i.i.d as ${\bf X}|{\bf Y}$. Then, we expand the non-linear
function $\mathrm{f}$ around $x_{0}$ chosen as an upper bound of ${\bf X}|{\bf
Y}$. Finally, the infinite power series is statistically formulated thanks to
the discrete random variable $N$ the order of Taylor expansion:

$N_{\bf Y}=B_{{\bf X}_{1}|{\bf Y}}+\sum_{i=2}^{+\infty}\,i\,B_{{\bf
X}_{i}|{\bf Y}}\,\prod_{q=1}^{i-1}(1-B_{{\bf X}_{q}|{\bf Y}})$ (21)

where the $B_{{\bf X}_{i}|{\bf Y}}$ are Bernoulli random variables:

$\begin{array}[]{ll}B_{{\bf X}_{i}|{\bf
Y}}=\left\\{\begin{array}[]{l}\displaystyle 1\hskip 28.45274pt\mbox{with
probability}\hskip 11.38109ptP_{{\bf X}_{i}|{\bf Y}}=\frac{{\bf X}_{i}|{\bf
Y}}{x_{0}}\vspace{3mm}\\\ \displaystyle 0\hskip 28.45274pt\mbox{with
probability}\hskip 11.38109pt1-P_{{\bf X}_{i}|{\bf
Y}}\end{array}\right.\end{array}$ (22)

In the end, the reformulation is:

$C={\cal E}_{{\bf Y},({\bf X}_{1},{\bf X}_{2},...,{\bf X}_{{\bf N}_{\bf
Y}})|{\bf Y}}\left(w(N_{\bf Y})\right)$ (23)

where the Monte Carlo weight function $w$ is

$w(u)=\left\\{\begin{array}[]{l}\displaystyle\exp\left(C_{4}\,x_{0}^{-\alpha}\right)\displaystyle\sum_{m=1}^{3}C_{m}\,x_{0}^{(m\,\alpha-1)}\;\;\mbox{if}\;u=1\vspace{6mm}\\\
w(1)+\exp\left(C_{4}\,x_{0}^{-\alpha}\right)\displaystyle\sum_{k=0}^{u}(-1)^{k}\displaystyle\sum_{p=0}^{k}\frac{1}{p!(k-p)!}\displaystyle\sum_{m=1}^{3}C_{m}\,x_{0}^{(m\alpha-1+p)}\\\
\hskip
85.35826pt\times\displaystyle\prod_{j=1}^{k-p}(m\alpha-j)\displaystyle\sum_{r=1}^{p}(-1)^{r}\,C_{4}\,x_{0}^{-(\alpha+r)}\prod_{l=0}^{r-1}(\alpha+l)\;\;\mbox{if}\;u>1\end{array}\right.$
(24)

with $u\in\mathbb{N}$ and $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$, $\alpha$ and
$x_{0}$ known constants.

###  1.3 Algorithm

Step 1

Uniform sampling of a moment ${\bf y}$ of the year and index initialisation
$n=1$.

Step 2

Sampling of the $n$-th optical path contributing to thermal-power collection
at instant ${\bf y}$ and computation of its contribution $x_{n}$ (as presented
in[5, 6]).

Step 3

Computation of the probability $P_{n}$ in Eq. 22 and sampling of a realisation
$b_{n}$ of the Bernoulli random variable $B_{{\bf X}_{n}|{\bf Y}}$.

Step 4

If $b_{n}=0$, the algorithm loops to step 2 after incrementation of $n$. Else,
the procedure is terminated and the weight ${\hat{w}}(n)$ computed according
to Eq. 24.

###  1.4 Simulated configuration

Results shown in Extended Data Figure 1 are obtained for a 1 MW solar plant
with a 80 m high central receiver tour and 160 heliostats arranged in a radial
stagered layout (nueen method). Values of the conversion-rate parameters
$C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$ and $\alpha$ are given in [3, 4]. Meteonorm
DNI database.

##  2 SI2 — Wave scattering

This example is fully detailed in Charon et al. 2016[7].

###  2.1 The problem

Here we address the solution of Schiff’s approximation[8], also known as the
anomalous diffraction approximation[9], for an incident plane wave with
propagation direction ${\bf e}_{i}$ and wave number $k$ scattered by the
scattering potential $U$ that takes value $U^{in}$ inside a large domain
compared to $1/k$ and $0$ outside (this domain corresponds to the shape of the
scatterer), with $|U^{in}|\ll k^{2}$. Let $\Gamma_{\bf Y}$ be the projected
surface of the scatterer seen from ${\bf e}_{i}$, given the scatterer
orientation ${\bf Y}$. Straight rays along ${\bf e}_{i}$ cross $\Gamma_{\bf
Y}$ at locations ${\bf P}|{\bf Y}$; the domain of definition of ${\bf P}|{\bf
Y}$ is then ${\cal D}_{{\bf P}|{\bf Y}}\equiv{\bf P}|{\bf Y}$. These rays are
attenuated and phase shifted over the crossing length $l({\bf P}|{\bf Y})$.
The complex random variable $X|{\bf Y}=X_{{\bf P}|{\bf Y}}({\bf e}_{s})$ is
the contribution of one of these rays to the scattered field in a direction
${\bf e}_{s}$ deviating from ${\bf e}_{i}$ by a small scattering angles
$\theta$[7]:

$X_{{\bf P}|{\bf Y}}({\bf e}_{s})=A_{\bf
Y}\frac{k}{2\,\pi}\exp(ik\,\theta\,{\bf b}\cdot{\bf P}|{\bf
Y})\left\\{1-\exp\left(\frac{i}{2k}U^{in}\,l({\bf P}|{\bf Y})\right)\right\\}$
(25)

where ${\bf b}$ is the unit vector along the projection of ${\bf e_{s}}$ on
the plane containing ${\cal D}_{{\bf P}|{\bf Y}}$ and $A_{\bf Y}$ is the area
of ${\cal D}_{{\bf P}|{\bf Y}}$. Given the orientation ${\bf Y}$ of the
scatterer, the conditional expectation $S_{\bf Y}({\bf e}_{s})={\cal
E}_{X|{\bf Y}}(X|{\bf Y})={\cal E}_{{\bf P}|{\bf Y}}(X_{{\bf P}|{\bf Y}}({\bf
e}_{s}))$ is the complex scattering amplitude in direction ${\bf e}_{s}$[8, 9,
10, 7]. For the scatterer orientation ${\bf Y}$, the far-field scattering
diagram is given by the _differential scattering cross-section_[9, 12]:
${\hat{W}}_{\bf Y}({\bf e}_{s})=|S_{\bf Y}({\bf e}_{s})|^{2}=|{\cal E}_{{\bf
P}|{\bf Y}}(X_{{\bf P}|{\bf Y}}({\bf e}_{s}))|^{2}$.

We address the Monte Carlo computation of $W({\bf e}_{s})={\cal E}_{\bf
Y}({\hat{W}}_{\bf Y}({\bf e}_{s}))$ which is ${\hat{W}}_{\bf Y}({\bf e}_{s})$
averaged over the statistics of the scatterer orientation ${\bf Y}$ [9, 12,
10]. The full expression is then

$\begin{array}[]{lcl}W({\bf e}_{s})&=&{\cal E}_{\bf Y}\left(|S_{\bf Y}({\bf
e}_{s})|^{2}\right)\\\ &=&{\cal E}_{\bf Y}\left(|{\cal E}_{{\bf P}|{\bf
Y}}(X_{{\bf P}|{\bf Y}}({\bf e}_{s}))|^{2}\right)\\\ &=&{\cal E}_{\bf
Y}\left(\left(\Re\,{\cal E}_{{\bf P}|{\bf Y}}(X_{{\bf P}|{\bf Y}}({\bf
e}_{s}))\right)^{2}+\left(\Im\,{\cal E}_{{\bf P}|{\bf Y}}(X_{{\bf P}|{\bf
Y}}({\bf e}_{s}))\right)^{2}\right)\end{array}$ (26)

where the configuration spaces are the orientation vectors of the scattering
potential with respect to ${\bf e}_{i}$ (${\cal D}_{\bf Y}$) and the projected
surface of the scattering potential seen from the incident direction ${\bf
e}_{i}$ (${\cal D}_{{\bf P}|{\bf Y}}$).

###  2.2 Non-Linear Monte Carlo formulation

$W({\bf e}_{s})$ is reformulated based on the definition of two independent
and identically distributed location random-variables ${\bf P}_{1}|{\bf Y}$
and ${\bf P}_{2}|{\bf Y}$:

$\begin{array}[]{lcl}W({\bf e}_{s})&=&{\cal E}_{\bf Y}\left(\Re\,{\cal
E}_{{\bf P}_{1}|{\bf Y}}(X_{{\bf P}_{1}|{\bf Y}}({\bf e}_{s}))\;\Re\,{\cal
E}_{{\bf P}_{2}|{\bf Y}}(X_{{\bf P}_{2}|{\bf Y}}({\bf e}_{s}))+\Im\,{\cal
E}_{{\bf P}_{1}|{\bf Y}}(X_{{\bf P}_{1}|{\bf Y}}({\bf e}_{s}))\;\Im\,{\cal
E}_{{\bf P}_{2}|{\bf Y}}(X_{{\bf P}_{2}|{\bf Y}}({\bf e}_{s}))\right)\\\
&=&{\cal E}_{{\bf Y},({\bf P}_{1},{\bf P}_{2})|{\bf Y}}\left(w({\bf
P}_{1},{\bf P}_{2},{\bf e}_{s})\right)\end{array}$ (27)

where the Monte Carlo weight function $w$ is:

$w({\bf P}_{1},{\bf P}_{2},{\bf e}_{s})=\Re\,X_{{\bf P}_{1}|{\bf Y}}({\bf
e}_{s})\;\Re\,X_{{\bf P}_{2}|{\bf Y}}({\bf e}_{s})+\Im\,X_{{\bf P}_{1}|{\bf
Y}}({\bf e}_{s})\;\Im\,X_{{\bf P}_{2}|{\bf Y}}({\bf e}_{s})$ (28)

with $X_{{\bf P}_{q}|{\bf Y}}({\bf e}_{s})$ defined in Eq. 25.

###  2.3 Algorithm

The sampling procedures of the Monte-Carlo algorithm are then:

Step 1

Isotropic sampling of an orientation ${\bf y}$ of the scattering potential.

Step 2

Uniform sampling of the first location ${\bf p}_{1}$ on the projected surface
defined by ${\bf y}$ and computation of the corresponding crossing length
$l({\bf p}_{1})$: the realisation $x_{1}$ of $X_{{\bf P}|{\bf Y}}({\bf
e}_{s})$ is computed according to Eq. 25.

Step 3

Uniform sampling of the second location ${\bf p}_{2}$ on the projected surface
defined by ${\bf y}$ and computation of the corresponding crossing length
$l({\bf p}_{2})$: the realisation $x_{2}$ of $X_{{\bf P}|{\bf Y}}({\bf
e}_{s})$ is computed according to Eq. 25.

Step 4

Computation of the weight ${\hat{w}}=w({\bf p}_{1},{\bf p}_{2},{\bf e}_{s})$
according to Eq. 28: ${\hat{w}}=\Re\,x_{1}\;\Re\,x_{2}+\Im\,x_{1}\;\Im\,x_{2}$

Codes for the implementation of this algorithm in the case of spheroidal and
cylindrical scattering potentials, as well as validation against reference
solution, are provided in Charon et al. 2016 [7] (see also
http://edstar.lmd.jussieu.fr/codes).

###  2.4 Simulated configuration

Results shown in Extended Data Figure 2 are obtained for $k=2\pi/\lambda$,
with wavelength $\lambda=500nm$ and the scattering potential
$U^{in}=-2k(m-1)$, with $m=1.2-i\,4\,.10^{-3}$ (scatterer relative refractive
index[10]). The shape of the scatterer is helical with length $L=50\mu m$,
pitch $P=15\mu m$, helix diameter $D=20\mu m$, cylinder diameter $d=3.5\mu m$.

##  3 SI3 — Phytoplankton growth in light-limited environments

###  3.1 The problem

Here we address the production of reference solutions for a photobioreactor
model. This model is based on a radiative transfer approach presented in[11].
The rate of photon absorption ${\bf A}_{\bf Y}$ by micro-algae at location
${\bf Y}$ within the culture volume is solution of the Radiative Transfer
Equation. The present study is based on the standard linear transport MC
algorithm presented in[15] for the estimation of ${\bf A}_{\bf Y}$ at any
location within the photobioreactor. In this algorithm, realisations of the
optical path random variable ${\bf\Gamma}(C)|{\bf Y}$ are sampled backward
from the absorption location ${\bf Y}$ to the light emitting surface (the
surface of the 979 light-diffusing optical fibres, see EDF3). Since micro-
algae are scattering visible light, ${\bf\Gamma}(C)|{\bf Y}$ depends on their
concentration $C$ within the suspension: the scattering and absorption
coefficient of the suspension are proportional to $C$ (independent
scattering). The random variable ${\bf X}(C)|{\bf Y}={\bf
X}_{{\bf\Gamma}(C)|{\bf Y}}(C)$ is the contribution of one of these optical
paths to the photon absorption rate at ${\bf Y}$ (see[15] for the detailed
expression of ${\bf X}_{{\bf\Gamma}(C)|{\bf Y}}(C)$). This contribution ${\bf
X}_{{\bf\Gamma}(C)|{\bf Y}}(C)$ depends on the biomass concentration $C$
because the absorption coefficient of the suspension determines the fraction
of incident light flux that is transmitted along the optical path
${\bf\Gamma}(C)|{\bf Y}$ (according to Beer law). In the end, the conditional
expectation ${\bf A}_{\bf Y}(C)={\cal E}_{{\bf X}|{\bf Y}}({\bf X}|{\bf
Y})={\cal E}_{{\bf\Gamma}(C)|{\bf Y}}({\bf X}_{{\bf\Gamma}(C)|{\bf Y}}(C))$ is
the rate of photon absorption at location ${\bf Y}$, for the biomass
concentration $C$.

Photons absorbed by a micro-algae are converted within the photosynthetic
units and the Z-scheme[17, 18], leading to the local biomass growth-rate[14,
11]

$\begin{array}[]{lcl}r_{\bf Y}(C)&=&\mathrm{f}({\bf A}_{\bf Y}(C),C)\\\
&=&C\left({\alpha\;\frac{K\,{\bf A}_{\bf Y}(C)}{K+{\bf A}_{\bf
Y}(C)}-\beta\;\frac{K_{r}}{K_{r}+{\bf A}_{\bf Y}(C)}}\right)\\\
&=&C\left({\alpha\;\frac{K\,{\cal E}_{{\bf\Gamma}(C)|{\bf Y}}({\bf
X}_{{\bf\Gamma}(C)|{\bf Y}}(C))}{K+{\cal E}_{{\bf\Gamma}(C)|{\bf Y}}({\bf
X}_{{\bf\Gamma}(C)|{\bf Y}}(C))}-\beta\;\frac{K_{r}}{K_{r}+{\cal
E}_{{\bf\Gamma}(C)|{\bf Y}}({\bf X}_{{\bf\Gamma}(C)|{\bf
Y}}(C))}}\right)\end{array}$ (29)

where $\alpha$, $\beta$, $K$ and $K_{r}$ are constant parameters that depend
on the studied microorganism (see[14] for values in the case of Chlamydomonas
Reinhardtii).

Due to light absorption and scattering by micro-algae, the field ${\bf A}_{\bf
Y}$ is heterogeneous within the volume (less light farther from the fibres and
for higher concentrations) and so is the local photosynthetic rate $r_{\bf
Y}(C)$. We address the Monte Carlo computation of $R(C)={\cal E}_{\bf
Y}(r_{\bf Y}(C))$, the local photosynthetic rate $r_{\bf Y}(C)$ averaged over
locations ${\bf Y}$ of the culture volume. The full expression is then

$R(C)={\cal E}_{\bf Y}\left(C\left({\alpha\;\frac{K\,{\cal
E}_{{\bf\Gamma}(C)|{\bf Y}}({\bf X}_{{\bf\Gamma}(C)|{\bf Y}}(C))}{K+{\cal
E}_{{\bf\Gamma}(C)|{\bf Y}}({\bf X}_{{\bf\Gamma}(C)|{\bf
Y}}(C))}-\beta\;\frac{K_{r}}{K_{r}+{\cal E}_{{\bf\Gamma}(C)|{\bf Y}}({\bf
X}_{{\bf\Gamma}(C)|{\bf Y}}(C))}}\right)\right)$ (30)

###  3.2 Non-Linear Monte Carlo formulation

$R(C)$ is reformulated based on the Taylor expansion of $\mathrm{f}$ around
$x_{0}$. First, independent and equally distributed optical-path random
variables ${\bf\Gamma}_{i}(C)|{\bf Y}$ are introduced, defining independent
and equally distributed contributions ${\bf X}_{i}(C)|{\bf Y}={\bf
X}_{{\bf\Gamma}_{i}(C)|{\bf Y}}(C)$. Then, we expand the non-linear function
in Eq. 29 around $x_{0}$ chosen as an upper bound of ${\bf X}_{i}(C)|{\bf Y}$.
Finally, the infinite power series is statistically formulated thanks to the
discrete random variable $N$ the order of Taylor expansion:

$N_{\bf Y}(C)=B_{{\bf X}_{1}(C)|{\bf Y}}+\sum_{i=2}^{+\infty}\,i\,B_{{\bf
X}_{i}(C)|{\bf Y}}\,\prod_{q=1}^{i-1}(1-B_{{\bf X}_{q}(C)|{\bf Y}})$ (31)

where the $B_{{\bf X}_{i}(C)|{\bf Y}}$ are Bernoulli random variables:

$\begin{array}[]{ll}B_{{\bf X}_{i}(C)|{\bf
Y}}=\left\\{\begin{array}[]{l}\displaystyle 1\hskip 28.45274pt\mbox{with
probability}\hskip 11.38109ptP_{{\bf X}_{i}(C)|{\bf Y}}=\frac{{\bf
X}_{i}(C)|{\bf Y}}{x_{0}}\vspace{3mm}\\\ \displaystyle 0\hskip
28.45274pt\mbox{with probability}\hskip 11.38109pt1-P_{{\bf X}_{i}(C)|{\bf
Y}}\end{array}\right.\end{array}$ (32)

In the end, the reformulation is:

$R(C)={\cal E}_{{\bf Y},({\bf\Gamma}_{1},{\bf\Gamma}_{2},...,{\bf\Gamma}_{{\bf
N}_{\bf Y}})|{\bf Y}}\left(w(N_{\bf Y})\right)$ (33)

where the Monte Carlo weight function $w$ is

$w(u)=\left\\{\begin{array}[]{l}\displaystyle
C\left[\alpha\;\frac{K\,x_{0}}{K+x_{0}}-\beta\;\frac{K_{r}}{K_{r}+x_{0}}\right]\;\;\mbox{if}\;u=1\vspace{6mm}\\\
w(1)-\displaystyle\sum_{i=1}^{u-1}C\left[\alpha\;\frac{K^{2}}{K+x_{0}}\left(\frac{x_{0}}{K+x_{0}}\right)^{i}+\beta\;\frac{K_{r}}{K_{r}+x_{0}}\left(\frac{x_{0}}{K+x_{0}}\right)^{i}\;\right]\;\;\mbox{if}\;u>1\end{array}\right.$
(34)

with $u\in\mathbb{N}$ and $\alpha$, $K$, $\beta$, $K_{r}$, $x_{0}$ and $C$
known constants.

###  3.3 Algorithm

Step 1

Uniform sampling of a location ${\bf y}$ within the culture volume and index
initialisation $n=1$.

Step 2

Sampling of the $n$-th optical path contributing to absorption at ${\bf y}$:
the realisation $x_{n}$ of ${\bf X}_{{\bf\Gamma}_{n}(C)|{\bf Y}}(C)$ is
computed according to[15].

Step 3

Computation of the probability in Eq. 32, $P=x_{n}/x_{0}$, and sampling of a
realisation $b_{n}$ of the Bernoulli random variable $B_{{\bf X}_{n}(C)|{\bf
Y}}$.

Step 4

If $b_{n}=0$, the algorithm loops to step 2 after incrementation of $n$. Else,
the procedure is terminated and the weight ${\bf w}(n)$ computed according to
Eq. 34.

###  3.4 Simulated configuration

Results shown in Extended Data Figure 3 are obtained for kinetics parameters
$\alpha=1.785\,.10^{-9}$ $kg/\mu mol_{h\nu}$, $\beta=4.057\,.10^{-6}$
$s^{-1}$, $K=32000$ $\mu mol_{h\nu}/kg/s$ and $K_{r}=7500$ $\mu
mol_{h\nu}/kg/s$. The photobioreactor is a 25L DiCoFluV[13]: reactor diameter
$16.5\,cm$, reactor height $1\,m$, optical fibres diameter $1.2\,mm$, fibres
height $1\,m$, $979$ fibres, hexagonal lattice fibre arrangement with distance
$4.8\,mm$ between two fibres axis (see[15, 13]), homogeneous surface flux
density $25\,\mu mol_{h\nu}/m^{2}/s$ emitted at fibres surface. The radiative
properties (absorption and scattering by micro-algae) are those presented
in[10] for Chlamydomonas Reinhardtii.

##  4 SI4 — Atmospheric radiative transfer: top-of-atmosphere specific
intensity (from earth toward the outer space)

This example is fully detailed in Galtier et al. 2015[19].

###  4.1 The problem

Photons are emitted either by a surface (here the ground) or by the volume
(the atmosphere). The radiation $I(\Delta\nu)$ perceived at the observation
location TOA is the fraction of all these photons that have a frequency inside
the observation-band $\Delta\nu$ and reach the sensor, i.e. which are not
absorbed by the atmosphere in between. When no scattering occurs (clear sky),
photon-paths are straight lines and the fraction of the photons of frequency
$\nu_{y}$ that travel (let say along the strict vertical) from altitude
$H_{y}$ to TOA is

$exp(-\tau_{y})$

where

$\tau_{y}=\int_{H_{y}}^{TOA}dx\sum_{j_{x}=1}^{N_{t}}h_{\nu_{y},j_{x}}(x)$

This exponential extinction is that of Beer law where $\tau_{y}$ is the
monochromatic optical thickness, i.e. the sum of all transition cross-sections
$h_{\nu_{y},j_{x}}(x)$ at all intermediate altitudes $x$. In this example, the
random variable $\bf Y$ is a vector that gathers the altitude $H_{\bf Y}$ of
emission, the photon-path $\bf\Gamma_{Y}$ from $H_{\bf Y}$ to TOA and
frequency $\nu_{\bf Y}$ of emission. Similarly, $\bf X$ gathers all the
description of an absorption event, altitude $H_{\bf X}$ and index of
absorption-transition $j_{\bf X}$.

In standard Monte Carlo approaches, $\tau_{y}$ is either precomputed and
tabulated, or is easily computable from tabulated values of the absorption
coefficient ($k=\sum_{j_{x}=1}^{N_{t}}h_{\nu_{y},j_{x}}$). Then the Monte
Carlo algorithm deals only with the sampling of $\bf Y$ and the weight
function has the form

$w(y)={\cal S}_{y}exp(-\tau_{y})$

${\cal S}_{\bf Y}$ is the source associated with the emission $\bf Y$
(Planck’s function at the local temperature times an emissive power that
depends on concentration, pressure and temperature). The sum over the $N_{t}$
transitions $j_{x}$ is thus not handled by standard Monte Carlo algorithms,
despite the fact that this sum has all the features inviting to make use of
statistical approaches: $N_{t}$ is huge, typically of the order of $10^{6}$,
and the deterministic pre-calculations can be computationally very demanding
(and are to be re-performed when testing each new spectroscopic assumption).
No attempt has been made so far to address this sum statistically, together
with the photon-path statistics because these two statistics are combined via
the non-linearity of the exponential extinction.

###  4.2 Non-Linear Monte Carlo formulation

So the present question is to design a Monte Carlo algorithm performing the
summation $\sum_{j=1}^{N_{t}}h_{\nu,j}$ together with the integrals over
frequency, emission location and absorption location, despite of the
nonlinearity of Beer extinction. This starts by looking at $\tau_{y}$ in
statistical terms : the integral over altitudes $H_{\bf X}$ is transformed
into the sampling of absorption locations, and the sum
$\sum_{j_{x}=1}^{N_{t}}h_{\nu_{y},j_{x}}$ into the sampling of absorption-
transitions. This double sampling is formally summarised into the random
variable $\bf X$. To get a sample $x$ of $\bf X$, one first samples $H_{\bf
X}$ and $j_{\bf X}$. Then $x$ is the value that the optical thickness
$\tau_{y}$ would have if the atmosphere was homogeneous at the thermodynamic
conditions of altitude $H_{x}$ and if all transitions were identical to
$j_{x}$. This random variable depends on $\bf Y$ and its expectation, knowing
${\bf Y}=y$, is ${\cal E}({\bf X}|{\bf Y}=y)=\tau_{y}$. We therefore get
$I\left(\Delta\nu\right)={\cal E}_{\bf Y}\left({\cal S}_{\bf
Y}\;exp\right(-{\cal E}_{{\bf X}|{\bf Y}}({\bf X}|{\bf Y})\left)\right)$

###  4.3 Algorithm

A strict application of the solution described in Methods implies to first
sample $\bf Y$ (as in a standard Monte Carlo algorithm), then sample the
degree $n$ of one monomial in the Taylor expansion of the negative-exponential
function $\mathrm{f}$, and draw $n$ independent samples $x_{1},x_{2}...x_{n}$
of $\bf X$, i.e $n$ absorption locations $H_{x_{q}}$ (between emission
location and TOA) and $n$ transition indexes $j_{x_{q}}$:

Step 1

sample a frequency,

Step 2

sample an emission-altitude,

Step 3

sample an order $n$ of the development of the exponential,

Step 4

sample $n$ paired absorption-altitudes, transitions.

###  4.4 Null-collision reformulation

Instead of implementing this solution in this straightforward manner, we chose
to use null-collisions[20]. The same quantities are sampled, but the order is
different, leading to quite intuitive physical pictures. The first step is
still to sample a frequency, but then we implement a backward tracking
algorithm, starting from the observation altitude. We sample a first
collision-altitude as if photons were emitted at the observation location in
the downward direction in an homogeneous atmosphere. Then we sample a
transition and a statistical test is made to determine whether the algorithm
stops at this altitude or continues (a Bernoulli test to determine whether the
collision is a true one or a null-collision). If it continues, this first
collision-altitude and first transition have the status of one of the $n$
absorption-altitudes and transitions in the above presented algorithm. Then a
next collision-location is sampled, together with a next transition, etc. When
the algorithm stops after $n$ collisions, the final altitude is interpreted as
the emission-altitude. So $n$ and the emission-altitude are not sampled first:
they are sampled by the successive Bernoulli tests as in a standard backward-
tracking multiple-scattering algorithm.

###  4.5 The final algorithm

Step 1

Initialisation of the current altitude at $H_{x}\leftarrow TOA$.

Step 2

Uniform sampling of a frequency $\nu_{y}$ in the considered infrared band (or
the whole infrared).

Step 3

Exponential sampling of a travelled distance $d$ before absorption for a
photon of frequency $\nu_{y}$ travelling from $H_{x}$ in the backward
direction within a virtual homogeneous atmosphere of absorption coefficient
$\hat{k}_{\nu_{y}}$.

Step 4

If the travelled distance leads the photon to hit the surface,
$H_{x}\leftarrow 0$ and the algorithm jumps to Step 7. Otherwise
$H_{x}\leftarrow H_{x}-d$.

Step 5

Sampling of a state-transition $j_{x}$ according to ${\cal P}(1),{\cal
P}(2)...{\cal P}(N_{t})$ (see[19]).

Step 6

Bernoulli trial of probability
$P_{a}=\frac{h_{\nu_{y},j_{x}}}{\hat{k}_{\nu_{y}}{\cal P}(j_{x})}$ to decide
whether the algorithm jumps to Step 7 or Step 3.

Step 7

$H_{y}\leftarrow H_{x}$ and the algorithm stops with $\hat{w}=\Delta\nu B_{y}$
where $B_{y}$ is the Planck function at frequency $\nu_{y}$ for the
atmospheric temperature at $H_{y}$.

###  4.6 Simulated configuration

The simulation results in Extended Data Fig. 4 have been obtained using the
HITRAN spectroscopic database. All other parameters are given in Galtier et
al. 2015[19].

###  4.7 Computational performance

As far as computational costs are concerned, our computations did not require
to first scan HITRAN and pre-compute absorption coefficients at all altitudes
and all frequencies before running the Monte Carlo code for transfer since our
Monte Carlo code handles both simultaneously despite the nonlinearity of the
exponential extinction. The computational benefit is very significant when
studying the effects of new spectroscopic data or new line-shape assumptions.
Otherwise, for earth applications and fixed spectroscopic assumptions, there
is no problem associated to the pre-computation of absorption coefficients and
to their tabulation as function of molecular composition and thermodynamic
state. But for combustion or astrophysics applications, thinking in particular
of exoplanets, the diversity of compositions is extremely wide, temperatures
can by high, imply the use of much larger spectroscopic databases, including
hot lines, and pre-computation/tabulation of absorption coefficients is a
today’s challenge in itself. Our nonlinear Monte Carlo suppresses this
need[19].

##  5 SI5 — Gas kinetics

###  5.1 The problem

We consider a gas of interacting particles, with collisions following Maxwell
model, and a cross section
$\sigma(\vec{c},\vec{c^{*}})=\frac{\kappa}{4\pi\|\vec{c}-\vec{c^{*}}\|}$.
Particles are confined within an harmonic static trap of pulsation $\omega$ so
that acceleration at position $\vec{r}$ is spring-like:
$\vec{a}(\vec{r})=-\omega^{2}\vec{r}$.

We follow the distribution function $f(\vec{r},\vec{c},t)$ at location
$\vec{r}$, velocity $\vec{c}$ at time $t$. We consider it known at some time
$t_{I}$, either constrained at local equilibrium, according to:

$f_{{}_{\bf LEQ}}/n_{I}=p_{{\cal N}(\vec{u}_{I},c_{q,I}^{2})}$ (35)

or, out of equilibrium, according to:

$f_{{}_{\bf
BKW}}/n_{I}=\frac{1}{3}\frac{(\vec{c}-\vec{u}_{I})^{2}}{\frac{c_{q,I}^{2}}{5/3}}\
p_{{\cal N}(\vec{u}_{I},\frac{c_{q,I}^{2}}{5/3})}$ (36)

where $n_{I}$ denotes density, $\vec{u}_{I}$ denotes the mean velocity and
$c_{q,I}$ denotes the mean square speed at $\vec{r}$ at time $t_{I}$ ;
$p_{{\cal N}(\mu,\sigma^{2})}$ denotes probability density function of a
Gaussian random variable, with expectation $\mu$ and variance $\sigma^{2}$.

#### For the case in fig. 3a,

we set $\omega=0$, and $n_{I}$, $\vec{u}_{I}$ and $c_{q,I}$ are set
homogeneous, corresponding to the case studied par Krook and Wu[21, 22].
Starting from $f_{{}_{\bf BKW}}$, the gas relaxes to equilibrium $f_{{}_{\bf
LEQ}}$.

Adimensional time is $\kappa n_{I}t$ and adimensional distribution function is
$\frac{1}{n_{I}}\left(\sqrt{2\pi}c_{q,I}\right)^{3}f$.

#### For the case in fig. 3b,

we set $\omega=2\pi$, and $n_{I}$, $\vec{u}_{I}$ and $c_{q,I}$ are set
heterogeneous to correspond to one state of the undamped oscillation
(_breathing mode_[23]):

$\begin{array}[]{lcl}n_{I}&=&p_{{\cal
N}(\vec{0},\frac{c_{q,EQ}^{2}}{\omega^{2}})}\vspace{3mm}\\\
\vec{u}_{I}&=&\epsilon\omega\vec{r}\vspace{3mm}\\\
c_{q,I}&=&\sqrt{1-\epsilon^{2}}\ c_{q,EQ}\end{array}$ (37)

with $\epsilon=\frac{\Delta c_{q}^{2}}{c_{q,EQ}^{2}}$ where $\Delta c_{q}^{2}$
is the maximal deviation of $c_{q}^{2}$ from its equilibrium value
$c_{q,EQ}^{2}$.

With these initial values, and starting from local equilibrium $f_{{}_{\bf
LEQ}}$, the local equilibrium remains unbroken and the gas displays the
undamped oscillation for any value of the cross section. Conversely, starting
from out of equilibrium $f_{{}_{\bf BKW}}$, any positive value of cross
section will lead to dampened oscillations.

Adimensional time is $\frac{\omega}{2\pi}t$, adimensional distribution
function is $\left(2\pi\frac{c_{q,EQ}^{2}}{\omega}\right)^{3}f$ and
adimensional parameter for cross-section is
$\frac{\omega^{2}}{c_{q,EQ}^{3}}\kappa$.

Simulation results are given for $\epsilon=0.2$ and
$\frac{\omega^{2}}{c_{q,EQ}^{3}}\kappa=3$.

###  5.2 Non-Linear Monte Carlo formulation

The distribution function obeys the Boltzmann equation:

$\displaystyle\frac{\partial f}{\partial t}+$
$\displaystyle\vec{c}\cdot\vec{\text{grad}}_{\cal R}(f)+\text{div}_{\cal
C}(f\vec{a})$ (38) $\displaystyle=$ $\displaystyle\int_{\cal
C}d\vec{c}_{*}\int_{4\pi}d\vec{u}\ \left\|\vec{c}-\vec{c}_{*}\right\|\
\sigma(\vec{c},\vec{c}_{*})\left(-ff^{*}+f^{\alpha}f^{\beta}\right)$

where $f\equiv f(\vec{r},\vec{c},t)$, $f^{*}\equiv f(\vec{r},\vec{c}_{*},t)$,
$f^{\alpha}\equiv f(\vec{r},\vec{c}_{\alpha},t)$ and $f^{\beta}\equiv
f(\vec{r},\vec{c}_{\beta},t)$.

Velocities $\vec{c}_{\alpha}$ and $\vec{c}_{\beta}$ are functions of
$\vec{c}$, $\vec{c}_{*}$ and $\vec{u}$ :
$\vec{c}_{\alpha}=\frac{1}{2}(\vec{c}+\vec{c}_{*}+\left\|\vec{c}-\vec{c}_{*}\right\|\vec{u})$,
$\vec{c}_{\beta}=\frac{1}{2}(\vec{c}+\vec{c}_{*}-\left\|\vec{c}-\vec{c}_{*}\right\|\vec{u})$.

The usual PDE expression of this model, given above, can be translated into
its Fredholm counterpart, following:

$\displaystyle{f(\vec{r},\vec{c},t)}=$
$\displaystyle\int_{-\infty}^{t}dt^{\prime}\ \hat{\nu}(t^{\prime})\
\exp\left(-\int_{t^{\prime}}^{t}dt^{\prime\prime}\
\hat{\nu}(t^{\prime\prime})\right)\times$ (39) $\displaystyle\Bigg{[}{\cal
H}(t_{I}-t^{\prime})\;{f(\vec{r}_{b}(t_{I}),\vec{c}_{b}(t_{I}),t_{I})}$
$\displaystyle+{\cal
H}(t^{\prime}-t_{I})\left(1-\frac{{\nu(t^{\prime})}}{\hat{\nu}(t^{\prime})}\right){f\left(\vec{r}_{b}(t^{\prime}),\vec{c}_{b}(t^{\prime}),t^{\prime}\right)}+\frac{{s(t^{\prime})}}{\hat{\nu}(t^{\prime})}\Bigg{]}$

with

$\displaystyle\nu(t^{\prime})=\int_{\cal C}d\vec{c}_{*}\int_{4\pi}d\vec{u}$
$\displaystyle\left\|\vec{c}_{b}(t^{\prime})-\vec{c}_{*}\right\|\sigma(\vec{c}_{b}(t^{\prime}),\vec{c}_{*}){f(\vec{r}_{b}(t^{\prime}),\vec{c}_{*},t^{\prime})}$
(40) $\displaystyle s(t^{\prime})=\int_{\cal
C}d\vec{c}_{*}\int_{4\pi}d\vec{u}$
$\displaystyle\left\|\vec{c}_{b}(t^{\prime})-\vec{c}_{*}\right\|\sigma(\vec{c}_{b}(t^{\prime}),\vec{c}_{*}){f(\vec{r}_{b}(t^{\prime}),\vec{c}_{\alpha}(t^{\prime}),t^{\prime})}{f(\vec{r}_{b}(t^{\prime}),\vec{c}_{\beta}(t^{\prime}),t^{\prime})}$

where $\vec{r}_{b}$ and $\vec{c}_{b}$ are location and velocity corresponding
to the ballistic path through $\vec{r}$ at time $t$ with velocity $\vec{c}$,
such as:

$\left\\{\begin{array}[]{l}\partial_{t^{\prime}}\vec{r}_{b}(t^{\prime})=\vec{c}_{b}(t^{\prime})\\\
\partial_{t^{\prime}}\vec{c}_{b}(t^{\prime})=\vec{a}(\vec{r}_{b}(t^{\prime}))\end{array}\right.$

and $\vec{r}_{b}(t)=\vec{r}$, $\vec{c}_{b}(t)=\vec{c}$.

The product
${f(\vec{r}_{b}(t^{\prime}),\vec{c}_{\alpha}(t^{\prime}),t^{\prime})}{f(\vec{r}_{b}(t^{\prime}),\vec{c_{\beta}}(t^{\prime}),t^{\prime})}$
has been treated following the NLMC expansion exposed in Methods.

The exponential term is handled using null collisions technique, as exposed in
SI4. The upper bound for cross section is set at the value it takes at the
center of the gas cloud at time of maximal contraction of the undamped
oscillation :

$\hat{\nu}=\frac{\omega^{3}}{\left(2\pi
c_{q,EQ}^{2}(1-\epsilon)\right)^{3/2}}\ \kappa$

###  5.3 Algorithm

According to eq. 39, three random variables are defined for the Monte Carlo
estimates: $T^{\prime}$ for $t^{\prime}$, $\vec{U}$ for $\vec{u}$ and
$\vec{C}_{*}$ for $\vec{c}_{*}$, with:

$\begin{array}[]{lcll}p_{T^{\prime}}(t^{\prime})&=&\hat{\nu}\
\exp\left(-\hat{\nu}(t-t^{\prime})\right)&\text{over}\;]-\infty,t]\\\
p_{\vec{U}}(\vec{u})&=&\frac{1}{4\pi}&\text{over the unit sphere}\\\
p_{\vec{C^{*}}}(\vec{c}_{*})&=&p_{{\cal
N}(\vec{u}_{Boltz},c_{q,Boltz}^{2})}&\text{over velocity space}\end{array}$
(41)

where $\vec{u}_{Boltz}$ et $c_{q,Boltz}$ are respectively the mean velocity
and the mean square speed at the considered time, and set to the values
predicted in the case of undamped oscillation.

${f(\vec{r},\vec{c},t)}$ can be estimated by this recursive algorithm:

Initialisation

Sample a date $t^{\prime}$.

Recursion termination:

If $t^{\prime}\leqslant t_{I}$ return ${f(\vec{x}_{I},\vec{c}_{I},t_{I})}$.

Recursion

Sample a velocity $\vec{c}_{*}$ and a unit vector $\vec{u}$.

Compute the ballistic solution $\vec{x}_{b}(t^{\prime})$ and
$\vec{c}_{b}(t^{\prime})$.

Estimate
${\sf{f}_{*}}\leftarrow{f(\vec{x}_{b}(t^{\prime}),\vec{c}_{*},t^{\prime})}$.

Set

$\begin{array}[]{lcl}{\sf
n}&\leftarrow&{\left(\left\|\vec{c}_{b}(t^{\prime})-\vec{c}_{*}\right\|\
\sigma(\vec{c}_{b}(t^{\prime}),\vec{c}_{*},\vec{u})\
{\sf{f}_{*}}\right)}\Big{/}{\left(p_{\vec{C}_{*}}(\vec{c}_{*})p_{\vec{U}}(\vec{u})\right)\vspace{3mm}}\\\
{\sf Q}&\leftarrow&{\sf n}/{\hat{\nu}}\end{array}$ (42)

If ${\sf Q}\in[0,1]$, then set ${\sf P}\leftarrow{\sf Q}$ else set ${\sf
P}\leftarrow\frac{{\sf Q}}{2{\sf Q}-1}$.

Sample a standard uniform $r$.

If $r>{\sf P}$, then estimate
${\sf{f}_{b}}\leftarrow{f(\vec{x}_{b}(t^{\prime}),\vec{c}_{b}(t^{\prime}),t^{\prime})}$
and return $\frac{1-{\sf Q}}{1-{\sf P}}\ {\sf{f}_{b}}$,

Else set

$\begin{array}[]{lcl}\vec{c}_{\alpha}(t^{\prime})&\leftarrow&\frac{1}{2}(\vec{c}_{b}(t^{\prime})+\vec{c}_{*}+\left\|\vec{c}_{b}(t^{\prime})-\vec{c}_{*}\right\|\vec{u})\\\
\vec{c}_{\beta}(t^{\prime})&\leftarrow&\frac{1}{2}(\vec{c_{b}(t^{\prime})}+\vec{c}_{*}-\left\|\vec{c}_{b}(t^{\prime})-\vec{c}_{*}\right\|\vec{u})\end{array}$
(43)

Estimate
${\sf{f}}_{\alpha}\leftarrow{f(\vec{x}_{b}(t^{\prime}),\vec{c}_{\alpha}(t^{\prime}),t^{\prime})}$.

Estimate
${\sf{f}}_{\beta}\leftarrow{f(\vec{x}_{b}(t^{\prime}),\vec{c}_{\beta}(t^{\prime}),t^{\prime})}$.

Return $\left({\sf Q}/{\sf P}\right)\left({{\sf{f}}_{\alpha}\
{\sf{f}}_{\beta}}/{\sf{f}_{*}}\right)$

## References

  * [1] Steinfeld, A. Solar thermochemical production of hydrogen–a review Sol. Energy. 78, 603 (2005)
  * [2] Müller, R. et al. Transient heat transfer in a directly-irradiated solar chemical reactor for the thermal dissociation of ZnO. Appl. Therm. Eng. 28, 524 (2008)
  * [3] Schunk, L.O. et al. Heat transfer model of a solar receiver-reactor for the thermal dissociation of ZnO-Experimental validation at 10 kW and scale-up to 1 MW. Chem. Eng. J. 150 502 (2009)
  * [4] Pitz-Paal, R. et al. Heliostat field layout optimization for high-temperature solar thermochemical processing. Sol. Energy. 85, 334 (2011)
  * [5] Delatorre, J. et al. Monte Carlo advances and concentrated solar applications. Sol. Energy. 103, 653 (2014)
  * [6] Farges, O. et al. Life-time integration using Monte Carlo Methods when optimizing the design of concentrated solar power plants. Sol. Energy 113, 57-62 (2015).
  * [7] Charon, J. et al. Monte Carlo implementation of Schiff’s approximation for estimating radiative properties of homogeneous, simple-shaped and optically soft particles: Application to photosynthetic micro-organisms. J. Quant. Spectrosc. Radiat. Transf. 172, 3-23 (2016).
  * [8] Schiff, L. I. Approximation Method for High-Energy Potential Scattering. Physical Review. (1956)
  * [9] van de Hulst, H. C. Light scattering by small particles. Dover Publication, Inc.. (1981)
  * [10] Dauchet, J. and al. Calculation of the radiative properties of photosynthetic microorganisms. J. Quant. Spectrosc. Radiat. Transfer. (2015)
  * [11] Dauchet, J., Cornet, J.-F., Gros, F., Roudet, M. and Dussap, C.-G. Chapter One — Photobioreactor Modeling and Radiative Transfer Analysis for Engineering Purposes. Adv. Chem. Eng. 48, 1-106 (2016).
  * [12] Bohren, C. F. and Huffman, D. R. Absorption and scattering of light by small particles. New York, Wiley-Interscience. (1983)
  * [13] Cornet J.F. Calculation of optimal design and ideal productivities of volumetrically lightened photobioreactors using the constructal approach. Chem. Eng. Sci. 65 985 (2010)
  * [14] Takache, H. et al. Kinetic modeling of the photosynthetic growth of Chlamydomonas reinhardtii in a photobioreactor. Biotechnol. Prog. 28 681 (2012)
  * [15] Dauchet, J. et al. The practice of recent radiative transfer Monte Carlo advances and its contribution to the field of microorganisms cultivation in photobioreactors. J. Quant. Spectrosc. Radiat. Transfer. 128 52 (2013)
  * [16] Delatorre, J. et al. Monte Carlo advances and concentrated solar applications. Sol. Energy. 103, 653 (2014)
  * [17] Hatch, M.D. et al. Photosynthesis and photo-respiration. Wiley-Interscience (1971)
  * [18] Cornet J.F. & Dussap C.G. A Simple and reliable formula for assessment of maximum volumetric productivities in photobioreactors Biotechnol. Prog. 25 424 (2009)
  * [19] Galtier, M. et al. Radiative transfer and spectroscopic databases: a line-sampling Monte Carlo approach. J. Quant. Spectrosc. Radiat. Transf. (2015).
  * [20] Galtier, M. et al. Integral formulation of null-collision Monte Carlo algorithms. J. Quant. Spectrosc. Radiat. Transf. 125, 57–68 (2013).
  * [21] Krook, M. & Wu, T. T. Formation of Maxwellian Tails. Phys. Rev. Lett. 36, 1107-1109 (1976).
  * [22] Krook, M. & Wu, T. T. Exact solutions of the Boltzmann equation. Phys. Fluids 20, 1589-1595 (1977).
  * [23] Guéry-Odelin, D., Muga, J. G., Ruiz-Montero, M. J. & Trizac, E. Nonequilibrium Solutions of the Boltzmann Equation under the Action of an External Force. Phys. Rev. Lett. 112, 180602 (2014).

## Extended Data Figures