Supplementary Figures

Supplementary Figure 1: Random Light Induced and DLVO interaction potentials. Calculated interaction potential for charge stabilized colloidal particles of size 2R = 2 μm suspended in water in the presence and the absence of the random light (rl) induced interactions. The combination of van der Waals attractions (vdW red line), Hamaker constant A = 0.1 kBT , and electrostatic double layer repulsions (dl-blue line), Debye length λD = 6 nm, contact potential Udl(D = 0) = 44 kBT leads to the DLVO potential (black line). The attractive potential due to random light forces (green line), monochromatic random illumination with a contact potential of -2 kBT , equation (6) is superimposed leading to the total potential (magenta line). The repulsive part dominates at very short distances and in turn this allows us to probe the superimposed light induced attractions, without particles sticking together irreversibly.


D (µm)
Double-layer repulsion U dl Van der Waals attraction U vdW DLVO potential U dl +U vdW Random light potential U rl Total potential U dl +U vdW +U rl  inside two identical optical traps are positioned at a distance r 0 . For an isolated particle the distribution of positions is Gaussian with a standard deviation σ, which is set by the trap stiffness. Thus 2σ is a measure for the typical distances r − r 0 probed by thermal motion. Interactions between particles lead to a characteristic change in the distribution of particle positions. Attractive interactions increase the probability for the particles to approach. Precise measurements of f pair (r) are therefore a sensitive tool to determine the particle-particle interaction potential U (D = r − 2R).

Supplementary Notes
Supplementary Note 1: Field-Field correlations and cross-spectral density in a stationary random field We consider a fluctuating electric field, E(r, t) in a transparent and non-dispersive homogeneous medium with real refractive index n h = √ h . For a stationary field 1 , the spatiotemporal fluctuations, are characterized by the the cross-spectral density tensor W ij (r, r , ω) given by 2 where G ij (r, r , ω) are the matrix elements of the Green tensor, (I 3 is the identity tensor) and u E (r, ω) is related to the time-averaged electric energy per unit volume, Supplementary Note 2: Multiple scattering between two compact bodies We consider a particle A, centered at the origin of coordinates, and a particle B displaced a distance r along the positive z-axis. From a physical point of view, instead of a continuous approach, each particle can be seen as made of discretized, N A and N B , identical cubic elements of volume v. This is also known as a discrete dipole approach (DDA) 3 . In the presence of an external polarizing field, E inc (r, ω), each volume element acts as an induced dipole proportional to the polarizing field, i.e.
where α(ω) is the polarizability which, for cubic or spherical elements of volume v, is given by 4 The polarizing field on a given element in particle B, E inc (r B n , ω), is given by the solution of the multiple scattering problem: (and an equivalent equation for particle A). For simplicity in the notation, we do not include here the explicit ω-dependence. These are a set of 3N A + 3N B equations that can be written in compact matrix form as Introducing the T-matrix, defined as where I 3N B is the 3N B × 3N B identity matrix (and an equivalent expression for T A ), the formal solution of the scattering problem can be written as Our approach can be seen as the DDA-like version of the well known T-matrix approach of multiple scattering of electromagnetic waves by two different objects usually described in terms of a basis of multipolar vector wave functions (see for example refs 5,6).

Supplementary Note 3: Optical interactions between two compact bodies induced by random light fields
In the presence of a random (stationary) field, E 0 (r, t), both the dipoles and the polarizing fields are, in general, fluctuating quantities and the time averaged force along the z-axis may be written as the sum of two different terms (see for example, ref. 7) where the first term describes the force induced by the fluctuating (external) field, E fluc 0 with the corresponding induced dipole p ind as discussed in the main text. The second involves the (spontaneous and thermal) fluctuations of the dipole p fluc .
We focus on lossless particles and discard the second contribution (in absence of absorption, there are no spontaneous and thermal fluctuations of the dipoles). From equations (5) and (13), the total time-averaged force on particle B is then given by where E inc (r n , ω) is given in equation (7) and the gradient of the incoming field is the sum of three different terms These three terms give three different contributions to the total force on particle B. The first term, F B1 z , can be seen as coming from the homogeneous radiation field on particle B (which after arriving at B, suffers multiple scattering events with particle A). The second, F B2 z , comes from the radiation first scattered by particle A. The last term, arising from the multiple scattering interactions inside the particle, does not contribute to the total force on B since these interactions cancel out when summing over all the dipoles in B after averaging over the random field. Taking into account that where "Tr" stands for the trace of the 3N B × 3N B matrix. After some algebra and taking into account that in absence of absorption (i.e. (ω) andα 0 are real) the second contribution can be shown to be given by Adding equations (18) and (20) we finally obtain that, in absence of absorption, the total force is conservative F B z = −∂U (r)/∂r with an interaction potential given by The dependence of the interaction on distance is completely contained in ↔ G B,A whereas all the shape and material dependence is contained in the T-matrices.
In the case of equilibrium thermal blackbody radiation the electric energy density, U E (ω), is which, at zero temperature gives u E (ω) = k 3 /(4π 2 ). For absorbing (emitting) particles in equilibrium, we can include in equation (13) the contribution of the fluctuating dipoles and the corresponding radiated fields 7 (linked through the fluctuation-dissipation theorem). Interestingly, in equilibrium, this additional contribution conspires with the force due to the field fluctuations to give a total interaction potential which is exactly given by equation (21) distance D becomes of the order of the particle radius; a detailed analysis of these interactions will be described elsewhere). When the optical response of the particles can be described by their electric and magnetic polarizabilities, α e n (ω) and α m n (ω) respectively (n = A, B). The presence of an external polarizing field induce both electric, p and magnetic, m, dipoles, i.e. p(r n , ω) = 0 h α e n (ω)E inc (r n , ω) where Z ≡ µ 0 /( 0 h ) is the impedance of the homogeneous medium. The polarizabilities are simply related to the first electric, a 1 , and magnetic, b 1 Mie coefficients: When we just consider dipolar particles we can write being G E,x (r) = G E,y (r) = 1 + i kr − 1 k 2 r 2 g(r) (28) and g(r) = e ikr /(4πr) the scalar Green function.
In absence of absorption the trace formula for the interaction potential [equation (21)] can be calculated in closed form as: As shown in Fig. 2 in the main text, for two identical particles near the first Mie dipolar magnetic resonance the interaction force can be repulsive in analogy with the repulsive interactions between resonant atoms 12 .
Supplementary Note 5: Random light forces between dipolar electric particles: gravitational like interactions If we consider the long wavelength limit, where the magnetic polarizability is negligible, equation (21) takes the simple form which can be shown to be equivalent to equation (11) in ref. 13 in absence of absorption.
A remarkable prediction concerning optically induced interactions between atoms, molecules or small dipolar particles 14 is that, after averaging over all orientations of the inter-atomic axis with respect to the incident beam, the interaction is an isotropic long-range, "gravitational-like", 1/r potential in the near field. It was suggested 15 that this averaging could be experimentally achieved by an isotropic external illumination by means of multiple incoherent beams which, for atomic systems, could give rise to stable Bose-Einstein condensates with unique static properties 15 . An alternative is to average over all orientations and polarizations of the incoming, uncorrelated, plane waves 13 : In the weak scattering limit, expanding equation (32) leads to which in the short distance limit gives the above mentioned long-range 1/r dependence of the optical interaction potential in agreement with previous results 13,14 . It is worth noticing that similar ideas were considered in the earlier proposal by Spitzer 16 of the so-called mock gravity, gravitylike interactions between matter in the universe due to background isotropic radiation pressure.

