Mégane Bati
email: megane.bati@irit.fr

Richard Fournier
email: richard.fournier@laplace.univ-tlse.fr

Stéphane Blanco
email: stephane.blanco@laplace.univ-tlse.fr

Jacques Gautrais
email: jacques.gautrais@univ-tlse3.fr

Nicolas Mellado
email: nicolas.mellado@irit.fr

Mathias Paulin
email: mathias.paulin@irit.fr

Christophe Coustet
email: christophe.coustet@meso-star.com

Vincent Eymet
email: vincent.eymet@meso-star.com

Vincent Forest
email: vincent.forest@meso-star.com

Benjamin Piaud
email: benjamin.piaud@meso-star.com.

Coupling Conduction, Convection and Radiative Transfer in a Single Path-Space: Application to Infrared Rendering

Keywords: Computing methodologies → Simulation types and techniques, Rendering, Physical simulation, • Applied computing → Physics, Mathematics and statistics

lack of coupled Monte Carlo algorithms allowing to solve all the physics at once, in the same path space, rather than combining several independent MC estimators, a combination that would make the global solver critically sensitive to the complexity of each simulation space. This brings to our proposal: a coupled, single path-space, Monte Carlo algorithm for efficient multi-physics problems solving.

In this work, we combine our understanding and knowledge of Physics and Computer Graphics to demonstrate how to formulate and arrange different simulation spaces into a single path space. We define a tractable formalism for coupled heat transfer simulation using Monte Carlo, and we leverage the path-space construction to interactively compute multiple simulations with different conditions in the same scene, in terms of boundary conditions and observation time. We validate our proposal in the context of infrared rendering with different thermal simulation scenarios: e.g., room temperature simulation, visualization of heat paths within materials (detection of thermal bridges), heat diffusion capacity of thermal exchanger. We expect that our theoretical framework will foster collaboration and multidisciplinary studies. The perspectives this framework opens are detailed and we suggest a research agenda towards the resolution of coupled PDEs at the interface of Physics and Computer Graphics.

In the past decades, Monte Carlo methods have shown their ability to solve PDEs, independently of the dimensionality of the integration domain and for different use-cases (e.g. light transport, geometry processing, physics simulation). Specifically, the path-space formulation of transport equations is a key ingredient to define tractable and scalable solvers, and we observe nowadays a strong interest in the definition of simulation systems based on Monte Carlo algorithms. We also observe that, when simulating combined physics (e.g. thermal rendering from a heat transfer simulation), there is a

INTRODUCTION

In the past decades, the simulation of radiative transfer has been extensively studied in Computer Graphics in the context of Physically Based Rendering (PBR), leading to the development of powerful models, algorithms and rendering systems [START_REF] Wenzel | Mitsuba 3 renderer[END_REF][START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]]. Such systems are nowadays tailored to cope with the ever-increasing complexity of the scenes (e.g. in terms of geometry and materials) and the simulation needs (e.g. spectral complexity, coupling several physics). This is achieved by leveraging efficient data structures and data models to manage and access the scene data, and also by exploiting the scalability of Monte Carlo methods to solve the Radiative Transfer Equation in path-space.

In the present work, we advocate that these rendering approaches, which are originally motivated by movies and games use-cases, will lead to potential breakthroughs in other fundamental fields, as illustrated in climate science by [START_REF] Villefranque | The "teapot in a city": A paradigm shift in urban climate modeling[END_REF]. We propose to demonstrate this potential for the simulation of heat transfer, a process that involves radiative transfer (similarly to standard rendering), but also conduction and convection. In order to couple these three phenomena, we build upon the proposal of Tregan et al. [2023], which introduces a single coupled path-space for modeling the heat-transfer, and which is expected to be solved by Monte Carlo. Starting from this theoretical model, we formulate and demonstrate a practicable, tractable, and scalable Monte Carlo simulation framework. Unlike discretized methods or naive coupling, our framework allows to render infrared images without having to resolve the temperature field first. To generate a thermal image (see Figure 1(b)), rendering paths start from the camera, traverse the physics of coupled heat transfer in all its temporal dimensions, and lead back to the "sources": not the objects surface temperature, but the heating sources within the scene, or the boundary conditions.

A notable property of our approach is its scalability: it provides the capability to simulate coupled physical processes at very large scale with complex geometry, for instance a thermal simulation at the city scale (see Figure 4). In addition, the paths information can be collected and stored by running one single steady-state Monte Carlo simulation (Figure 1(a)), and can then be reused to simulate any transient state at no cost, e.g. to produce transient infrared movies rendering the energy propagation in capacitive objects (Figure 1(c)). The same path information can also be reused for instant rendering of any other boundary and initial conditions, providing interactive feedback (see the supplemental video). We believe that, combined with the visual analysis of the thermal paths throughout the scene, interactive thermal simulation is a game changer for use cases related to energy-transition, e.g. the rehabilitation of buildings.

From the very origin, the computer graphics community is strongly connected with physics and physicists have always heavily benefited of computer graphics. However, both communities have been historically structured upon their own agendas, and decisive advances had to first pervade one community through domain-specific communications before they became visible by the other. In this sense, even when they address closely similar questions, the two communities advance separately. For instance, physicists could import in their own field the usage of geometrical libraries and acceleration grids [START_REF] Villefranque | A Path-Tracing Monte Carlo Library for 3-D Radiative Transfer in Highly Resolved Cloudy Atmospheres[END_REF]] only once they became a major reference in computer graphics [START_REF] Wald | Embree: a kernel framework for efficient CPU ray tracing[END_REF]. In mirror, computer graphics researchers have imported null-collision algorithms only recently [START_REF] Raab | Unbiased Global Illumination with Participating Media[END_REF]] despite very early references in particle transport [START_REF] Coleman | Mathematical Verification of a Certain Monte Carlo Sampling Technique and Applications of the Technique to Radiation Transport Problems[END_REF]]. Even more recent was the convergence towards common theoretical foundations of differential light transport in computer graphics and domain deformation sensitivities in transport physics [START_REF] John | The past and future of the Monte Carlo method in thermal radiation transfer[END_REF][START_REF] Lapeyre | Monte-Carlo and sensitivity transport models for domain deformation[END_REF][START_REF] Li | Differentiable Monte Carlo Ray Tracing through Edge Sampling[END_REF][START_REF] Zhang | A Differential Theory of Radiative Transfer[END_REF][START_REF] Zhao | Physics-based differentiable rendering: from theory to implementation[END_REF].

Our approach focuses on coupling multiple physics in a single path space. The idea of coupling is precisely one of these interface questions that has recently emerged in both Computer Graphics (multiscale coupling [START_REF] Bitterli | A Position-Free Path Integral for Homogeneous Slabs and Multiple Scattering on Smith Microfacets[END_REF][START_REF] Guo | Position-Free Monte Carlo Simulation for Arbitrary Layered BSDFs[END_REF][START_REF] Heitz | Multiplescattering microfacet BSDFs with the Smith model[END_REF][START_REF] Wang | Position-Free Multiple-Bounce Computations for Smith Microfacet BSDFs[END_REF]) and Physics, in fields such as electromagnetic, solid physics, photochemistry [START_REF] Dauchet | Addressing nonlinearities in monte carlo[END_REF][START_REF] Dauchet | The practice of recent radiative transfer Monte Carlo advances and its contribution to the field of microorganisms cultivation in photobioreactors[END_REF][START_REF] Galtier | Radiative transfer and spectroscopic databases: A line-sampling Monte Carlo approach[END_REF][START_REF] Gattepaille | Integral formulations of multi-scale models for the optimization of solar photo-catalytic processes[END_REF]]. We believe that computer graphics and physics research would strongly benefit from a common theoretical framework to foster collaboration on that topic. This leads us to share our understanding and knowledge of Physics and Computer Graphics to propose the following contributions:

• building upon the proposal of [START_REF] Marc Tregan | Coupling radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm: A general theoretical framework for linear situations[END_REF], we define a tractable formalism for heat transfer simulation using Monte Carlo, and present how we leverage the path-space information to interactively compute multiple simulations with different conditions in the same scene, in terms of boundary conditions and observation time (see Section 4), • we present some results that highlight the practical significance of this proposal (see Section 5), • the perspectives we expect from this framework are detailed in Section 6 and, beyond heat transfer, we explain how its core principles open avenues towards the resolution of coupled PDEs in a broader sense.

RELATED WORKS

Radiative transfer and Physically Based Rendering. In Computer Graphics, the radiative transfer is modeled either using the rendering equation [START_REF] Kajiya | The Rendering Equation[END_REF]] on surfaces, or by the Radiative Transfer Equation [START_REF] Chandrasekar | Radiative Transfer[END_REF]] (RTE) in volumes. In order to efficiently render images, these equations have been reformulated to model the simulation as a path space sampling problem [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF]], which can be accelerated using dedicated data structures [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]]. The transfer on surfaces is modeled using Bidirectional Reflectance Distribution Functions (BRDF), e.g. modeling the microgeometry of rough surfaces [START_REF] Bitterli | A Position-Free Path Integral for Homogeneous Slabs and Multiple Scattering on Smith Microfacets[END_REF][START_REF] Wang | Position-Free Multiple-Bounce Computations for Smith Microfacet BSDFs[END_REF] or diffractive ones [START_REF] Holzschuch | A Two-Scale Microfacet Reflectance Model Combining Reflection and Diffraction[END_REF], glint paints [START_REF] Xavier Chermain | Procedural Physically-based BRDF for Real-Time Rendering of Glints[END_REF], or volumic micro-structures when the medium is diffuse and/or multi-layer [START_REF] Randrianandrasana | Transfer Matrix Based Layered Materials Rendering[END_REF]]. These models might require a spectral [START_REF] Wilkie | Hero Wavelength Spectral Sampling[END_REF] or wave-optics [START_REF] Steinberg | Towards Practical Physical-Optics Rendering[END_REF] rendering. In case of volumes, macroscopic transport can be modeled in participating media with spatial and spectral heterogeneities [START_REF] Galtier | Integral formulation of nullcollision Monte Carlo algorithms[END_REF][START_REF] Kutz | Spectral and Decomposition Tracking for Rendering Heterogeneous Volumes[END_REF] efficiently stored in dedicated data structures [START_REF] Villefranque | A Path-Tracing Monte Carlo Library for 3-D Radiative Transfer in Highly Resolved Cloudy Atmospheres[END_REF], and can be computed using Monte Carlo, even if being non-linear. Considering that light has a finite speed, movies of light propagation can be generated using transient rendering [START_REF] Jarabo | A Framework for Transient Rendering[END_REF][START_REF] Marco | Progressive Transient Photon Beams[END_REF].

Monte Carlo beyond radiative transfer. In order to model realistic scenes, the light transfer needs to be simulated, but it also requires some phenomena to be described by other physical rules, e.g. fluids dynamics for clouds or smokes [START_REF] Rioux-Lavoie | A Monte Carlo Method for Fluid Simulation[END_REF]]. Hence, even if they focus on radiative transfer, contributions proposed by the computer graphics community are not restricted to this unique physic. For instance, the Walk-on-Sphere (WoS) algorithm, introduced in physics by [START_REF] Muller | Some Continuous Monte Carlo Methods for the Dirichlet Problem[END_REF], has been recently implemented by [START_REF] Sawhney | Monte Carlo Geometry Processing: A Grid-Free Approach to PDE-Based Methods on Volumetric Domains[END_REF] to solve PDE-based geometry processing tasks. This algorithm has been extended to spatially varying coefficients [START_REF] Sawhney | Grid-free Monte Carlo for PDEs with spatially varying coefficients[END_REF]] and bi-directional resolution from the sources [START_REF] Qi | A bidirectional formulation for Walk on Spheres[END_REF]. These works open new perspective to process a wide range of physical models using Monte Carlo. However, in these works, when simulating infrared pictures, the WoS is used to compute temperatures using Monte Carlo, and only then, this information is used to render the final image, in a postprocessing step [START_REF] Sawhney | Grid-free Monte Carlo for PDEs with spatially varying coefficients[END_REF]. By contrast, we focus in the present paper on the idea of coupling several physics, pushing the idea that multiple models can be simulated at once, processing physics and rendering in a single step, using a single path space.

Coupling several physics. Rather than relying on a precomputation, we focus on solving the physical simulation at the same time as rendering. This idea of coupling has already been introduced for light transport in rendering, even if not explicitly stated. Specifically, several works have proposed to model light interaction at different scales for materials with micro-geometries, e.g. rough [START_REF] Bitterli | A Position-Free Path Integral for Homogeneous Slabs and Multiple Scattering on Smith Microfacets[END_REF][START_REF] Heitz | Multiplescattering microfacet BSDFs with the Smith model[END_REF][START_REF] Wang | Position-Free Multiple-Bounce Computations for Smith Microfacet BSDFs[END_REF]], multi-layer [Guillén et al. 2020;[START_REF] Guo | Position-Free Monte Carlo Simulation for Arbitrary Layered BSDFs[END_REF]] of granular diffuse [START_REF] Meng | Multi-Scale Modeling and Rendering of Granular Materials[END_REF] materials. This scale coupling appears when evaluating the BRDF in the scene (macroscopic scale) using a microscopic Monte Carlo estimation. Usually, the goal is not to evaluate this microscopic estimation per se, but rather to build a model in pre-process, e.g., to minimize variance. For these reasons, most models use position-free approximations and ignore lateral transport in the material.