Supplementary Note 6: Weak scattering approximation
In the weak scattering limit, we can expand equation (21) It is easy to see that this is the result obtained at lowest order in [the so-called Born (or Rayleigh-Gans-Debye) approximation]. For two identical spheres of radius R, their centers being a distance r apart, in a quasi monochromatic random field, the interaction energy in the weak scattering limit can be seen as a Hamaker's integral 17

Laser trapping experiment and data treatment
The focused beam of a Topica DL 100 diode laser operating at a wavelength of 785 nm is time shared between two points using a galvano mirror (Galvoline G1432) driven by a square-wave oscillation at a frequency of 500 Hz. A telescope (Thorlabs, 3x Galilean optical beam expander, The thermal motion of the two particles in the adjacent traps is illustrated in Supplementary   Figure 2. For a given mean separation distance r 0 of the time shared optical traps we perform two experiments (see Fig. 1 of main text). In a first experiment we acquire a movie of 4000 images at a frame rate of 90 Hz under the influence of a random light field. In a subsequent reference experiment the random light field is turned off and the measurement is repeated under otherwise identical conditions. The recorded sequence of images is analysed using a standard particle tracking algorithm 20 to obtain distributions of the center-to-center separations of the trapped particles for both experiments. We follow the approach of Grier and coworkers 22 to obtain an autocali-brated measurement of the colloidal interaction potential by analysing the differences between the distributions in the presence and absence of optically induced forces.
where f rl and f 0 are the corresponding distributions of center-to-center separations.
From n measurements of the center-to-center separation r we compute the pair distribution f pair (r) using the technique of nonparametric density estimation 22,23 : where r i reflects the separation distance determined from one image i (one measurement); h opt is a smoothing parameter. The estimator's kernel J[(r − r i )/h opt ] can be any smooth function that satisfies the following conditions: (i) continuous and symmetric around zero (ii) integrable with its maximum J max at zero and (iii) normalized and non-negative 23 . For convenience we choose a Gaussian function of the form: The smoothing parameter h opt reflects the kernel's bandwidth. A proper choice of h opt is crucial.
A too large width obscures features in the pair distribution f pair (r) whereas a too small width yield noisy results. A good trade-off is given by Silverman's rule 23 : h opt = [4/(3n)] 1/5 σ r where σ r is the standard deviation of all separation distances r i . The benefit of nonparametric density estimation over histograms is (i) the convergence speed; for n data points the statistical error in histograms decreases as n −1/2 whereas for nonparametric density estimation the error improves as n −4/5 (see refs 22,24). More importantly the nonparametric density estimation does not rely on the choice of discrete bins.

Total interaction potential of the charge stabilized microspheres
Particles suspended in water involve both van der Waals and double layer (dl) electrostatic repulsive interactions which in combination can be described by the well known DLVO theory 25 .
The latter is dominantly repulsive for stable suspensions and thus prevents particle coagulation.
Equally, in our measurements, this repulsive part dominates at very short distances and in turn this allows us to probe the superimposed light induced attractions, without particles sticking together irreversibly. For illustration we show in Supplementary Figure 1 a typical DLVO potential for micron sized particles (2R = 2 µm) and the superimposed attraction due to random light fields corresponding to the case shown in Fig. 3a in the main text. Exact values for the Hamaker constant and the contact potential U dl (D = 0) are not known and we have chosen reasonable estimates consistent with the observed stability of the melamine particles.

Amorphous turbid layer
The turbid layer at the entry of the light filled cavity is composed of a dense amorphous assem- is then projected on a screen, imaged and analyzed with the digital camera. From this we obtain an estimate for the diffuse total transmission of T diff ∼ 1/3.