Except this coupling of the RTE at different scales, we are not aware of other work solving multiple physics coupling using Monte Carlo. For solving problems involving different physics, the definition of a single path space seems to be promising: there is no need then to use different Monte Carlo algorithm for going from one physical model to another, but rather to follow a unique path going from one space to another.

Paper organization. Coupled heat simulation (i.e. radiative transfer, convection and conduction, as reminded in Sec. 3) can indeed be expressed in a single path space (Sec 4.1), as demonstrated in theoretical physics by [START_REF] Marc Tregan | Coupling radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm: A general theoretical framework for linear situations[END_REF]. We build upon this result to adapt its formulation to computational and practical constraints (Sec. 4.2). Our proposed coupled model is not only spatial: it is also temporal so that it handles unstationary states (Sec. 4.3). We then show how a single simulation allows to build propagation animations, as well as to change boundary conditions, using path-replay ideas [Mul 2016;Vicini et al. 2021] (Sec. 4.4).

THE PHYSICAL MODEL AND ITS DIDACTIC VERSION

In this section, we review the physical model of heat-transfer, i.e. convection, conduction and radiative transfer. First, we present the model in its general formulation. We recall the RTE with some thermal sources associated to matter in a state of local thermodynamic equilibrium (Sec. 3.1) and the advecto-diffusion model (Sec. 3.2) (i.e. conduction and convection), and we explain their coupling (Sec. 3.3). Since our concern lies essentially on the coupling, we then present a simplified version of the model (Sec. 3.4), which will be used hereafter. A synthetic summary of the fundamentals of thermodynamics is given in App. A.1. Heat-transfer specialists can directly jump to Section 4.

Thermal RTE

The model describing the transport of radiative energy is classically given by a) a linear single-velocity Boltzmann equation in the volume (Eq. 1), and b) the Rendering Equation on the opaque boundaries of the system (Eq. 2):

ì 𝜔 • ì ∇𝐿 𝜆 = -𝜅 𝑎 [𝐿 𝜆 -𝐿 𝑒𝑞 𝜆 (𝜃 )] -𝜅 𝑠 𝐿 𝜆 + ∫ 4𝜋 𝜅 𝑠 𝑝 𝑠 ( ì 𝜔 | ì 𝜔 𝑖 )𝐿 𝜆,𝑖 𝑑 ì 𝜔 𝑖 (1) 𝐿 𝜆 = 𝜖𝐿 𝑒𝑞 𝜆 (𝜃 ) + (1 -𝜖) ∫ 2𝜋 𝑝 𝑟 ( ì 𝜔 | ì 𝜔 𝑖 )𝐿 𝜆,𝑖 𝑑 ì 𝜔 𝑖 (2)
where 𝐿 𝜆 ≡ 𝐿 𝜆 ( ì 𝑥, ì 𝜔, 𝑡) is the radiance (or specific intensity), ì 𝑥, ì 𝜔, 𝑡, 𝜆 denote respectively the position, direction, time and wavelength, with 𝐿 𝜆,𝑖 ≡ 𝐿 𝜆 ( ì 𝑥, -ì 𝜔 𝑖 , 𝑡). The absorption and scattering coefficients (𝜅 𝑎 and 𝜅 𝑠 ), the emissivity 𝜖 as well as the volumic and surface phase functions (𝑝 𝑠 and 𝑝 𝑟 ) are all quantities that potentially depend on ì 𝑥, ì 𝜔, 𝜆. The Planck function 𝐿 𝑒𝑞 𝜆 (𝜃 ) describes the equilibrium radiance, and depends on the temperature 𝜃 and the wavelength.

Note that radiative transfer is considered in its stationary version, due to the time scales of photon propagation compared to the time scales of other transfer modes. Nevertheless, radiance keep time dependency as soon as we consider coupled situations where the temperature 𝜃 evolves over time. Also note that, the variation of the index of refraction with the temperature being negligible in our context, the refractive form of the RTE [START_REF] Ament | Refractive Radiative Transfer Equation[END_REF]] is not considered.

In radiative heat transfer, Equations 1 and 2 are strictly compatible with the equilibrium constraints. This imposes the specific form of the source terms (involving 𝐿 𝑒𝑞 𝜆 ) in the equations, regarding the emission and absorption parts.

Advection-diffusion

The conductive and convective processes that we find respectively in solid media and fluids are typically modeled by an advecto-diffusive

𝜌𝑐 𝜕𝜃 𝜕𝑡 = -ì ∇ • -𝑘 ì ∇𝜃 + 𝜌𝑐 ì 𝑉 𝜃 + 𝑝
which can be written as:

𝜕𝜃 𝜕𝑡 = -ì ∇ • -𝛼 ì ∇𝜃 + ì 𝑉 𝜃 + 𝑆 (3)
where

𝜃 ≡ 𝜃 ( ì 𝑥, 𝑡), ì 𝑉 ≡ ì 𝑉 ( ì 𝑥, 𝑡), 𝛼 ≡ 𝛼 ( ì 𝑥) = 𝑘 ( ì 𝑥) 𝜌𝑐 and 𝑆 ≡ 𝑆 ( ì 𝑥, 𝑡) = 𝑝 ( ì 𝑥,𝑡 )
𝜌𝑐 . The properties 𝑘, 𝜌 and 𝑐 are respectively the thermal conductivity, the density and the heat capacity of the material. The term -𝑘 ì ∇𝜃 represents the surface density vector of conductive flux given by Fourier's law. The term 𝜌𝑐 ì 𝑉 𝜃 is the temperature advection flux (transport of energy by the fluid movement). It vanishes in solids ( ì 𝑉 = ì 0). For fluids, the velocity field ì 𝑉 is prescribed, and how it can be computed by fluid mechanics models is not discussed here (see Perspectives Sec. 6). Finally, in 𝑆, 𝑝 is a source term of voluminal power.

Advecto-diffusive phenomenology appears in numerous fields of physics, well beyond thermal transfers. Classically, numerical grid methods are employed to obtain the solutions. These have the usual limitation in terms of scalability [START_REF] Sawhney | Monte Carlo Geometry Processing: A Grid-Free Approach to PDE-Based Methods on Volumetric Domains[END_REF] and statistical alternatives have emerged in the physics community as early as in the 1960s. In line with the theoretical work of Feynman- Kac [Feynman and Brown 2005;[START_REF] Itô | Diffusion processes and their sample paths: Reprint of the 1974 edition[END_REF][START_REF] Kac | Random walk and the theory of Brownian motion[END_REF][START_REF] Kac | On some connections between probability theory and differential and integral equations[END_REF], who writes the solution of these equations as an expectation on drifted stochastic Brownian processes, several Monte Carlo computational alternatives have emerged. The seminal work on walk-on-sphere proposal [START_REF] Muller | Some Continuous Monte Carlo Methods for the Dirichlet Problem[END_REF]] introduces an elegant algorithm with outstanding convergence properties. Recent work by [START_REF] Sawhney | Walk on Stars: A Grid-Free Monte Carlo Method for PDEs with Neumann Boundary Conditions[END_REF] extends the initial WoS to cases with Neumann-like boundary conditions. Nevertheless, handling heterogeneous situations in diffusivity and velocity field remains an open question and this research area is active, even for stationary cases [Lejay andMaire 2010, 2013;[START_REF] Sawhney | Grid-free Monte Carlo for PDEs with spatially varying coefficients[END_REF]. The instationarity brings an additional difficulty for the management of the events ending inside the sphere. Some bridges have been established between these proposals and the literature of first-pass Green methods [START_REF] Deaconu | The walk on moving spheres: a new tool for simulating Brownian motion's exit time from a domain[END_REF]Deaconu andLejay 2006, 2007;[START_REF] Hwang | Analysis and comparison of Green's function first-passage algorithms with "Walk on Spheres[END_REF]Karl 2019;[START_REF] Karl | Monte Carlo tracking drift-diffusion trajectories algorithm for solving narrow escape problems[END_REF], which open interesting perspectives but, until now, do not allow to scale up to practical geometric and temporal complexity.

Coupling

In most cases, the energy transfers within a system are based on the coupling of the aforementioned modes. The radiative transfer model is naturally coupled to the temperature field through the Planck function, which appears in both the volume and boundaries of the system (Eq. 1 and 2). The temperature field is given by the advecto-diffusive heat transfer equation (Eq. 3) which is affected itself by radiative exchanges. The coupling with radiation occurs in two different ways: through sources in the volume and at the system boundaries. Regarding the first point, the radiative source 𝑆 = 𝑝 𝜌𝑐 in the volume (see Eq. 3) is no longer prescribed but expressed from the radiative flux:

𝑝 = -ì ∇ • ∫ 4𝜋 𝑑 ì 𝜔 ∫ +∞ 0 𝑑𝜆𝐿 𝜆 ì 𝜔 , (4)
where the dependency on the radiance 𝐿 𝜆 becomes explicit. Regarding the second point, exchanges through system boundaries, we use the flux continuity property. The connection condition between a fluid domain and an opaque solid is expressed by the continuity of the temperature (𝜃 is unique at the interface) and the continuity of energy fluxes. The conductive energy flux on the solid side of the interface 𝜑 𝑐𝑜𝑛𝑑,𝑆 is strictly equal to the sum of the radiative and convective fluxes on the fluid side. Since the velocity of the fluid is null at the boundary (no-slip condition), the convective flux is restricted to the conductive part. Hence:

𝜑 𝑐𝑜𝑛𝑑,𝑆 [𝜃 ] = 𝜑 𝑐𝑜𝑛𝑑,𝐹 [𝜃 ] + 𝜑 𝑟𝑎𝑑 [𝐿 𝜆 ], (5)
where [•] denotes the operator dependence on temperature and radiance. Equation 5shows the coupling at the surface between the variables 𝜃 and 𝐿 𝜆 (the temperature field and the radiance at the boundary).

A simplified version of the model to enlighten the coupling

Now, we expose the simplified version of the heat transfer model that will be used in Sec. 4 to design a single path-sampling Monte Carlo algorithm for the coupling of conduction, convection and radiation. For didactic reasons, we will stick to the very basis of coupled heat transfer, making strong simplifications of the initial model presented above. Each simplification made here will be closely reconsidered in Perspectives (Section 6).

Opaque diffuse-gray solids and transparent fluids. The scene is considered to be only made of either opaque solids (no radiation inside the solids, only conduction) or transparent fluids. Radiative transfer reduces therefore to straight line propagation across fluids (neither diffusion nor absorption/emission in the fluid) so it is driven by the standard Rendering Equation (Eq. 2) at all locations along the solid-fluid interfaces.

The radiative properties of solid surfaces are independent of temperature. Their emissivity 𝜖 is independent of direction and wavelength (gray), and reflection can be either diffusive (cosine distributed) or specular (or a mix of both

𝐿 = 𝜖 𝜎𝜃 4 𝑆 𝜋 + (1 -𝜖) ∫ 2𝜋 𝐿 𝑖 𝑝 𝑟 ( ì 𝜔 | ì 𝜔 𝑖 ) 𝑑 ì 𝜔 𝑖 ( 6
)
where 𝜎 is the Boltzmann constant and 𝜃 𝑆 the solid surface temperature. The net radiative flux density at a solid surface of unit normal ì 𝑛 is: If this radiance temperature is unknown, a new radiative sub-path is initiated as when starting from the camera. (c) In the conductive case, after injection within the solid, when the solid temperature is unknown, a walk-on-sphere path is followed until the solid border is found at a location ì 𝑥 𝑤𝑜𝑠 (Eq. 24). The Monte Carlo weight is then 𝜃 𝑆 ( ì 𝑥 𝑤𝑜𝑠 ). (d) When the reached location is at the interface between two solids and when the temperature of this interface is unknown, one of the two solids is randomly selected and the Monte Carlo weight is the temperature of the selected solid 𝛿-away from the interface, i.e. 𝜃 𝑆 ( ì . 26). (e) In the convective case, when 𝜃 𝐹 is unknown, a location ì 𝑥 𝐹 is sampled at the surface of the solid surrounding the fluid cell. The Monte Carlo weight is the solid temperature 𝜃 𝑆 ( ì 𝑥 𝐹 ) at this location (Eq. 28). (f) The algorithm alternates between different sub-path modes until it reaches a known temperature.

𝜑 𝑟𝑎𝑑 = ∫ 2𝜋 𝜖 𝜎𝜃 4 𝑆 𝜋 emitted -𝜖𝐿 𝑖 absorbed ì 𝜔 𝑖 • ì 𝑛 𝑑 ì 𝜔 𝑖 (7)
𝑥 𝛾 -𝛿 1 ì 𝑛) or 𝜃 𝑆 ( ì 𝑥 𝛾 + 𝛿 2 ì 𝑛) (Eq
Uniform solids by parts. Solids are rigid (no deformation velocity field) and are divided into sub-parts where conductivity and specific heat capacity are uniform. For each sub-part, Eq. 3 therefore turns into pure diffusion:

𝜕𝜃 𝑆 𝜕𝑡 = 𝛼 Δ𝜃 𝑆 , (8)
where 𝜃 𝑆 ≡ 𝜃 𝑆 ( ì 𝑥, 𝑡) is the solid temperature. At any location at an interface of unit normal ì 𝑛 between two solids (denoted 𝑠1 and 𝑠2), the conductive flux density on both sides is defined as

𝜑 𝑐𝑜𝑛𝑑 = -𝑘 ì ∇𝜃 𝑆 • ì 𝑛, (9)
while flux continuity implies that

𝜑 𝑐𝑜𝑛𝑑,𝑠1 ( ì 𝑥, 𝑡) = 𝜑 𝑐𝑜𝑛𝑑,𝑠2 ( ì 𝑥, 𝑡). ( 10
)
Perfectly mixed fluid cells. All connected fluid cells are assumed isothermal except for thin boundary layers along the solid-fluid interfaces. Across these layers, the heat flux density 𝜑 𝑐𝑜𝑛𝑣 is modeled as proportional to the difference between the temperature 𝜃 𝐹 of the isothermal fluid core and the surface temperature 𝜃 𝑆 of the facing solid (Newton law):

𝜑 𝑐𝑜𝑛𝑣 = ℎ (𝜃 𝑆 -𝜃 𝐹 ) . ( 11
)
The proportionality coefficient ℎ (the convective exchange coefficient), via the Nusselt number, is a direct translation of the fluid flux inside the boundary layer and is assumed to be independent of the temperature. Under this assumption, as detailed in App. A.1, Equation 3 can be rewritten:

𝜌𝑐V 𝑑𝜃 𝐹 𝑑𝑡 = ∫ 𝒮 ℎ( ì 𝑥) (𝜃 𝑆 -𝜃 𝐹 ) 𝑑 ì 𝑥, (12)
where 𝜃 𝐹 ≡ 𝜃 𝐹 (𝑡), V is the volume of the fluid cell and 𝒮 is its boundary. At any locations along a solid-fluid interface, there is conduction on the solid side and both convection and radiation on the fluid side. Here, energy conservation implies that

𝜑 𝑐𝑜𝑛𝑑 ( ì 𝑥, 𝑡) = 𝜑 𝑐𝑜𝑛𝑣 ( ì 𝑥, 𝑡) + 𝜑 𝑟𝑎𝑑 ( ì 𝑥, 𝑡). ( 13
)
We note that in spite of the simplifications made here, this model is the one that is most commonly used in many application areas.

COUPLING THE THREE HEAT TRANSFER MODES

In this section, we detail how to extend the model described in Sec. 3.4 in order to make it tractable for real-world scenarios. The key physical interactions modeled in our proposal (i.e. transport in medium and interactions at boundaries) define the main steps of our algorithm. The following section is structured around these interactions, the corresponding equations are not only highlighted but also illustrated in Figure 2 in order to simplify the reading. The pseudocode is provided in App. B.

Overview. In Section 4.1, we define how we model backward paths (from camera to sources) propagating through the entire system. These paths combine the three heat transfer modes, and transport temperature information from locations where temperature is imposed (sources) to locations where temperature is sought (camera pixels). Section 4.2 shows how to address the geometric/physical information using only line-surface intersection operators. For this, we replace the standard walk-on-sphere algorithm with a deltasphere random walk algorithm for thermal diffusion inside solids. In Section 4.3 we go beyond the stationary case and extend our model to transient heat-transfer. Finally, we analyze in Section 4.4 the physical pictures associated to the resulting path statistics, and how this information can be used to reconstruct new infrared images very rapidly when changing the sources. The question of further using this propagation information is addressed in Section 5. . If we temporarily assume that the surface temperature of the solid is known, then the Monte Carlo algorithm can be summarized as the following expectation:

Paths combining the three heat transfer modes

𝐿 = E 𝜎𝜃 4 𝑆 ( ì 𝑥 Γ ) 𝜋 . ( 14
)
Defining the radiance temperature as

𝜃 𝑅 = 1 𝜎 𝜋𝐿 1/4 , (15)
and assuming that all temperatures are close to a known reference temperature 𝜃 𝑟𝑒 𝑓 , it is standard practice to linearize radiance as function of radiance temperature using the 1-order development

𝐿 = 𝜎𝜃 4 𝑅 𝜋 ≃ 𝜎𝜃 4 𝑟𝑒 𝑓 𝜋 + 4𝜎𝜃 3 𝑟𝑒 𝑓 𝜋 (𝜃 𝑅 -𝜃 𝑟𝑒 𝑓 ). (16)
In the Monte Carlo context, this linearization can be repeated inside the expectation in Eq. 14:

𝐿 ≃ E       𝜎𝜃 4 𝑟𝑒 𝑓 𝜋 + 4𝜎𝜃 3 𝑟𝑒 𝑓 𝜋 𝜃 𝑆 ( ì 𝑥 Γ ) -𝜃 𝑟𝑒 𝑓       (17) = 𝜎𝜃 4 𝑟𝑒 𝑓 𝜋 + 4𝜎𝜃 3 𝑟𝑒 𝑓 𝜋 E [𝜃 𝑆 ( ì 𝑥 Γ )] -𝜃 𝑟𝑒 𝑓 (18)
and the right sides of Eq. 16 and 18 can be equaled to give the radiance temperature expressed itself as an expectation:

radiance temperature

𝜃 𝑅 = E [𝜃 𝑆 ( ì 𝑥 Γ )] . (19)
As detailed in supplementary material (Sec. 3), the linearization approximation can be removed, in which case Equation 19expresses 𝜃 𝑅 as the expectation of a random variable whose definition involves 𝜃 𝑅 . We denote this situation as non-linear coupling, further discussed in Section 6. With these equations in hand, we have a Monte Carlo algorithm where the radiance temperature is modeled as an expectation of solid-surface temperatures. But in a coupled heat transfer problem, there are surfaces where the temperature is not known. Instead, it results from the combined presence of conduction on the solid side, and convection plus radiation on the fluid side. Equation 13translates this combination ; combined with the previous definitions and the same linearization of (𝜎𝜃 4 𝑆 -𝜎𝜃 4 𝑅 ) as for Eq. 16, it leads to

-𝑘 ì ∇𝜃 𝑆 • ì 𝑛 = ℎ (𝜃 𝑆 -𝜃 𝐹 ) + ∫ 2𝜋 ℎ 𝑅 (𝜃 𝑆 -𝜃 𝑅 ) ì 𝜔 • ì 𝑛 𝜋 𝑑 ì 𝜔 ( 20
)
where ℎ 𝑅 = 4𝜎𝜖𝜃 3 𝑟𝑒 𝑓 . The term ì ∇𝜃 𝑆 • ì 𝑛 can be approximated using finite difference:

-ì ∇𝜃 𝑆 • ì 𝑛 ≈ 𝜃 𝑆 ( ì 𝑥 -𝛿 ì 𝑛) -𝜃 𝑆 ( ì 𝑥) 𝛿 , (21)
where 𝛿 must be small compared to the solid thickness. Thanks to this formulation, we can express 𝜃 𝑆 as an expectation:

𝜃 𝑆 ( ì 𝑥) = 𝑃 𝑐𝑜𝑛𝑑 𝜃 𝑆 ( ì 𝑥 -𝛿 ì 𝑛) + 𝑃 𝑐𝑜𝑛𝑣 𝜃 𝐹 + 𝑃 𝑟𝑎𝑑 ∫ 2𝜋 𝜃 𝑅 ì 𝜔 • ì 𝑛 𝜋 𝑑 ì 𝜔, (22)
with

𝑃 𝑐𝑜𝑛𝑣 = ℎ 𝑘 𝛿 +ℎ+ℎ 𝑅 , 𝑃 𝑟𝑎𝑑 = ℎ 𝑅 𝑘 𝛿 +ℎ+ℎ 𝑅 and 𝑃 𝑐𝑜𝑛𝑑 = 1-(𝑃 𝑐𝑜𝑛𝑣 +𝑃 𝑟𝑎𝑑 ).
Note that, contrary to the others simplifications and approximations we made previously, approximation ( 21) is required to express 𝜃 𝑆 as an expectation, and can not be released.

To represent the test with three branches, we define the random vector (H 𝑐𝑜𝑛𝑑 , H 𝑐𝑜𝑛𝑣 , H 𝑟𝑎𝑑 ) on support {(1, 0, 0), (0, 1, 0), (0, 0, 1)} with respective probabilities 𝑃 𝑐𝑜𝑛𝑑 , 𝑃 𝑐𝑜𝑛𝑣 and 𝑃 𝑟𝑎𝑑 , and defining ì Ω as a random direction cosine-distributed around ì 𝑛, then solid temperature at a solid-fluid interface

𝜃 𝑆 ( ì 𝑥) = E H 𝑐𝑜𝑛𝑑 𝜃 𝑆 ( ì 𝑥 -𝛿 ì 𝑛) + H 𝑐𝑜𝑛𝑣 𝜃 𝐹 + H 𝑟𝑎𝑑 𝜃 𝑅 ( ì 𝑥, ì Ω) . ( 23
)
Reporting Equation 19into Equation 23we can use double randomization [Maire and Nguyen 2016]: expectation of an expectation is an expectation and therefore linearly combining two path-sampling Monte Carlo algorithm leads to one single path-sampling Monte Carlo.

So we can now estimate the radiance temperature at an infrared camera location by sampling a multiple reflection radiative path until it stops at a solid-fluid interface. Then, in case the solid surface temperature is not known, the algorithm simply continue the same path by

• randomly choosing the heat transfer mode according to 𝑃 𝑐𝑜𝑛𝑑 , 𝑃 𝑐𝑜𝑛𝑣 and 𝑃 𝑟𝑎𝑑 ; • in case conduction is selected, retaining the solid temperature 𝜃 𝑆 ( ì 𝑥 -𝛿 ì 𝑛) at a distance 𝛿 away from the interface; • in case convection is selected, retaining the fluid temperature 𝜃 𝐹 ; • in case radiation is selected, first sampling the direction ì 𝜔 and retaining the radiance temperature 𝜃 𝑅 ( ì 𝑥, ì 𝜔). The same logic of double randomization is used to continue the algorithm for the selected heat transfer mode when the temperature to be retained is unknown.

The conduction term in Equation 23 requires computing 𝜃 𝑆 ( ì 𝑥 -𝛿 ì 𝑛). This can be translated as the sampling of a walk-on-sphere random path 𝑊 𝑂𝑆 within the corresponding uniform solid sub-part, the path ending at a location ì 𝑥 𝑊 𝑂𝑆 along the sub-part boundary:

walk-on-sphere

𝜃 𝑆 ( ì 𝑥 -𝛿 ì 𝑛) = E [𝜃 𝑆 ( ì 𝑥 𝑊 𝑂𝑆 )] . (24)
If at ì 𝑥 𝑊 𝑂𝑆 the solid temperature is known, the algorithm stops retaining this temperature as the Monte Carlo weight. If the solid temperature is not known there are two cases. Either the sub-part boundary is a solid-fluid interface, in which case we come back to Equation 23 and continue the algorithm by selecting again one heat transfer mode. Or the sub-part boundary is a solid-solid interface, in which case we need to make a finite difference approximation on , Vol. 42, No. 4, Article . Publication date: August 2023. both sides. Equation 10becomes

𝑘 1 𝜃 𝑆 ( ì 𝑥 -𝛿 1 ì 𝑛) -𝜃 𝑆 ( ì 𝑥) 𝛿 1 ≈ 𝑘 2 𝜃 𝑆 ( ì 𝑥) -𝜃 𝑆 ( ì 𝑥 + 𝛿 2 ì 𝑛) 𝛿 2 , (25)
and solid temperature at a solid-solid interface

𝜃 𝑆 ( ì 𝑥) = E H 𝑐𝑜𝑛𝑑,1 𝜃 𝑆 ( ì 𝑥 -𝛿 1 ì 𝑛) + H 𝑐𝑜𝑛𝑑,2 𝜃 𝑆 ( ì 𝑥 + 𝛿 2 ì 𝑛) , (26)
where the Bernoulli random variables H 𝑐𝑜𝑛𝑑,1 has probability

𝑃 𝑐𝑜𝑛𝑑,1 = 𝑘 1 𝛿 1 /( 𝑘 1 𝛿 1 + 𝑘 2 𝛿 2
) and H 𝑐𝑜𝑛𝑑,2 = 1 -H 𝑐𝑜𝑛𝑑,1 . So the algorithm continues with either a 𝛿 1 jump inside solid 1 or a 𝛿 2 jump inside solid 2, and starts a new walk-on-sphere algorithm from the reached location.

The convection term in Equation 23requires computing 𝜃 𝐹 . The stationary version of Equation 12(i.e. 𝑑𝜃 𝐹 𝑑𝑡 = 0) leads to

𝜃 𝐹 = ∫ 𝒮 ℎ( ì 𝑥)𝜃 𝑆 ( ì 𝑥)𝑑 ì 𝑥 ∫ 𝒮 ℎ( ì 𝑥)𝑑 ì 𝑥 , (27)
and therefore fluid temperature

𝜃 𝐹 = E 𝜃 𝑆 ( ì 𝑋 𝒮 ) , (28)
where ì 𝑋 𝒮 is a random location along the fluid cell boundary with probability density proportional to the convective exchange coefficient. So when 𝜃 𝐹 is needed, a location is sampled on the boundary, with location probability as a function of the intensity of convective exchanges, and retains the solid temperature at this location.

Finally, the radiation term in Equation 23 requires computing 𝜃 𝑅 ( ì 𝑥, ì 𝜔), and this is identical to our initial question of estimating 𝜃 𝑅 at a camera location.

Altogether, as summarized in Figure 2, Equations 19, 23, 26 and 28, combined with a walk-on-sphere algorithm for pure conduction within the uniform solid sub-parts, provide a closed Monte Carlo algorithm. This algorithm samples paths across the system, switches between heat transfer modes at the solid-fluid interfaces, ends when reaching a location where the temperature is known and retains this temperature value as the Monte Carlo weight. At this stage, although the core of the algorithm is defined, its efficient implementation impacts the choices made on the physical model itself in order to benefit from complex geometry scalability.

Walk-on-delta-sphere

According to our hypothesis, the solids are divided into uniform sub-parts and the fluids are transparent, which means that all the parts properties can be accessed from their envelop. Same thing for the heat transfer parameters: conductivity and diffusivity on both sides for solid-solid interfaces; conductivity and diffusivity on the solid side, emissivity and convective heat transfer coefficient on the fluid side for solid-fluid interfaces. Fast access to these data boils down to the fast sampling of surface locations for each of the three heat transfer modes.

As far as radiation is concerned, standard ray-tracing strategies are immediately available to optimize the computational effort associated to line-surface intersections [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]. For convection, sampling ì 𝑋 𝒮 is trivial: there is no dependence on location and cumulated distributions can be pre-constructed for the geometrical primitives (e.g. triangles) of the envelope, weighting them by their size and the local value of the convective exchange coefficient. For conduction, walk-on-sphere algorithms open a new set of questions. If the discussion was strictly restricted to the basic heat transfer model of the present section, then we could simply use the acceleration strategies that have recently been proposed for finding the triangle corresponding to the smallest tangent sphere as function of the sphere-center location [START_REF] Shellshear | Fast distance queries for triangles, lines, and points using SSE instructions[END_REF][START_REF] Wald | Embree: a kernel framework for efficient CPU ray tracing[END_REF]]. But our question is broader: we want to set up pathsampling algorithms flexible enough for extension to heterogeneous diffusivity, heterogeneous velocity (when advection is considered), volume-coupling with radiation within semi-transparent media (fluids or solids), and nonlinear coupling with additional physical models (e.g., diffusivity and velocity depend on temperature). These questions are at least partially addressed in the walk-on-sphere literature [Sabelfeld 2019;[START_REF] Karl | A global random walk on grid algorithm for second order elliptic equations[END_REF], but they remain very open. We then take a practical decision: considering that we already made the choice of using a finite difference approximation for the handling of solid-solid and solid-fluid interfaces, we choose to replace the walk-on-sphere algorithm with a walk-on-𝛿-sphere algorithm [START_REF] Marc Tregan | Coupling radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm: A general theoretical framework for linear situations[END_REF], i.e. sampling the next location on a sphere of fixed radius 𝛿 instead of the tangent sphere. For stationary diffusion inside uniform solids, exactly like walk-on-sphere, the algorithm is exact far from the boundary. Starting from a location ì 𝑥, we handle the presence of the boundary with the following procedure:

• We sample an isotropically distributed random direction ì Ω.

• We test if a triangle is intersected before the distance 𝛿 in direction ì Ω. If an intersection is found we store the intersection distance 𝑠 + , otherwise 𝑠 + = 𝛿.

• We do the same, defining 𝑠 -for the opposite directionì Ω. • We define δ = min(𝑠 -, 𝑠 + ), and we jump at ì 𝑥 + δ ì Ω, i.e. on the sphere of radius δ (which will be of radius 𝛿 if no triangle was intersected). This is a first order spatial approximation (i.e. exact when the temperature profile is trilinear). The pseudo-code of this algorithm, as well as the full Monte Carlo procedure are given in App. B.

Extension to transient heat transfer

We model transient conduction and convection by keeping the temporal terms in Equations 8 and 12, respectively. During the path sampling, we also track time, such that the backward construction of the path may stop because the initial time 𝑡 𝐼 is reached. As the initial conditions of the system are known, all temperatures are known at 𝑡 𝐼 , whatever the reached location of the Monte Carlo weight. Convection. Because of the kinetic nature of Eq. 12, time tracking can be directly computed for convection. The formal solution writes

𝜃 𝐹 (𝑡) = exp - 𝑡 -𝑡 𝐼 𝜏 𝑐𝑜𝑛𝑣 𝜃 𝐹 (𝑡 𝐼 )+ ∫ 𝑡 -𝑡 𝐼 0 1 𝜏 𝑐𝑜𝑛𝑣 exp - 𝑡 𝑏 𝜏 𝑐𝑜𝑛𝑣 𝜃 * 𝐹 (𝑡 -𝑡 𝑏 ) 𝑑𝑡 𝑏 (29)
with

𝜏 𝑐𝑜𝑛𝑣 = 𝜌𝑐 V ℎ𝒮 and ℎ = 1 𝒮 ∫ 𝒮 ℎ( ì 𝑥)𝑑 ì 𝑥,
and where 𝜃 * 𝐹 (𝑡) is the same average surface temperature as in the stationary case. It was already formulated as an expectation in Eq. 28, leading to:

𝜃 * 𝐹 (𝑡) = 1 ℎ𝒮 ∫ 𝒮 ℎ( ì 𝑥)𝜃 𝑆 ( ì 𝑥, 𝑡)𝑑 ì 𝑥 = E 𝜃 𝑆 ( ì 𝑋 𝒮 , 𝑡) . ( 30
)
Observing that exp

-𝑡 -𝑡 𝐼 𝜏 𝑐𝑜𝑛𝑣 = ∫ +∞ 𝑡 -𝑡 𝐼 1 𝜏 𝑐𝑜𝑛𝑣 exp -𝑡 𝑏
𝜏 𝑐𝑜𝑛𝑣 𝑑𝑡 𝑏 , Equation 29 can be written as an expectation by defining an exponentially distributed random variable 𝑇 𝑏,𝑐𝑜𝑛𝑣 of mean 𝜏 𝑐𝑜𝑛𝑣 , itself defining H 𝐼,𝑐𝑜𝑛𝑣 = 𝐻 (𝑇 𝑏,𝑐𝑜𝑛𝑣 -𝑡 + 𝑡 𝐼 ), where 𝐻 is the Heaviside function. H 𝐼,𝑐𝑜𝑛𝑣 is 1 if the initial condition is retained, 0 otherwise:

𝜃 𝐹 (𝑡) = ∫ +∞ 0 1 𝜏 𝑐𝑜𝑛𝑣 exp - 𝑡 𝑏 𝜏 𝑐𝑜𝑛𝑣 𝐻 (𝑡 𝑏 -𝑡 + 𝑡 𝐼 )𝜃 𝐹 (𝑡 𝐼 ) + 𝐻 (-𝑡 𝑏 + 𝑡 -𝑡 𝐼 )𝜃 * 𝐹 ( ì 𝑥, 𝑡 -𝑡 𝑏 ) 𝑑𝑡 𝑏 = E H 𝐼,𝑐𝑜𝑛𝑣 𝜃 𝐹 (𝑡 𝐼 ) + (1 -H 𝐼,𝑐𝑜𝑛𝑣 )E 𝜃 𝑆 ( ì 𝑋 𝒮 , 𝑡 -𝑇 𝑏,𝑐𝑜𝑛𝑣 ) = E H 𝐼,𝑐𝑜𝑛𝑣 𝜃 𝐹 (𝑡 𝐼 ) + (1 -H 𝐼,𝑐𝑜𝑛𝑣 )𝜃 𝑆 ( ì 𝑋 𝒮 , 𝑡 -𝑇 𝑏,𝑐𝑜𝑛𝑣 ) (31)
The resulting transient convective sub-path is sampled the following way:

• 𝑇 𝑏,𝑐𝑜𝑛𝑣 is sampled and is interpreted as an exponentially distributed backward jump in time; • when this jump leads to a time 𝑡 -𝑇 𝑏,𝑐𝑜𝑛𝑣 before 𝑡 𝐼 , the initial fluid temperature is retained; • otherwise ì 𝑋 𝒮 is sampled to get a position at the surface of the solid surrounding the fluid (as for the stationary path) and the solid temperature is retained, at this location at time 𝑡 -𝑇 𝑏,𝑐𝑜𝑛𝑣 .

Conduction. It is handled exactly the same way as convection, as detailed below. Up to now we used a first order approximation in space for diffusion when estimating the temperature of a solid-solid or a solid-fluid interface (re-injection inside the solid at distance 𝛿 when conduction is selected) and when adapting the walk-on-𝛿spheres algorithm in the vicinity of the boundary. We use the same level of spatial approximation to write the following approximate version of transient diffusion equation (see Eq. 8). As detailed in App. A.2, this is done by using a finite difference to transform the Laplace operator in an arbitrarily chosen Cartesian coordinate system and then averaging over all the possible orientations of the coordinate system:

𝜕𝜃 𝑆 𝜕𝑡 = 𝛼 - 6 δ2 𝜃 𝑆 + 6 δ2 𝜃 * 𝑆 , (32)
where 𝜃 * 𝑆 ( ì 𝑥, 𝑡) is the same average temperature over the 𝛿-sphere as for the stationary case. It can be translated as an expectation using the isotropic random direction ì Ω:

𝜃 * 𝑆 ( ì 𝑥, 𝑡) = E 𝜃 𝑆 ( ì 𝑥 + δ ì Ω, 𝑡) . (33)
The formal solution of Equation 32 is

𝜃 𝑆 ( ì 𝑥, 𝑡) = exp - 𝑡 -𝑡 𝐼 𝜏 𝑐𝑜𝑛𝑑 𝜃 𝑆 ( ì 𝑥, 𝑡 𝐼 )+ ∫ 𝑡 -𝑡 𝐼 0 1 𝜏 𝑐𝑜𝑛𝑑 exp - 𝑡 𝑏 𝜏 𝑐𝑜𝑛𝑑 𝜃 * 𝑆 ( ì 𝑥, 𝑡 -𝑡 𝑏 )𝑑𝑡 𝑏 (34)
which is consistent with the solution of Eq. 8 as detailed in supplementary material (Sec. 2). As for Eq. 29, we directly write

𝜃 𝑆 ( ì 𝑥, 𝑡) = E H 𝐼,𝑐𝑜𝑛𝑑 𝜃 𝑆 ( ì 𝑥, 𝑡 𝐼 ) + (1 -H 𝐼,𝑐𝑜𝑛𝑑 )𝜃 𝑆 ( ì 𝑥 + δ ì Ω, 𝑡 -𝑇 𝑏,𝑐𝑜𝑛𝑑 ) , (35)
and the corresponding sampling algorithm is exactly the same as for transient convection.

Storing and displaying path information

For each camera pixel, the above described algorithm is used to sample 𝑛 paths starting from the camera at the observation time. The path traverses the system backward in time, switching from one heat transfer mode to the other, until either the initial time is reached (initial condition) or it crosses a location where the temperature is known (imposed temperature). These are the two only ways for a path to find its end, which means that the Monte Carlo weight is always the value of the temperature at the end of the path. As for standard rendering, generating the image does not require to store the complete history of the path in addition to the Monte Carlo weights. But in the context of thermal simulation, there are several practical use-cases where storing the entire path space can have a very strong impact. For instance, in the context of energy transition, engineers and citizens can use infrared cameras to understand the influences of a fluctuating environment, get insight about their heat control choices, or about the insulation properties of building. From this perspective, statistically analyzing the heat transfer paths gives more information than a simple average of the associated weight (i.e. the radiative flux at the pixel surface). Displaying the paths within any scientific viewer has already a strong practical significance by itself. The same statement is made to draw an analysis from the visualization of light paths [Kopylov and Dmitriev 2000], we stick here to simple visualizations although more advanced methods exist (e.g. [START_REF] Zirr | Visualization of Coherent Structures of Light Transport[END_REF]).

We illustrate this point in Section 5, together with two other useful practices that emerge as straightforward consequences of the present algorithmic proposal due to the linearity of the physical model. First, storing only the end of the path allows to very rapidly replay the rendering for new initial conditions or new imposed temperatures (see Figure 3). This idea is similar to relighting [Mul 2016] where sources can be modified after the simulation. Note that the term path-replay comes from the work of [START_REF] Vicini | Path Replay Backpropagation: Differentiating Light Paths using Constant Memory and Linear Time[END_REF]. Second, storing all surface encounters allows to change the observation time. Indeed, when replaying the rendering at any date between the initial condition and the observation time for which the paths were constructed, the initial condition will be found somewhere between two known surface encounters. Specifically, it will be located inside an identified uniform solid sub-part or an identified fluid cell. Using this information, the Monte Carlo weight can be directly recomputed with a known initial temperature value. In practice, this means that a single Monte Carlo run computes all the required information to generate a sequence of infrared images showing the temporal response evolving over time, for a fixed camera location and a fixed geometry (see Figure 1 and supplementary video).

RESULTS

In this section, we show the usefulness of coupling several physics and manipulating the propagator on some concrete applications, mainly related to building construction and rehabilitation. The theoretical tools enable to analyze infrared renderings in the thermal context. Here, we show some results on infrared images, but it is much more general: a flux, an electric consumption, etc. could also be computed the same way.

Implementation details. We developed our prototype implementation using Stardis [Sta 2018] for Monte Carlo simulations and evaluations, restricting ourselves to the didactic model (i.e. having uniform conditions). This assumption is further discussed in supplementary material (Sec. 4). We implemented the path-replay of an infrared image and the time propagator, the later being restricted to a uniform initial temperature. The surfacic phase function 𝑝 𝑟 is composed of a specular part and a Lambertian part, sampled accordingly to the specular fraction parameter 𝑠 (the probability that reflection is specular otherwise diffuse). We set the reinjection distance 𝛿 about 1/20 of the object's characteristic size. For all the presented infrared images, we provide in table 1 either their simulation or evaluation times, performed on an Intel Xeon Silver 4208 machine with 32 CPUs. As our implementation is not optimized, these timings might Table 1. Scene configurations and associated computation times, for either full Monte Carlo simulation (s) or replay (r). For teaser, replay times are only two orders of magnitude smaller than MC times. This is because the replay code anticipates for further internal sources (see supplemental). For farm simulations, the replay code has been simplified for engineering needs and the replay to Monte Carlo ratio is of the order of 5 • 10 -7 (Fig. 6 replay times estimated on average over a large number of run replications). be further reduced, (e.g., using vectorization). In any case, the purpose of our prototype implementation is to demonstrate that the didactic model already enables interesting results and analysis of the coupled heat transfer.

Validation. As there is no other code against which to validate the infrared images, the Stardis code is validated on temperature field computation and all the references are available at the webpage https://www.meso-star.com/projects/stardis/validation.html. These validations are made against analytical models on simple geometries and cross-comparisons with commercial finite difference codes (ANSYS Fluent [Flu 2017] and Comsol [Com 2020]) on more complex geometries. Note that the open access Syrthes [Syr 2021] code embeds Stardis, which allows systematic comparisons on the exact same CAD material. These validations are provided for cases that remain tractable by deterministic approaches (requiring the discretization of the whole volume) but cannot be made for complex geometries such as cities where solving on mesh is impractical.

The propagator implementation is validated against the classical Monte Carlo simulation, both on a temperature computation (see Fig. 7) and infrared renderings (see Fig. 1 -right).

Infrared imaging using coupled-paths : the case of the Heat exchanger

We first illustrate how we use coupled-paths to generate an infrared image, here of a heat exchanger (in Fig. 1) made of a Kelvin cells foam

(𝑘 = 237 W m -1 K -1 , 𝜌 = 2700kg m -3 , 𝑐 = 890 J K -1 kg -1
) with emissivity 𝜖 𝑒 = 0.9 and specular fraction 𝑠 𝑒 = 0.5 for the reflective part, placed on a floor with emissivity 𝜖 𝑓 = 0.8 and a specular fraction 𝑠 𝑓 = 0.8. The whole is plunged within a uniform radiative background maintained at 680K. For the sake of clarity, we only consider coupled conductive and radiative paths in this illustration. At initial time, the temperature of the system is set at 𝜃 𝐼 = 700K everywhere, the left and right walls of the exchanger are set at 850K and 750K respectively, and keep being imposed. The temperature of the floor remains imposed at 𝜃 𝐼 = 700K. Left panel in Fig. 1 shows the infrared rendering at infinite time (𝑡 𝑜𝑏𝑠 = +∞), i.e. the stationary state: imposed temperatures at the walls propagate towards the core of the exchanger, as expected, but not very deeply (for the given set of parameters) because the central part is conversely cooled by radiative exchange with background and floor. Middle panel illustrates how this rendering has been built using coupling between conductive and radiative paths. Starting from pixels of the virtual camera, some radiative paths propagate back to the structure up to an intersection point. There, they can switch to conduction within the material or be reflected. Then, they chain conductive paths through the material (blue subpaths) and radiative paths between internal faces (red subpaths), until they either meet a left/right wall or the floor where the temperature is known (e.g. blue dot on the left wall), or exit to radiative background. Using such paths, each pixel value can be estimated in a single procedure integrating the fully coupled thermal transfer, without having to compute the temperature field prior the rendering.

Visualizing paths helps to interpret the results. For instance, the probability that a path ends at the known temperature of a wall is higher for paths that first the structure near either wall. Hence such paths will convey more often the wall information (850K or 750K) than the background radiative value (680K) or floor (700K). Conversely, paths that meet the structure near the core will more often end by radiative exits to background or floor, and meet a wall only rarely so that the exchanger core remains darker. This contrast between core and walls will depend on the probability exiting to background or floor before meeting either wall, a probability which is driven by the material properties (e.g. the higher the conductivity, the lower this probability) and the foam density (the higher the density, the lower this probability). We also note that the floor having a high specular part, paths that hit it as first intersection from the camera can bounce back to the structure and then retrieve its reflected image, albeit attenuated because the floor emissivity is high so that its own imposed temperature prevails : the floor here acts as an attenuated infrared mirror.

As mentioned in Section 4, storing information about the paths that have been realized in this first rendering can be reused to hugely accelerate rendering when other conditions apply (see the supplementary video). We illustrate this point with the right panel of Fig. 1, where this information is used to render transient states at 𝑡 𝑜𝑏𝑠 = 30s and 𝑡 𝑜𝑏𝑠 = 150s. Contrasting with stationary state, paths can now meet the initial condition back in time (in addition to ending at walls, floor or background), so that some of the paths that would have ended at the wall at stationary state now meet the initial condition before reaching the wall. As a consequence, they retrieve 𝜃 𝐼 = 700K instead of the higher wall temperature, resulting overall into a lower temperature. Similarly, some of the paths exiting to background from within the core at stationary state (retrieving then 680K) now meet the initial condition before exiting, and retrieve 700K, resulting overall into a higher temperature. Of course, the probability of meeting the initial condition before intersecting a wall / floor or exiting to background vanishes as time grows and the full history to the stationary state can be rebuilt interactively (see the supplementary video) with no need to sample again the full path space for each snapshot.

Scalability

In Figure 4 we reproduce an image that was displayed in [START_REF] Villefranque | The "teapot in a city": A paradigm shift in urban climate modeling[END_REF] to illustrate the ability to synthesize urban infrared images. It includes heat transfer within the buildings at each scale, as well as transfer between buildings (via the ground or via radiative exchanges), a feature typically required for the modeling of urban heat islands. In our scene, the environment is cold and each room within each building is heated with an air conditioning system. As a typical output, this image highlights the thermal bridges associated to concrete slabs and balconies, as well as a strong heat loss due to an insufficient thermal insulation of the roof.

Figure 4 shows that our objective of a full geometric scalability is reached: computation times are independent of the number of building in the scene. With the choices we made in Section 4, data-access acceleration is exactly the same as for todays standard rendering and we can indeed expect that computation times are insensitive to the number of geometric primitives. But in this scene, interactions are even more complicated, as the buildings are interacting with each other. Despite this, the computation times are the same when imaging 1, 9, 81, or 961 buildings. Under the hood, it means that the average computation time required to construct the paths is nearly the same if the building is surrounded by a known environment (the path stops when it exits the building), or if path can reach the neighbor buildings and continue their way inside these buildings. Note that inside a building, a path can have a wide range of interactions, e.g., cross the walls by conduction, interact with internal air by convection inside the rooms, or jump from one wall location to another by radiation. We believe that the stable simulation times are due to the fact that cities are dominantly extended in two dimensions. Making the scene more complex does not change significantly the average time a heat transfer path travels through the housing before it exits by convection toward the external surrounding air or by radiation toward the upper atmosphere, both here assumed at known temperatures. Computation times would on the contrary increase when the geometry becomes more complex if this increases the length that paths must travel before exit, i.e. before encountering a known temperature or an initial condition (e.g. an extended three-dimensional structure at stationary regime with known temperatures only at the outside). So computation times may change when changing the geometry because of heat transfer features changing the structure of the paths (as it is the case for multiple reflection and multiple scattering in standard rendering), but not because more triangles were used to accurately describe the scene.

Application: a thermal insulation study

We now use a case study of thermal insulation to illustrate how the Monte Carlo resolution of the coupled heat-transfer model may be used. We provide a deeper description in supplementary material (Sec. 1). Analysis of an infrared image. How do we read an infrared image to thermally understand a building and optimize its design (esp. to minimize energy loss)? We take as an example a simplified building inspired from a typical farm from central France made of a room for living topped by an empty hayloft, and we question the difference of its thermal behavior when insulation is put inside versus outside. A wall-mounted heater is turned on in the room and we take, at the stationary state, a nightly infrared picture of the whole building from outside (Fig. 5, upper panel). At first sight, we obviously diagnose the thermal bridge at the ground level in both cases (which calls for primary action to prevent this). To go further, interpreting the shallow temperature difference in external surfaces temperature at hayloft and room levels is less intuitive. To this end, we can take a closer look at the paths themselves, starting here from a probe point located on the inside of the room (Fig. 5, lower panels. Note that the 3D paths are projected on the plane).

To illustrate how reading paths can help, we first focus on three of them, each one starting in convective mode from the probe point (gray paths starting from the gray sphere in Fig. 5, middle panel). With internal insulation (left), we see one of them ending at the outside soil temperature, a second one ending at the heater, and the third one ending at the external air temperature near the roof. Changing the configuration to external insulation (right), these three paths start the same but undergo very different histories: one reaches external air by conduction through the gable, and the two other ones end at the heater by convection. Most importantly, we visually note that paths now chain conduction inside the walls more often. Confirming these "statistics with the eyes" with 60 paths (Fig. 5, lower panel), we understand that a larger fraction of paths can reach the hayloft, thereby transporting heat upside, such heat being then lost for the living room. Reading paths then shows that the inside walls are forefront: when unconstrained under external insulation, they act as thermal bridges towards the empty hayloft and homogenize temperature over the whole building.

Varying conditions. As a consequence of the computational cost, storing as much information as we can about the paths statistics can be of a huge interest since they are unaffected by time or boundaries conditions (such as different heating commands, alternative weather conditions, or even for integrating day/night alternation over days, ...). Computing the propagator only once for each insulation choice and camera settings is information enough to explore sensitivity to boundary conditions. In supplementary material (Sec. 1), we show a typical engineering study which compares the responses of the internally/externally-insulated buildings when changing the altitude and the weather condition from a foggy atmosphere to a clear sky. After one Monte Carlo run has been made, predicting temperature under either of these conditions is instantaneous. Such predictions can be made using only four propagator values (four scalars to be multiplied by their associated boundary temperature values). They represent the four imposed temperatures in our scene: the yearly average temperature at 3 meters deep into the ground (related to the altitude), the heater temperature, the outside air temperature (related to the altitude) and the radiance temperature in directions coming from the upper atmosphere (related to weather condition).

Transient states of thermal comfort. In the same spirit as for the heat exchanger, the propagator can also be used to replay the whole dynamics of heating the room. As an illustration, we computed the two propagators (internal / external insulations) for the temperature of the room at stationary state. From these, we can estimate this temperature any sooner in time (Fig. 6). Here again, with only one single Monte Carlo run for each insulation mode, we capture the information required for any transient dynamic. In our example, with the same heating temperature, the external insulation leads to a lower stationary temperature of about 2 K, but more strikingly, the time needed to reach it can be up to one month, while the internal insulation brings thermal comfort within one hour after the heating is turned on. The retarding effect of exterior insulation depends of course on the properties of the insulation materials and the walls and can be adjusted, for instance to phase the retarded effect of solar sources with the house occupation times. The walls of the considered farm are made of very thick and dense rocks that explain the extremely long duration of transient regimes, but modern buildings have walls that are thermally less capacitive.

Computation time, variance and memory

We have shown above the immediate applicability of our simple physical model of heat-transfer. Here, we give an overview of the performance characteristics of our method, typical of Monte Carlo approaches. In Fig. 8, the computation time and variance are provided as functions of the number of samples, when computing the temperature of the lower room of the farm. As mentioned in Section 5.2, the computation time depends on the underlying physical phenomena, especially the ability of the paths to quickly reach the source. In our case, conductive paths are longer to compute than radiative and convective ones, because of the number of walk steps required to reach the boundary. Also, adverse probabilities at interfaces (in Eq. 23) can trap trajectories within a geometry and negatively impact the computation time as described in the next Section 5.5.

As for classical rendering, point sources are difficult to reach, resulting in a higher variance. Indeed, we compare in Figure 8 the variance estimates when the probe is put inside vs outside the farm. The variance for the outside probe is larger because less than 2.5% of the paths actually reach the heater due to the wall insulation which prevents access to the hot source (more details in the supplementary material (Sec. 1).

As the construction of the paths is time-consuming, their replay with different boundary or initial conditions as well as observation time is of high interest. This storage is linear wrt the image resolution and the number of samples per pixel to be generic: for spatio-temporal replays and error bar recomputation. A path-replay Fig. 9. Path entrapment due to high conductivity ratios. Sampling heat transfer paths in electronics is a well representative example of how computation times may be dominated by entrapment. Electric connections conduct heat efficiently and tend to be quasi-isothermal. In path sampling terms, this is translated by very long paths inside the connectors because each time the path reaches the boundary of a connector the probability is high that it is reinjected inside the connector. Here the connectors are spheres with only a very small contact area with other solids of similar conductivity levels. In the displayed example, one single conductive path (in blue, from entry point in red to exit point in white) very longly visits the sphere before it finds an exit. Mesh courtesy of DeepConcept.

to change any uniform boundary condition in a scene only requires to store the source identifier for each sample. When these conditions are non uniform, the spatio-temporal impact of the path on the source is needed. To replay the simulation for a previous observation time, we need to additionally store the lifetime of each path if we assume a homogeneous initial temperature. Our prototype currently stores 116 bytes per path, reusing a Stardis existing structure which was originally neither designed to replay images nor observation times. This implementation, which could be optimized, is already practical. However, considering heterogeneous initial temperatures would drastically change the storage cost because time passage between medium would have to be saved.

A limit case with path-entrapment

We now detail a practical limitation we have identified, namely a case with path-entrapment. In the example presented in Fig. 9, a highly conductive solid which is embedded inside another less conductive solid, when a conductive path reaches the highly conductive solid, it tends to stay inside it. When a path reaches the interface, the probability that it is re-injected inside the very conductive part is high. So it revisits this part over and over until it finally exits. But at this exit, it is injected in the other solid, still in the close vicinity of the interface and has therefore a strong probability to be trapped again, etc, see Figure 9. Of course, this is indeed physically meaningful: a highly conductive solid is quasi-isothermal and this is here translated by the fact that wherever the path starts, it visits the whole solid several times and looses the information about its initial position. The expectation of the Monte Carlo weight is the same whatever the starting location, i.e. the temperature is the same. Although trivial, this feature leads to computation times that can be disastrous. Of course we can think of changing the model, stipulating the solid is isothermal from the start, as we did for perfectly mixed fluid cells, but this is not an easy task and we should rather retain that, at present stage, even such a very simple limit is enough to damage our proposal.

PERSPECTIVES

We have shown how heat transfer rendering can be addressed by using one single path statistics to make the infrared image of an object with no need to compute its temperature. For this, we have developed the theoretical framework to couple the three modes of heat transfer, starting from a model suitable for a large number of scenes, still with some constraints. The latter are our starting points for drawing perspectives:

Constraints:

(1) Solids are opaque and fluid are transparent (2) Fluid cells are perfectly mixed (3) Properties are uniform and constant in each sub-volume (4) All parameters fields are prescribed (5) Linearization assumptions have been made (e.g. radiative exchange coefficient ℎ 𝑅 assumed independent of temperature)

What should we do to extend the application range to more advanced heat transfer or to other physics? With our framework, a key observation is that the problem is expressed as a set of kinetic equations. Solutions of such PDEs are left with two kinds of parameters: sources and (kinetic) coefficients, as detailed in Sec. 6.1. Hence, the coupling can pass either by a source or by a kinetic coefficient (see Figure 10).

By posing the constraints of Sec. 3.4, the coupling can pass only by the sources. Moreover, in our radiative transfer model, the sources were linear functions of solid temperature and, for the conduction and convection model, the convective, conductive and radiative fluxes appeared linearly in the boundary conditions. Then, we have essentially designed our algorithms based on the well-known double randomization technique: when a Monte Carlo algorithm exists for one phenomenon, and when the corresponding Monte Carlo weight is a linear function of the solution of a Monte Carlo algorithm for another phenomenon, then we can branch the two algorithms to form a new Monte Carlo algorithm for the coupled problem, with paths switching from one phenomenon to the other.

Hence, our question becomes more specifically: what should we do to extend the application range, still keeping the benefits of double randomization, whether the coupling is via a source or via a kinetic coefficient? For coupling via more complex sources, the answer is straightforward: the proposal of Sec. 4 is a self-consistent framework that leads to the chaining of sub-paths and paves the way to easily incorporating additional phenomena to our heat transfer problem or addressing fully distinct physics (Sec. 6.2). By contrast, coupling by kinetic coefficient will call for additional algorithmic design, starting from existing preliminary work (Sec. 6.3).

Set of coupled kinetic equations

With the simplified model exposed in Sec. 3.4, an essential point is that, in the present paper, we illustrated the linear coupling of linear kinetics models.

In a kinetic model, the addressed quantity 𝜂 is transported at a given velocity ì 𝑐 (by a spatio-temporal transport operator T [𝜂] = 𝜕𝜂 𝜕𝑡 + ì 𝑐 • ∇𝜂, typically like in transient rendering) and reacts with an external material that pulls it back to a known kinetic source 𝜂 * with a characteristic time 1/𝜅, where 𝜅 is named here the kinetic Piaud In Sect. 6.1 The problem is expressed as a set of coupled kinetic equations, and we define 𝜂 * as the source and 𝜅 as the kinetic coefficient. Coupling can pass by the sources or by the kinetic coefficient.

(a) Sect. 6.2 Coupling via the sources. A path is sampled as if estimating 𝜂 alone (e.g. a radiative path, in red). At the end of the path the Monte Carlo weight should be 𝜂 * (e.g. the temperature at the absorption location). But 𝜂 * is only known as the solution of another physics (e.g. conduction). A new path (in blue) is then sampled as if estimating this other physics. At the end of this second path the Monte Carlo weight should be the local value of the source for this other physics. Again this source may be only known as the solution of another physics, etc. This is the chaining of sub-paths proposed in Sec. 4. In Sect. 6.3 Situation 1: uncoupled prescribed heterogeneous 𝜅-fields are considered, e.g. volume radiation in a heterogeneous cloud, to introduce null-collisions techniques. Here, locations of null-collision are represented by empty dots along the path. In a standard null-collision algorithm, collision locations are exponentially sampled using an upper bound κ of the kinetic coefficient and the acceptance probability is 𝜅/ κ. (b) Sect. 6.3 Situation 2: coupling with an independent physics via the kinetic coefficient. When 𝜅 is only known as the solution of another physics, at each collision along the main path (red dots), the collision type cannot be decided before launching a new path (in black), as if estimating 𝜅 itself. This new path produces a Monte Carlo weight 𝑤 𝜅 . The main path is then continued with this gathered information: the collision is accepted with probability 𝑤 𝜅 / κ. As 𝑤 𝜅 is random, this acceptance probability is random.

(c) Sect. 6.3 Situation 3: nonlinear coupling via the kinetic coefficient. The path launched for 𝜅 may itself deal with unknown kinetic coefficients and therefore require the further launching of new paths before it is continued, etc, leading to branching paths (three recursive branching levels are displayed, in red, black and brown). When 𝜅 depends itself upon the first level quantity (e.g. ℎ 𝑅 depends on the temperature), branching paths open the way to tackle non linear kinetics. In the particular case of one single nonlinear kinetic equation (e.g. the original Boltzmann equation), all paths would address the same quantity (and would be represented with the same color here).

coefficient:

T [𝜂] = -𝜅𝜂 + 𝜅𝜂 * .
(36) Diffusive or advecto-diffusive models are not of the kinetic type, but using the walk-on-𝛿-sphere approximation, each step is modeled using a kinetic equation (see Eq. 32 where T [𝜂] ≡ 𝜕𝜃 𝑆 𝜕𝑡 is only temporal, as for transport in a uniform field, spatialization being recovered because 𝜂 * ≡ 𝜃 * 𝑆 is an integral over the sphere). Note that when using a convective exchange coefficient to simplify the convective model, this model also gets a kinetic form. So, in practice, we were dealing with the coupling of kinetic equations.

One major feature of kinetic equations is that their solutions are always accessible as sums of sources (in the sense of Green theory, i.e. comprising boundary fluxes and initial conditions) multiplied by their propagators, each propagator involving negative exponentials of the integral of the kinetic coefficient over time or space. Formally, the solution of each kinetic equation, considered independently, writes

𝜂 (𝑡) = L exp - ∫ 𝜅 (𝑡)𝑑𝑡 𝜂 * . ( 37
)
The coupling passes via a source when 𝜂 * depends on the solution of another kinetic equation (Sec. 6.2). The coupling passes via a kinetic coefficient when 𝜅 (𝑡) depends on the solution of another kinetic equation (Sec. 6.3). In some cases, both can be present.

Coupling via the sources

Linearly coupling via the sources is already sufficient to alleviate some of the constraints listed above. We explain in supplementary material (Sec. 4) how the algorithm of Section 4 can actually be extended to recover most of the full initial heat transfer model of Section 3. For instance, had we considered volume radiation, we would still have been dealing with kinetic equations because the Radiative Transfer Equation is kinetic by nature. This alleviates constraint (1), provided that we still leave aside the nonlinearity of radiation (otherwise, see below). True advection add neither specific difficulty as far as coupling is concerned if we still assume that the velocity field is prescribed and independent of temperature (otherwise, see below). This alleviates constraint (2).

Coupling via the sources can also be used to incorporate additional physics: as sketched in Fig. 10(a), it suffices that one of the sources be a linear function of this additional physics to justify that when the first path ends at this source, it is chained with a sub-path for the new physics. Let us take the example of a three-dimensional electric flow inside a solid, and that heats the solid because of the local resistance to the electric flow (Joule effect). Assuming that we have a Monte Carlo algorithm for this flow and its associated power density, can we couple it to a Monte Carlo heat transfer algorithm such as those we illustrated in Section 4? The answer is yes because this additional coupling with electrical dissipation is via a volume source in the physics of solid heat-conduction: along the walk-on-sphere algorithm, when the volume source is needed, an electric sub-path would be initiated to fetch the missing information. All such enrichments of the phenomenology are straightforward provided that the coupling is linear and passes via the sources and from this point of view, the proposal of Sec. 4, although depicted through a practical example, may be considered as a closed formal framework.

Coupling via the kinetic coefficient

As an introduction to the additional algorithmic work required for tackling coupling via the kinetic coefficient, let us consider first Eq. 37 with regards to the third constraint (3).

If we consider subvolume properties as uniform and constant by parts, then 𝜅 is uniform along a subpath within a subvolume. In this case, the reexpression of Eq. 37 to eliminate the time integral is immediate. This is how we constructed Eq. 29 and 32. If now the field of 𝜅 is continuously heterogeneous and time dependent, the Monte Carlo algorithm solving 𝜂 must account for the fact that, in Equation 37, 𝜅 is varying along the path when addressing the integral ∫ 𝜅 (𝑡)𝑑𝑡 over the propagation path. We can here face three situations:

(1) the heterogeneous field and its time dependence are prescribed (uncoupled situation) (2) the field has to be estimated from a model which does not depend upon the heat transfer (constraint (4), coupling with an independent physics), (3) the field has to be estimated from a model in which the solution itself depends upon the heat transfer (constraint (5), nonlinear coupling)

In order to tackle heterogeneous time-dependent kinetic coefficient, practical solutions are already available based on the concept of null-collisions. This concept is indeed originally introduced as a way to stick rigorously to the Monte Carlo framework when progressing through a heterogeneous field. To take an example in radiative transfer, let us consider a volume with a heterogeneous distribution of colliders (𝜅 field) so that the exponential of the integral can not be known analytically. Let us consider null colliders in addition to the true colliders. When a photon encounters a null collider, it continues its path as if no collision had occurred so that the transfer of radiation is unchanged.

6.3.1 Situation 1: uncoupled prescribed heterogeneous 𝜅-fields. In standard practice, if the field is prescribed, the amount of null colliders can be adjusted locally so that, in total (true colliders plus null colliders), the new apparent coefficient field becomes homogeneous over the subvolume, and the difficult question of integrating along the path within the exponential just vanishes. When the field of collisions is known, Monte Carlo practice is then: set an upper bound value from the known field, and at each collision location, deciding whether the collision is true or null is made using only the knowledge of the local value of the true extinction coefficient versus the upper bound value. That would allow to alleviate constraint (3). Null collisions are sketched in Fig. 10(a), inspired of multiple scattering in clouds for instance, by plotting collisions as dots along the path, some of them leading to the path being continued as if no collision had occurred. 6.3.2 Situation 2: coupling with an independent physics. In this case, 𝜅 depends on the solution of another physics, which can be solved on its own.

Even if L is linear in Equation 37, the dependence is inside the exponential. For discussion, let us assume that in Eq. 36 transport is only in time (∇𝜂 is null and T [𝜂] = 𝜕𝜂 𝜕𝑡 ). But now the kinetic coefficient is itself modeled as the solution of a second kinetic equation. For this second equation, all we posit is that a Monte Carlo algorithm is available that estimates 𝜅 (𝑡) as the expectation of a weight 𝑊 𝜅,𝑡 :

𝜕𝜂 𝜕𝑡

= -𝜅𝜂 + 𝜅𝜂 * 𝜅 (𝑡) = E(𝑊 𝜅,𝑡 ) (38)
Using the same derivation as for Equations 29 and 34, the formal solution of this coupled model would have the form of Equation 37, the expectation of 𝑊 𝜅 appearing in place of 𝜅 inside the temporal integral and inside the exponential. Yet, as this integral is inside the exponential, and the exponential is itself inside an expectation over the random sampling of the path, the problem does not write anymore as the expectation of an expectation: double randomization does not apply directly.

Basically, what plays the role of the true extinction coefficient 𝜅 becomes unknown: it is the solution of a second kinetic equation to be solved together with the first one. This question has been investigated in [START_REF] Galtier | Radiative transfer and spectroscopic databases: A line-sampling Monte Carlo approach[END_REF]] and [START_REF] Terrée | Addressing the gas kinetics Boltzmann equation with branching-path statistics[END_REF]] and although the presentations were very much domain-specific, a common formalism can be extracted form these works. They both start with a standard null-collision algorithm. Let us assume that there exists a known upper bound ŵ𝜅 to 𝑊 𝜅,𝑡 , meaning that for each realization 𝑤 𝜅,𝑡 of 𝑊 𝜅,𝑡 at any time, 𝑤 𝜅,𝑡 < ŵ𝜅 . This upper bound is therefore also an upper bound for 𝜅 (𝑡). We can use it to define the null-collision amount:

𝜕𝜂 𝜕𝑡 = -ŵ𝜅 𝜂 + ŵ𝜅 𝜅 ŵ𝜅 𝜂 * + 1 - 𝜅 ŵ𝜅 𝜂 (39)
At this stage, all we have is more extinction (the kinetic coefficient is replaced with its upper bound) and the surplus in extinction is compensated by the fact that the new source (inside the parenthesis) is the sum of the old source and the sought value of 𝜂 (a recursion), weighted by the true-collision probability 𝜅 ŵ𝜅 and its complementary null-collision probability.

After this first-step of introducing null-collisions, [START_REF] Galtier | Radiative transfer and spectroscopic databases: A line-sampling Monte Carlo approach[END_REF]] and [START_REF] Terrée | Addressing the gas kinetics Boltzmann equation with branching-path statistics[END_REF] propose that 𝜅 may be replaced by its expectation inside the new source and that true-collision and null-collision probabilities can be made random themselves. An additional sampling is then required to decide these probabilities. Once the probabilities have been sampled, they are used to sample the Bernoulli random variable:

𝜅 ŵ𝜅 𝜂 * (𝑡) + 1 - 𝜅 ŵ𝜅 𝜂 (𝑡) = E(𝑊 𝜅,𝑡 ) ŵ𝜅 𝜂 * (𝑡) + 1 - E(𝑊 𝜅,𝑡 ) ŵ𝜅 𝜂 (𝑡) = E 𝑊 𝜅,𝑡 ŵ𝜅 𝜂 * (𝑡) + 1 - 𝑊 𝜅,𝑡 ŵ𝜅 𝜂 (𝑡) = E H 𝑡𝑟𝑢𝑒 (P) 𝜂 * (𝑡) + H 𝑛𝑢𝑙𝑙 (P) 𝜂 (𝑡) (40) with P = 𝑊 𝜅,𝑡 ŵ𝜅
and H 𝑛𝑢𝑙𝑙 (P) = 1 -H 𝑡𝑟𝑢𝑒 (P). The test H 𝑡𝑟𝑢𝑒 (P) is a Bernoulli random variable of probability P and the fact that this probability is itself a random variable simply means that the Bernoulli test is made after P is realized. In summary, true-collision and null-collision probabilities are random: before deciding the collision type, a path is sampled corresponding to the physics of 𝜅, as if estimating 𝜅 alone, the resulting Monte Carlo weight is used in place of 𝜅 to compute the collision-type probabilities and this random probability is used to continue the algorithm exactly the same way it would do in a standard null-collision algorithm.

Null collisions offer us a practical way to set up a single path statistics for physics that are coupled via a kinetic coefficient: the path is sampled through a field where the local values of the kinetic coefficient are unknown, and, at each collision, a 𝜅-subpath is launched to fetch the missing information. This alleviates constraint (4). Figure 10(b) illustrate this first level of path-branching with a main path (for 𝜂) and secondary paths (for 𝜅) launched at each collision location. In [START_REF] Galtier | Radiative transfer and spectroscopic databases: A line-sampling Monte Carlo approach[END_REF], 𝜅 was a gaseous absorption cross-section and 𝜅-subpaths were randomly visiting molecular spectroscopic databases to gather information about millions of absorption lines.

As a limit to this approach, we need to assume that there exists an upper bound to the Monte Carlo weight. This is not always possible. For instance, when there is a volume source inside a solid, the Monte Carlo weight will not only carry the temperature at the end of the path but will also include a contribution of the volume sources accumulated along all the visited location. Since here the visiting time does not have an upper bound, nor does the Monte Carlo weight. 6.3.3 Situation 3: nonlinear kinetics. In [START_REF] Terrée | Addressing the gas kinetics Boltzmann equation with branching-path statistics[END_REF], the addressed question is fluid mechanics rarefied gases, 𝜂 is the distribution function (the equivalence of radiance for molecules instead of photons), 𝜅 is the molecular collision cross-section and the approach is essentially the same.

But, from the present point of view, this contribution opens a new perspective: in short, molecules collide with molecules and 𝜅 is 𝜂 itself. This has two implications:

• the kinetic model for 𝜂 is nonlinear (the original Boltzmann equation); • when a 𝜅-subpath is launched because of missing information about the local value of 𝜅, as 𝜅 is 𝜂, the new path is practically an 𝜂-path that may itself require the launching of a 𝜅-subpath, etc.

A nonlinear kinetic model is therefore translated into a branchingpath statistics. Such statistics remain widely unknown [START_REF] Dimov | Monte Carlo algorithm for solving integral equations with polynomial non-linearity. Parallel implementation[END_REF][START_REF] Rasulov | Monte Carlo Solution of Dirichlet Problem for Semi-linear Equation[END_REF][START_REF] Rioux-Lavoie | A Monte Carlo Method for Fluid Simulation[END_REF]], but there are already cases with effective application potentials where they can be deployed successfully.

In the present context, such branching paths techniques open avenues to alleviate constraint (5), so as to consider the effect of wider temperature ranges upon radiative transfer, i.e. to take into account situations where the approximation considering ℎ 𝑅 as independent of the temperature does not hold. Such branching paths are sketched in Fig. 10: along the main path (for 𝜂), sub-paths are launched at each collision location (for 𝜅), and these sub-paths themselves initiate sub-paths for unknown kinetic coefficients, etc. Three recursive branching levels are represented. The particular case of one single nonlinear kinetic equation is when each launched sub-path is of the same type as the primary path.

FUTURE WORKS

In the introduction and related works, we highlighted that there were, up to our knowledge, only few reported attempts to address the coupling of several physics in a single path space Monte Carlo algorithm. Figure 4 illustrates very well the expected benefits of our coupling proposal, and specifically the better scalability of Monte Carlo methods over deterministic approaches to solve complex problems.

Beyond the range of cases that can be readily handled on the basis of the work presented above, we have identified two limitations that would deserve future work: the risk of path-entrapment and the existence of an upper bound when using null-collision. In addition, we believe that several improvements should also be tackled as fundamental research topics.

Brownian motion is among the leading requirements. Pragmatically, we used a walk-on-𝛿-sphere algorithm to model diffusion, in association with heterogeneous advection and heterogeneous sources. By construction, walk-on-𝛿-sphere loses general guaranties about convergence, while walk-on-sphere ensures a logarithmic convergence under homogeneous hypothesis (thanks to the tangent sphere approach). Keeping such convergence guaranties in heterogeneous situations is a wide research question in the walkon-sphere literature, but also in an active literature about walks on other shapes than spheres [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF].

A very similar comment holds for branching path statistics. If the reported practical successes are convincing, there is still a lack of a general framework, and we believe that making a connection with the active literature opened by the work of [START_REF] Ivan T Dimov | Monte Carlo methods for applied scientists[END_REF] should help for nonlinear coupling and more specifically nonlinear kinetics.

In our proposal, we did not consider importance sampling nor next event estimation but we pointed-out convergence difficulties associated with images of a well isolated house, where paths starting from an outside camera and reaching the heating system were scarce. If similar questions are already well studied in the rendering community, they remain open questions for thick diffusive media, typically for radiative transfer in systems with huge optical thicknesses. The difficulty has been well identified in the nuclear shielding literature since the early ages of Monte Carlo simulations [START_REF] Hammersley | Monte carlo methods[END_REF]], but it remains open. Several links should be explored in this still-active literature.

Sensitivity modeling is another promising topic evoked in Section 4. We computed propagators for replay, and propagators are sensitivities to the sources. We also insisted on the information carried by the paths themselves: even simply plotting them is already useful as far as analysis is concerned. Obviously, this is true for other parameters than sources (conductivities, capacities, geometry, etc). Following the spirit of differential light transport, we believe that this path information could be translated into the computation of sensitivities to these other parameters. This would be highly useful in practical engineering terms, e.g., for inverse problem solving.

Along this paper, all concepts were introduced and discussed from the heat-transfer perspective. This is how the question was introduced in [START_REF] Villefranque | The "teapot in a city": A paradigm shift in urban climate modeling[END_REF], thinking of the energy transition and climate-change communities. But even in these communities, the research questions for which computer graphics Monte Carlo simulations could be useful widely extend heat transfer. Fluid mechanics, electrical engineering, biology modeling are only few of the phenomena addressed by a broad scientific community that has the same flexibility and scalability requirements as those addressed by computer graphics. We hope that the coupling of Monte Carlo path sampling algorithms will pave the way for new fruitful connections.

A BASICS OF HEAT TRANSFER A.1 Conduction and convection models for solid and fluid volumes

We detail how Equations 8 and 12 are derived from the basic laws of thermodynamics. For continuous matter with a uniform mass density, the first principle of thermodynamics (the law of energy conservation) relates the internal energy density 𝑢 ≡ 𝑢 ( ì 𝑥, 𝑡) to the energy flux density vector ì 𝑗 ≡ ì 𝑗 ( ì 𝑥, 𝑡) as

𝜕𝑢 𝜕𝑡 + ì ∇ • ì 𝑗 = 0 (41) with 𝜕𝑢 𝜕𝑡 = 𝜌𝑐 𝜕𝜃 𝜕𝑡 ( 42
)
where 𝜌 is the mass density and 𝑐 the specific heat capacity. For an opaque solid, ì 𝑗 reduces to the conduction of heat (no advection, no radiation) and is given by Fourier's law:

ì 𝑗 = -𝑘 ì ∇𝜃 ( 43
)
where 𝑘 is the thermal conductivity. Eq. 8 is obtained by replacing Eq. 42-43 into Eq. 41.

We now consider a cavity inside a solid. The cavity is filled with a fluid. At any location ì 𝑥 inside the cavity, the fluid temperature is 𝜃 𝐹 ≡ 𝜃 𝐹 ( ì 𝑥, 𝑡). At any location ì 𝑦 at the surface of the solid interface, the solid temperature is 𝜃 𝑆 ≡ ( ì 𝑦, 𝑡). Newton's law is used to replace the detailed description of the boundary layer of the fluid along the surface: 𝜃 𝐹 ( ì 𝑦, 𝑡) is assumed to be the temperature of the fluid outside the boundary layer ; therefore, 𝜃 𝐹 ( ì 𝑦, 𝑡) ≠ 𝜃 𝑆 ( ì 𝑦, 𝑡) (the boundary layer is replaced by a temperature discontinuity) ; then, the flux density at ì 𝑦 is :

𝜑 ( ì 𝑦, 𝑡) = ℎ( ì 𝑦)(𝜃 𝐹 ( ì 𝑦, 𝑡) -𝜃 𝑆 ( ì 𝑦, 𝑡)). ( 44
)
Choosing the value of the proportionality coefficient ℎ requires either experimental or numerical fluid mechanics (characterizing the boundary layer). We further assume that outside the boundary layer, the cavity is perfectly uniform and mixed which means that 𝜃 𝐹 is independent of ì 𝑦. We assume a constant pressure. The first principle for an open system is then the same as Eq. 41 with enthalpy instead of internal energy is where 𝒮 is the fluid boundary surface. As for Eq. 42, the enthalpy relates to the isobaric specific heat capacity 𝑐 𝑝 :

𝑑𝐻 𝑑𝑡 = ∫ V 𝜌𝑐 𝑝 𝜕𝜃 𝐹 𝜕𝑡 𝑑 ì 𝑥 = V𝜌𝑐 𝑝 𝑑𝜃 𝐹 𝑑𝑡 (46)
where V denotes the fluid volume. Eq. 12 is obtained by replacing Eq. 44 and 46 into Eq. 45.

A.2 Finite difference on the Laplace operator

We detail the derivation of Eq. 32 from Eq. 9. Each term of the Laplace operator Δ𝜃 𝑆 := 𝜕 2 𝜃 𝑆 𝜕𝑥 2 + 𝜕 2 𝜃 𝑆 𝜕𝑦 2 + 𝜕 2 𝜃 𝑆 𝜕𝑧 2 is approximated by two finite differences as in any cartesian coordinate system, e.g. As the Laplace operator is independent of the coordinate system, we can apply any rotation and then average the source on the 𝛿-sphere which we denote 𝜃 * 𝑆 .

B ALGORITHMS

We provide the pseudo-codes for the method presented in Sec. 4. The sub-paths construction for each physical mode are given in separate algorithms, respectively Algo. 1 for convection, 2 for radiative transfer and 3 for conduction. The selection of a physical mode occurs at boundaries, as shown in Algo. 4.

  Fig. 1. We propose a Monte Carlo approach to tackle multiple physics with a single algorithm, translating their coupling into a single path-space composed of randomly chained sub-paths for each physics. Application is exemplified with heat transfer. (a) An infrared image of a steady state thermal exchanger, with temperature imposed on the left and right walls. (b) Monte Carlo paths alternate between heat-transfer modes (here conduction and radiation). (c) A huge benefit is the fast production of transient simulations, at any time, using the information gathered in (a), i.e. from only one Monte Carlo run at steady state.

  Let us start with pure radiative transfer, typically addressing the Monte Carlo estimation of the spectrally integrated radiance 𝐿 ≡ 𝐿( ì 𝑥, ì 𝜔) at the location ì 𝑥 of an infrared camera inside the fluid, as in Fig. 2. A standard russian roulette algorithm leads to the backward sampling of a multiple reflection random path Γ ending at location , Vol. 42, No. 4, Article . Publication date: August 2023. ì 𝑥 Γ on the solid-fluid interface, the Monte Carlo weight being the infrared source at this location, i.e. 𝑊 Γ = 𝜎𝜃 4 𝑆 ( ì 𝑥 Γ ) 𝜋

,

  Vol. 42, No. 4, Article . Publication date: August 2023.

,

  Vol. 42, No. 4, Article . Publication date: August 2023.

,

  Vol. 42, No. 4, Article . Publication date: August 2023.

Fig. 4 .

 4 Fig. 4. Computation times insensitive to geometric complexity. Images of buildings in a cold environment with air-conditioning in each room. Although buildings thermally interact with each other via the ground or because of infrared exchanges between external walls, computation time is the same for 1, 9, 81 and 961 buildings (small differences are only due to the number of building pixels).

  Fig. 5. (top) Infrared renderings of a farm being insulated by (left) the interior or (right) the exterior. The boundary conditions are the heater at 300K inside the lower room, the ground at 286K at 3m depth, the surrounding air temperature at 280K and the incident radiance temperature at 280K for downward directions (coming from the upper atmosphere assuming a foggy environment i.e. equal to the air temperature). (below) Example paths: starting with convection from a probe inside the lower room, 3 paths in first picture, 60 in second picture. For each path, radiative subpaths are in red, conductive subpaths in blue and convective subpaths in gray. Conductive subpaths are successions of straight lines between 𝛿-sphere centers (see Sec. 4.2). As internal air is perfectly mixed, there is no explicit heat migration within the fluid and convective subpaths are only sketched as straight lines from the surface location where convection is started to the surface location where it ends.

Fig. 7 .

 7 Fig. 7. Validation of the propagator implementation (triangles) against a Monte Carlo simulation (circles) on the computation of the lower-room probe temperature and its uncertainty (as in Fig.6), which are similar in both cases.

Fig. 8 .

 8 Fig.8. Typical of a Monte Carlo resolution, the computation time (blue) of our coupled method is linear wrt the number of samples 𝑁 and the normalized variance (i.e. v 𝑎𝑟 = 𝑣𝑎𝑟 /(𝑚𝑒𝑎𝑛 -280) in red) wrt 1/ √ 𝑁 . Considering a farm insulated by the interior, two probes are set (left) inside the lower room (triangle) and at the external wall (circle). When the probe is outside the farm, the relative variance is far bigger because the source of heat is hard to reach (due to the insulation) and the computation time is shorter because most paths stop directly in the known radiative environment.

Fig. 10 .

 10 Fig.10. Illustrations of the different cases considered in Sec. 6. In Sect. 6.1 The problem is expressed as a set of coupled kinetic equations, and we define 𝜂 * as the source and 𝜅 as the kinetic coefficient. Coupling can pass by the sources or by the kinetic coefficient. (a) Sect. 6.2 Coupling via the sources. A path is sampled as if estimating 𝜂 alone (e.g. a radiative path, in red). At the end of the path the Monte Carlo weight should be 𝜂 * (e.g. the temperature at the absorption location). But 𝜂 * is only known as the solution of another physics (e.g. conduction). A new path (in blue) is then sampled as if estimating this other physics. At the end of this second path the Monte Carlo weight should be the local value of the source for this other physics. Again this source may be only known as the solution of another physics, etc. This is the chaining of sub-paths proposed in Sec. 4. In Sect. 6.3 Situation 1: uncoupled prescribed heterogeneous 𝜅-fields are considered, e.g. volume radiation in a heterogeneous cloud, to introduce null-collisions techniques. Here, locations of null-collision are represented by empty dots along the path. In a standard null-collision algorithm, collision locations are exponentially sampled using an upper bound κ of the kinetic coefficient and the acceptance probability is 𝜅/ κ. (b) Sect. 6.3 Situation 2: coupling with an independent physics via the kinetic coefficient. When 𝜅 is only known as the solution of another physics, at each collision along the main path (red dots), the collision type cannot be decided before launching a new path (in black), as if estimating 𝜅 itself. This new path produces a Monte Carlo weight 𝑤 𝜅 . The main path is then continued with this gathered information: the collision is accepted with probability 𝑤 𝜅 / κ. As

  42, No. 4, Article . Publication date: August 2023.

.

  + 𝛿) -𝜃 (𝑥) 𝛿 -𝜃 (𝑥) -𝜃 (𝑥 -𝛿) 𝛿 = 𝜃 (𝑥 + 𝛿) -2𝜃 (𝑥) + 𝜃 (𝑥 -𝛿) 𝛿 2 Summing these terms, Δ𝜃 𝑆 = 1 𝛿 2 (-6𝜃 𝑆 + 6 θ𝑆 ) where the source term is defined as: θ𝑆 := 1 6 [𝜃 (𝑥 + 𝛿, 𝑦, 𝑧) + 𝜃 (𝑥 -𝛿, 𝑦, 𝑧) +𝜃 (𝑥, 𝑦 + 𝛿, 𝑧) + 𝜃 (𝑥, 𝑦 -𝛿, 𝑧) +𝜃 (𝑥, 𝑦, 𝑧 + 𝛿) + 𝜃 (𝑥, 𝑦, 𝑧 -𝛿)].

  𝑥 𝛾 . The Monte Carlo weight is 𝜃 𝑆 ( ì 𝑥 𝛾 ) if the solid surface temperature is known (Eq. 19). (b) When the solid surface temperature is not known, one of the three heat transfer modes is randomly selected (Eq. 23). If conduction is selected, the Monte Carlo weight is 𝜃 𝑆 ( ì 𝑥 𝛾 -𝛿 ì 𝑛), i.e. the temperature of the solid at an injection location 𝛿-away from the interface. If convection is selected, the Monte Carlo weight is the fluid temperature 𝜃 𝐹 . If radiation is selected, a direction ì 𝜔 is sampled and the Monte Carlo weight is the incoming radiance temperature 𝜃 𝑅 ( ì 𝑥 𝛾 , ì 𝜔).

	𝜃𝑆 ( ì 𝑥𝛾 ) unknown
			𝜃𝑆 ( ì 𝑥𝑤𝑜𝑠 )	𝜃𝑆 ( ì 𝑥𝑤𝑜𝑠 -𝛿1 ì 𝑛)	𝛿1
					𝛿2
	𝛿 𝑥𝛾 -𝛿 ì 𝜃𝑆 ( ì 𝑛)	ì 𝑥𝛾 -𝛿 ì 𝑛	𝜃𝑆 ( ì 𝑥𝑤𝑜𝑠 + 𝛿2 ì 𝑛)	𝜃𝑆 ( ì 𝑥𝐹 )
	ì 𝜔					𝜃𝑆 ( ì 𝑥 ′ 𝐹 )
		𝜃𝐹
	𝜃𝑅 ( ì 𝑥𝛾 , ì 𝜔 )					𝜃 (known2 ) → end
	𝜃𝑆 ( ì 𝑥 ′ 𝛾 ) known → end						𝜃 (known1 ) → end
	(a) Eq. 19	(b) Eq. 23	(c) Eq. 24	(d) Eq. 26	(e) Eq. 28	(f) Summary
	Fig. 2. Construction of a coupled path. (a) The path starts with a radiative sub-path 𝛾 initiated at the camera, potentially reflecting at the solid surface until
	absorption at ì
	, Vol. 42, No. 4, Article . Publication date: August 2023.

1

  Monte Carlo at 𝑡 obs = +∞ Fig.3. A stationary Monte Carlo estimation of a probe (triangle) temperature (left) is a weighted sum of the sources, here the known exterior (dot) and heater (cross) temperatures. The propagator enables to evaluate the samples under different source temperatures (top) shown by different colors, without performing a new simulation. The observation time can also be changed (bottom) for an extra-memory cost, then the initial temperature also becomes a source (squares).

		change
		sources
		change
		𝑡 𝑜𝑏𝑠
	probe	𝜃 ext	𝜃 heat	𝜃 (𝑡 𝐼 )	𝜃 ′ ext	𝜃 ′ heat

  290K (circle) or 300K (triangle), the temperature is similar for both insulations at stationary states. However, it increases much faster when the insulation is inside, resulting in a better thermal comfort. One single Monte Carlo simulation is performed at the steady state with 10000 samples. The stored information then allows to very rapidly replay the temperature estimation (see Fig.1), now in transient regime at any time.

				1 min	1 hour	1 day	1 month
		294 296	Insulation inside outside	Heater at 290K 300K
		292
	temperature (K)	282 284 286 288 290	Turn on heater at t = 0				Compute propagator at t =
		280
			10 0	10 2	10 4 time (s)		10 6	10 8
	Fig. 6. Temporal evolution (in s) of the lower room temperature (in K) of a
	farm being insulated by inside (blue) or outside (orange). The heater being
	280 282 284 286 288 290 292 set to 10 0 temperature (K) Turn on heater at t = 0 294 296 Insulation inside outside	10 2 1 min Computation Monte-Carlo Propagator	10 4 time (s) 1 hour	1 day	10 6 1 month	10 8 propagator Compute at t =

, Vol. 42, No. 4, Article . Publication date: August 2023.

  𝜏 𝑏,𝑐𝑜𝑛𝑣 ← exponential sampling with a mean 𝜏 𝑐𝑜𝑛𝑣 ; 𝑡 ← 𝑡 -𝜏 𝑏,𝑐𝑜𝑛𝑣 ; /* Time rewind */ if 𝑡 ≤ 𝑡 𝐼 then return 𝜃 𝐹 (𝑡 𝐼 ) ; Conductive sub-path Input: A position ì 𝑥 and a time 𝑡 Output: A temperature 𝜃 Loop ì 𝜔 ← isotropic sampling of a direction; 𝑠 + ← Trace a ray from ì 𝑥 along + ì 𝜔; 𝑠 -← Trace a ray from ì 𝑥 along -ì 𝜔; δ ← 𝑚𝑖𝑛(𝑚𝑖𝑛(𝛿, 𝑠 + ), 𝑚𝑖𝑛(𝛿, 𝑠 -)); 𝜏 𝑏,𝑐𝑜𝑛𝑑 ← exponential sampling with a mean 𝜏 𝑐𝑜𝑛𝑑 ; 𝑡 ← 𝑡 -𝜏 𝑏,𝑐𝑜𝑛𝑑 ; /* Time rewind */ if 𝑡 ≤ 𝑡 𝐼 then return 𝜃 𝑆 (𝑡 𝐼 ) ; Vol. 42, No. 4, Article . Publication date: August 2023. 20 • Mégane Bati, Stéphane Blanco, Christophe Coustet, Vincent Eymet, Vincent Forest, Richard Fournier, Jacques Gautrais, Nicolas Mellado, Mathias Paulin, and Benjamin Piaud Algorithm 4: Chaining selection at boundaries Input: A position ì 𝑥 Output: A temperature 𝜃 if ì 𝑥 is on a solid-solid interface then 𝑃 𝑐𝑜𝑛𝑑,1 ← 𝑃 𝑐𝑜𝑛𝑑 ← 1 -𝑃 𝑐𝑜𝑛𝑣 -𝑃 𝑟𝑎𝑑 ; 𝑟 ← uniform sampling in [0, 1); if 𝑟 < 𝑃 𝑐𝑜𝑛𝑣 then return 𝜃 𝐹 ; /* Algo. 1 */ else if 𝑟 < 𝑃 𝑐𝑜𝑛𝑣 + 𝑃 𝑐𝑜𝑛𝑑 then return 𝜃 𝑆 ( ì 𝑥 -𝛿 ì 𝑛) ; /* Algo. 3 */ else if 𝑟 < 𝑃 𝑐𝑜𝑛𝑣 + 𝑃 𝑐𝑜𝑛𝑑 + 𝑃 𝑟𝑎𝑑 then ì 𝜔 ← cosine weighted sampling around ì

					Algorithm 2: Radiative sub-path
					Input: A position ì 𝑥 and a direction ì 𝜔
					Output: A temperature 𝜃
					Loop
	𝑘 1 𝛿 1 + 𝑟 ← uniform sampling in [0, 1); 𝑘 1 𝛿 1 𝑘 2 ; 𝛿 2	𝑠 ← Trace a ray from ì 𝑥 along ì 𝜔; if 𝑠 = ∞ then return 𝜃 𝑅,𝑎𝑚𝑏𝑖𝑒𝑛𝑡 ;
	if r < 𝑃 𝑐𝑜𝑛𝑑,1 then return 𝜃 𝑆 ( ì 𝑥 -𝛿 1 ì 𝑛) ; else	/* Algo. 3 */	else ì 𝑥 ← ì 𝑥 + 𝑠 ì 𝜔; 𝑟 ← uniform sampling in [0, 1);
	return 𝜃 𝑆 ( ì 𝑥 + 𝛿 2 ì 𝑛) ; end	/* Algo. 3 */	if 𝑟 < 𝜖 ( ì 𝑥) then return 𝜃 𝐵 ( ì 𝑥) ;	/* Algo. 4 */
	else if ì 𝑥 is on a solid-fluid interface then 𝑃 𝑐𝑜𝑛𝑣 ← ℎ ( ì 𝑥) 𝑘 𝛿 +ℎ ( ì 𝑥)+ℎ 𝑅 ;	else ì 𝜔 ← Sample a direction wrt 𝑝 𝑟 at ì 𝑥; end
	𝑃 𝑟𝑎𝑑 ←	ℎ 𝑅 𝛿 +ℎ ( ì 𝑘 𝑥)+ℎ 𝑅	;	end EndLoop
	return 𝜃 𝑅 ( ì 𝑥, ì 𝜔) ; end end	𝑛; /* Algo. 2 */	Algorithm 3: 𝜏 𝑐𝑜𝑛𝑑 ←	6𝑘 ; 𝜌𝑐 δ2
					/* Initial temperature */
					else
					ì 𝑥 ← ì 𝑥 + δ ì 𝜔;
	Algorithm 1: Convective sub-path	if ì 𝑥 is on boundary then
	Input: A time 𝑡	return 𝜃 𝐵 ( ì 𝑥) ;	/* Algo. 4 */
	Output: A temperature 𝜃	end
	𝜏 𝑐𝑜𝑛𝑣 ←	𝜌𝑐 V ℎ𝒮	;	end EndLoop
					/* Initial temperature */
	else
	ì 𝑥 𝒮 ← sampling of ì 𝑋 𝒮 proportionally to ℎ;
	return 𝜃 𝐵 ( ì 𝑥 𝒮 ) ;	/* Algo. 4 */
	end

,

, Vol. 42, No. 4, Article . Publication date: August 2023.

ACKNOWLEDGMENTS

This work was supported by the ANR under grand No.: ANR-18-CE46-0012 (https://anr.fr/Project-ANR-18-CE46-0012).