Jon Asier Bárcena-Petisco

Enrique Zuazua
email: enrique.zuazua@fau.de

Averaged dynamics and control for heat equations with random diffusion

Keywords: Averaged controllability, averaged observability, observability, random heat equation AMS subject classification: 35K05, 35R60, 93B07, 93E99

come

Introduction

We analyze the problem of controlling the averaged value of the heat equation with random diffusion. This problem is relevant in applications in which the control has to be chosen independently of the random value, in a robust way. This problem has been studied in the literature in bounded domains and with diffusivities independent of the space and time variables and with a strictly positive minimum common to almost every realization. Notably, in [START_REF] Zuazua | Stable observation of additive superpositions of Partial Differential Equations[END_REF] and [START_REF] Lü | Averaged controllability for random evolution Partial Differential Equations[END_REF] the authors consider diffusivities which follow the uniform and exponential probability distributions respectively, whereas a more general study is done in [START_REF] Coulson | On average controllability of random heat equations with arbitrarily distributed diffusivity[END_REF]. In those papers it is shown that, under those assumptions, the averaged dynamics inherits many properties from the dynamics of the heat equation (regularity, controllability, observability, etc.), with the only notable exception of the semi-group property. This is done by considering the Fourier representation of the averaged solutions.

In this paper we pursue the study to diffusivities which are allowed to take any positive value. In this scenario the averaged dynamics is still analytic (see Proposition 4.1 below), and we prove that the averaged dynamics is approximately controllable provided that we have a hierarchic decay in the time variable of the different frequencies. However, the averaged dynamics may acquire a fractional nature, or an even less diffusive one, so it may not be null controllable. What determines if we can control it is how fast the density of averaging decays when the diffusivity α vanishes. In the critical threshold, which is given by all the random variables whose density functions ρ(α) decay like e -Cρα -1 for some C ρ > 0 when α → 0, the dynamics of the average is similar to the 1 2 -fractional Laplacian, which is well-known to be critical in the context of controllability of fractional diffusion processes.

The mathematical model and main results

In this paper we treat the random heat equation described by the following system:

      
y t -α∆y = g, in (0, T ) × G, y = h, on (0, T ) × ∂G, y(0, •) = y 0 , on G,

for G a domain, g a source term, h the Dirichlet boundary conditions, y 0 the initial configuration and α the diffusivity coefficient, which is a positive random variable with density function ρ (the regime in which α is allowed to take negative values is studied in Section 8). The averaged solution of (1.1) is given by: ỹ(t, x; y 0 , g, h) := +∞ 0 y(t, x; α, y 0 , g, h)ρ(α)dα.

Moreover, we can model a control f located in G 0 ⊂ G or on Γ ⊂ ∂G by posing g = f 1 G0 or h = f 1 Γ respectively.

In order to study (1.1) we rely on being the eigenfunctions of -α∆ independent of α, as this allows us to work with the Fourier representation of the averaged solution. Thus, we may not use the same techniques in more general heat equations, like: y t -div(σ(x, α)∇y) + A(x, α) • ∇y + a(x, α)y = 0.

(1.2)

Indeed, even if we assume that w → -div(σ(x, α)∇w) + A(x, α) • ∇w + a(x, α)w can be diagonalised, its eigenfunctions depend on α and, in particular, the averages of the eigenfunctions with respect to α may not form an orthogonal set. However, by studying the dynamics and controllability of (1.1) we highlight some of the most fundamental phenomena involving (1.2). The techniques presented in this paper also work for heat equations of the type:

y t -αLy = 0,
for L any self-adjoint elliptic operator of order 2 with compact resolvent.

Since studying controllability with internal or boundary controls is almost equivalent, this paper is mainly devoted to controllability with an internal control and the few differences are explained in Section 8. In addition, to study the controllability properties of (1.1) we follow the classical duality approach (see Section 7.1 for further details) and focus on the observability properties of its adjoint system, which is given by:

       -ϕ t -α∆ϕ = 0, in (0, T ) × G, ϕ = 0, on (0, T ) × ∂G, ϕ(T, •) = φ, on G.
(1.3)

To lighten the notation we work, as usual, in its time-reversed system, which is given by:

       u t -α∆u = 0, in (0, T ) × G, u = 0, on (0, T ) × ∂G, u(0, •) = φ, on G.
(1.4)

We define the average of (1.4) as: ũ(t, x; φ) := +∞ 0 u(t, x; α, φ)ρ(α)dα.

The first property of ũ that we prove is its analyticity in the time variable from (0, +∞) to L 2 (G). Next, using this together with a hierarchic decay in the time variable of the different frequencies, we obtain some unique continuation results for (1.4). Finally we determine when the averaged dynamics of (1.4) is null observable by combining the Fourier representation of the solutions of (1.4) and the monotonicity of the solutions of (1.1)

with respect to the boundary conditions.

In order to illustrate the effect of averaging in the dynamics, let us study the dynamics of (1.4) when G = R d .

As averaging and the Fourier transform commute, we work on the Fourier transform of the fundamental solution of the heat equation, which is given by: exp(-α|ξ| 2 t).

Consequently, the Fourier transform of the average of the fundamental solutions is given by: +∞ 0 exp(-α|ξ| 2 t)ρ(α)dα;

i.e. the Laplace transform of ρ evaluated in |ξ| 2 t. In particular, for r ∈ (0, 1) if ρ(α) ∼ 0 + e -Cα -r 1-r we have that:

exp(-α|ξ| 2 t)ρ(α)dα ∼ exp(-C|ξ| 2r t r ), (1.5)
when |ξ| 2 t → +∞ as shown in (2.8) below. Thus, for those density functions the averaged dynamics in R d has a fractional nature. As we are going to prove, for G bounded this is also true and we have the usual controllability and observability results of fractional dynamics (see, for example, [START_REF] Fattorini | Exact controllability theorems for linear parabolic equations in one space dimension[END_REF][START_REF] Micu | On the controllability of a fractional order parabolic equation[END_REF][START_REF] Miller | On the controllability of anomalous diffusions generated by the fractional laplacian[END_REF][START_REF] Biccari | Local elliptic regularity for the Dirichlet fractional Laplacian[END_REF][START_REF] Biccari | Controllability of a one-dimensional fractional heat equation: theoretical and numerical aspects[END_REF]); that is, (1.5) implies that the averaged unique continuation is preserved, but (1.5) preserves the null averaged observability if and only if r > 1/2, being the threshold density functions those which satisfy:

ρ(α) ∼ 0 + e -Cα -1 . (1.6)
2 Quantification of the main results

In this section we introduce the precise definition of the previously introduced observability notions, quantify the main results and give some specific examples.

To start with, we recall the definitions of the introduced observability notions:

Definition 2.1. Let G ⊂ R d be a domain and G 0 ⊂ G be a subdomain. System (1.4) is null averaged observable or null observable in average in G 0 if for all T > 0 there is a constant C > 0 such that for any φ ∈ L 2 (G):

ũ(T, •; φ) L 2 (G) ≤ C ũ(•; φ) L 2 ((0,T )×G0) . (2.1)
If (1.4) is null averaged observable, we define the cost of the null averaged observability as:

K(G, G 0 , ρ, T ) = sup φ∈L 2 (G)\{0} ũ(T, •; φ) L 2 (G) ũ(•; φ) L 2 ((0,T )×G0) . (2.2)
Definition 2.2. Let G ⊂ R d be a domain and G 0 ⊂ G be a subdomain. System (1.4) satisfies the averaged unique continuation property in G 0 if for all T > 0 the equality ũ = 0 in (0, T ) × G 0 implies that φ = 0.

To continue with, we state the precise hypotheses on ρ. For that, we focus on the Laplace transform of ρ, which also appears naturally when G is a bounded domain (see (3.1) below). We use the asymptotic notation f (s) g(s), which means that there is C > 0 such that f (s) ≥ Cg(s) for s large enough.

• To have the unique continuation we need for some r > 0 that: (2.3)

-
• To have the null observability we need (2.3) for some r > 1 2 .

• To prove the lack of null observability we need for some C > 0 and r ∈ 0, 1 2 that:

+∞ 0 e -sα ρ(α)dα e -Cs r .

(2.4)

Let us now state the main results of this paper:

• The first main result of this paper is that in many cases we have the unique continuation property:

Theorem 2.3. Let G ⊂ R d be a Lipschitz domain, G 0 ⊂ G be
a subdomain, and ρ = 1 (0,1) or ρ be a density function which satisfies (2.3) for some r > 0. Then, system (1.4) satisfies the averaged unique continuation property in G 0 .

The proof of Theorem 2.3 is given in Section 4. For the uniform distribution it relies on explicit computations of the averaged solutions, whereas for the more diffusive case it relies on the analyticity of the averaged dynamics from t ∈ (0, +∞) to L 2 (G) (see Proposition 4.1 below) and on the fact that there is some hierarchy in how the frequencies decay, a technique dating back to [START_REF] Borel | Sur les zéros des fonctions entières[END_REF].

• The second main result of this paper concerns some cases in which we do not have averaged observability:

Theorem 2.4. Let G ⊂ R d be a Lipschitz domain, G 0 ⊂ G be a subdomain such that G 0 = G and ρ be a density function which satisfies (2.4) for some C > 0 and r ∈ 0, 1 2 . Then, system (1.4) is not null observable in average in G 0 .

We know from Theorem 2.3 that the lack of observability is not caused by a lack of unique continuation.

In fact, we prove Theorem 2.4 in Section 5 by giving a sequence φ N ∈ L 2 (G) such that: • The last main result of the paper concerns some cases in which we have averaged observability:

lim N →∞ ũ(T, •; φ N ) L 2 (G) T 0 G0 |ũ(t, x; φ N )| 2 dxdt 1/2 = +∞. ( 2
Theorem 2.6. Let G ⊂ R d be a Lipschitz locally star-shaped domain, G 0 ⊂ G be a subdomain, T > 0 and ρ be a density function which satisfies (2.3) for some r > 1 2 . Then, system (1.4) is null observable in average. In addition, there are T 0 , C > 0 such that for all T ∈ (0, T 0 ] we have that:

K(G, G 0 , ρ, T ) ≤ Ce CT -(2r-1) -1 . (2.6)
We recall that the locally star-shaped domains are defined in [2, Section 3] and include all the C 2 domains.

We prove Theorem 2.6 in Section 6 by adapting the ideas of [START_REF] Miller | A direct Lebeau-Robbiano strategy for the observability of heat-like semigroups[END_REF]; that is, we use the Fourier representation and the decay properties of the averaged dynamics.

Remark 2.7. The estimate (2.6) is an upper estimate for short-time horizons. Ideally, it would also be good to have a lower bound and to precise the constant of the exponential by some geometric terms as in the heat equation (see, for instance, [START_REF] Miller | Geometric bounds on the growth rate of null-controllability cost for the heat equation in small time[END_REF], [START_REF] Miller | The control transmutation method and the cost of fast controls[END_REF], [START_REF] Tenenbaum | New blow-up rates for fast controls of Schrödinger and heat equations[END_REF], [START_REF] Ervedoza | Sharp observability estimates for heat equations[END_REF] and [START_REF] Laurent | Observability of the heat equation, geometric constants in control theory, and a conjecture of Luc Miller[END_REF]), though this problem goes beyond the objective of this work.

Example 2.8. The density functions which satisfy the hypotheses of Theorems 2.3 and 2.6 are those which decay sufficiently fast when the diffusivity vanishes. Similarly, the density functions which satisfy the hypothesis of Theorem 2.4 are those which do not decay fast enough (including those which do not decay at all) when the diffusivity vanishes. Meaningful examples include:

1. Any density function ρ such that ρ(α) = 0 for all α ∈ (0, δ ρ ) for some δ ρ > 0 satisfies (2.3) for r = 1.

Hence, all these functions satisfy the hypotheses of Theorems 2.3 and 2.6.

2. If k ∈ (0, +∞), and p and q are two polynomials such that p(α) + q(α -1 ) = 0 for some α ∈ (0, 1), the density function: 3. The density functions given by (2.7) satisfy (2.4) for r = k k+1 . Thus, if k ∈ (0, 1) they satisfy the hypothesis of Theorem 2.4.

ρ(α) = p(α) + q(α -1 ) e -α -k 1 (0,1) (α) 1 0 (p(s) + q(s -1 )) e -s -
4. The density functions ρ(α) = 1 (0,1) (α) (that is, when α is a random variable with uniform distribution in (0, 1)) and ρ(α) = e -α 1 (0,+∞) (α) (that is, when α is a random variable with exponential distribution in (0, +∞)) satisfy (2.4) for all r > 0. Indeed, any continuous density function ρ such that ρ(0) > 0 does so.

Thus, all these functions satisfy the hypotheses of Theorem 2.4.

The proofs of items 1 and 4 are straightforward. As for items 2 and 3, we can prove them by considering the asymptotic result:

1 0 α r e -sα-α -k dα ∼ s -2+2r+k 2+2k e -c k s k k+1 , (2.8)
for all r ∈ R for some c k > 0 (independent of r) when s → +∞. These asymptotic similarities can be proved with the Laplace method. In fact, we have that:

1 0 α r e -sα-α -k dα ∼ +∞ 0 α r e -sα-α -k dα = s -1+r 1+k +∞ 0 t r e -s k k+1 (t+t -k ) dt,
where we have used the change of variables α = ts -1 k+1 . Next, we can show the equivalence by using the classical Laplace method. In fact, if φ is any convex function in (0, +∞) with minimum at t, and if f is a continuous function in a neighbourhood of t with subexponential growth it is well-known the limit:

lim θ→∞ ∞ 0 f (t)e -θφ(t) dt f t e -θφ(t) θφ t 2π = 1,
which is proved by considering that the mass of the integral is concentrated on a neighbourhood of t, by using a Taylor expansion of order 2 in the exponent and the continuity of f , then extending again the integral to (0, +∞) and finally explicitly computing the Gaussian function. For a more detailed proof of (2.8) one can consult, for instance, [4, (6.4.35) and Example 6.4.9].

The rest of the paper is organized as follows: in Section 3 we present some basic results, in Section 4 we prove Theorem 2.3, in Section 5 we prove Theorem 2.4, in Section 6 we prove Theorem 2.6, in Section 7 we resume the controllability problem, in Section 8 we comment some possible extensions, and in Appendix A we prove a technical result.

Preliminaries

In this section we introduce some basic facts and notation that we use later on. In particular, we study the spectral properties of the Dirichlet Laplacian, the size of the solutions of the heat equation and the decay implied by (2.3).

Some results on the spectral decomposition of the Dirichlet Laplacian

As usual, e i denotes (starting at i = 0) the eigenfunctions of the Dirichlet Laplacian, λ i their respective eigenvalues and Λ λ := {i : λ i ≤ λ}. In addition, for any λ > 0, P λ denotes the orthogonal projection of L 2 (G) into e i i∈Λ λ and P ⊥ λ the orthogonal projection of L 2 (G) into e i ⊥ i∈Λ λ .

To begin with, we recall that, as shown in [START_REF] Zuazua | Stable observation of additive superpositions of Partial Differential Equations[END_REF], the Fourier representation of the averaged solution is:

ũ(t, x; φ) := +∞ 0 u(t, x; α, φ)ρ(α)dα = i∈N +∞ 0 e -αλit ρ(α)dα φ, e i L 2 (G) e i (x). (3.1)
Next, we recall that the eigenvalues have a growth limited by Weyl's law:

Lemma 3.1 (Weyl's law). Let G ⊂ R d be a Lipschitz domain. We have:

lim λ→∞ |Λ λ | λ d/2 = Vol(B(0, 1)) Vol(G) (2π) d .
In particular, there is C > 0 such that for all λ ≥ λ 0 :

|Λ λ | ≤ Cλ d/2 . (3.2)
Weyl's law is proved for instance in [START_REF] Ivrii | 100 years of Weyl's law[END_REF].

Finally, we recall the following elliptic result proved in [2, Theorem 3]:

Lemma 3.2 ([2]
). Let G be a locally star-shaped domain and G 0 ⊂ G a subdomain. There exists a constant C > 0 such that for all λ > 0 and {c i } ⊂ R:

i∈Λ λ |c i | 2 1/2 ≤ Ce C √ λ i∈Λ λ c i e i L 2 (G0) . (3.3)
This result is a refinement of [27, Theorem 1.2], which was a refinement of the results proved in [START_REF] Lebeau | Contrôle exact de l'équation de la chaleur[END_REF].

Some results on the heat equation

In this subsection we state some properties of the solutions of the heat equation. We first recall that their time derivative can be estimated by using the analyticity and contraction of the semigroup of the heat equation (see [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]Sections 2.5 and 5.6]) and Cauchy's integration formula:

Lemma 3.3. Let G be a bounded domain. Then, there is C > 0 such that for all k ∈ N, s ∈ R + and φ ∈ L 2 (G)
we have that:

∂ k s v(s, •; φ) L 2 (G) ≤ C k k! s k φ L 2 (G) , (3.4)
for v the solution of:

       v t -∆v = 0, in (0, T ) × G, v = 0, on (0, T ) × ∂G, v(0, •) = φ, on G. (3.5)
It is interesting to consider the solutions of (3.5) because of the identity:

u(t, x; α, φ) = v(tα, x; φ), (3.6)
for u the solution of (1.4).

Another result that we need is that the propagation of the mass, when the initial value in some subdomain is null, is exponentially slow:

Lemma 3.4. Let G be a bounded domain and let Ĝ, G 0 ⊂ G be Lipschitz domains satisfying Ĝ ⊂⊂ G \ G 0 .
Then, there are c, C > 0 such that for all φ satisfying supp(φ) ⊂ Ĝ and all T, α > 0 we have that:

u(•; α, φ) C 0 ([0,T ];L 2 (G0)) ≤ Ce -c αT φ L 2 (G) , (3.7)
for u the solution of (1.4).

We recall that A ⊂⊂ B means that A is contained in a compact set contained in B. Lemma 3.4, whose originality we do not claim, is a consequence of the comparison principle. Indeed, following for example the ideas of [10, Lemma 4], we obtain Lemma 3.4 by comparing the solutions of (1.4) with initial value ±φ to the solution of the heat equation in R d with initial value |φ|1 G , a solution which can be estimated by using its representation with the kernel of the heat equation.

Decay properties implied by (2.3)

In this subsection we show that if the density function ρ satisfies (2.3), the averaged solutions of (1.4) have a decay similar to that of the solutions of the fractional heat equation. In particular, we prove the following result:

Lemma 3.5. Let ρ be a density function which satisfies (2.3) for some r ∈ (1/2, 1]. Then, there is c > 0 such that for all λ ≥ λ 0 and t 1 , t 2 ∈ [0, 1) satisfying t 1 < t 2 we have that:

+∞ 0 e -t2λα ρ(α)dα ≤ e -cλ r (t2-t1) +∞ 0 e -t1λα ρ(α)dα. (3.8)
We recall that λ 0 is the first eigenvalue of the Laplacian.

Proof. In order to prove Lemma 3.5 we first remark that for all s ≥ 0 we have that:

- d ds ln +∞ 0 e -sα ρ(α)dα = +∞ 0 e -sα αρ(α)dα +∞ 0 e -sα ρ(α)dα ≥ c min{s r-1 , 1}. (3.9)
Indeed, (3.9) follows from (2.3) and from the continuity of s → +∞ 0 e -sα αρ(α)dα

+∞ 0 e -sα ρ(α)dα in [0, +∞). Thus, from (3.9)

we obtain for all s 1 , s 2 ≥ 0 with s 1 < s 2 the estimate:

+∞ 0 e -s2α ρ(α)dα ≤ exp -c s2 s1 min{s r-1 , 1}ds +∞ 0 e -s1α ρ(α)dα. (3.10)
Next, we fix t 1 and t 2 satisfying t 1 , t 2 ∈ [0, 1) and t 1 < t 2 , and use different approaches depending on the value of λ:

• If λ ∈ [λ 0 , t -1 2 ], from (3.10) taking s 1 = λt 1 and s 2 = λt 2 , since [s 1 , s 2 ] ⊂ (0, 1]
, we obtain that:

+∞ 0 e -t2λα ρ(α)dα ≤ e -cλ(t2-t1) +∞ 0 e -t1λα ρ(α)dα. (3.11)
In addition, since λ ≥ λ 0 we have that:

-λ ≤ -λ 1-r 0 λ r . (3.12)
Thus, from (3.11) and (3.12) we obtain (3.8) for some c > 0 and all λ ∈ [λ 0 , t -1 2 ].

• If λ ∈ [t -1 1 , +∞), from (3.10) taking s 1 = λt 1 and s 2 = λt 2 , since [s 1 , s 2 ] ⊂ [1, +∞), we obtain that: +∞ 0 e -t2λα ρ(α)dα ≤ e -cλ r (t r 2 -t r 1 )
+∞ 0 e -t1λα ρ(α)dα.

(3.13) Moreover, we consider that:

-(t r 2 -t r 1 ) = (-t r 2 + t r 1 )(t 1-r 2 + t 1-r 1 ) t 1-r 2 + t 1-r 1 = t 1 -t 2 -t r 2 t 1-r 1 + t 1-r 2 t r 1 t 1-r 2 + t 1-r 1 ≤ t 1 -t 2 t 1-r 2 + t 1-r 1 ≤ - t 2 -t 1 2 . (3.14)
We have used in the first inequality of (3.14) that t 1 < t 2 and r ∈ (1/2, 1], and in the second one that t 1 -t 2 < 0 and t 1 , t 2 ∈ (0, 1]. Thus, from (3.13) and (3.14) we obtain (3.8) for some c > 0 and all

λ ∈ [t -1 1 , +∞). • If λ ∈ (t -1 2 , t -1 1 )
, we have that:

(t 1 , t 2 ) = (t 1 , λ -1 ] ∪ (λ -1 , t 2 ).
For the time interval (λ -1 , t 2 ) we may use the result of the second case for t 1 replaced by λ and obtain that:

+∞ 0 e -t2λα ρ(α)dα ≤ e -cλ r (t2-λ -1 )
+∞ 0 e -α ρ(α)dα.

(3.15)

In addition, for the time interval (t 1 , λ -1 ] we may use the result of the second case for t 2 replaced by λ and obtain that:

+∞ 0 e -α ρ(α)dα ≤ e -cλ r (λ -1 -t1) +∞ 0 e -t1λα ρ(α)dα. (3.16)
Thus, from (3.15) and (3.16), by taking smaller constants c if necessary, we get (3.8) for all λ ∈ (t -1 2 , t -1 1 ).

Remark 3.6. The first and third cases in the proof of Lemma 3.5 might be empty depending on the values of λ 0 , t 1 and t 2 . However, since we use this result for t 1 and t 2 arbitrarily small, we need to prove those cases.

In a similar way, we can also prove the following result:

Lemma 3.7. Let ρ be a density function which satisfies (2.3) for some r > 0. Then, there is c > 0 such that for all λ, λ such that λ > λ ≥ λ 0 and t ∈ [1, +∞] we have that:

+∞ 0 e -t λα ρ(α)dα ≤ e -ct r (( λ) r -λ r ) +∞ 0 e -tλα ρ(α)dα.

(3.17) Indeed, integrating both sides of (2.3) in (tλ, t λ) we find (3.17).

Finally, we underline the following result:

Remark 3.8. A consequence of (3.1) is that t → ũ(t, •; φ) L 2 (G)
is a decreasing function. Indeed, we have that:

ũ(t, •; φ) 2 L 2 (G) = i∈N +∞ 0 e -αλit ρ(α)dα 2 | φ, e i | 2 ,
which is a series of decreasing functions.

Unique continuation property for averaged solutions

In this section we prove the unique continuation property for averaged solutions (Theorem 2.3). We first study the uniform distribution and then the density functions which satisfy (2.3).

Proof of Theorem 2.3 for the uniform distribution

Let us compute the averaged solutions of (1.4) when α has the uniform distribution in (0, 1). For that, as in [42, Section 3] and [29, Section 3], we present ũ as the difference of two terms of known nature:

ũ(t, x; φ) = i∈N 1 0 e -λiαt φ, e i e i (x)dα = 1 t i∈N 1 λ i φ, e i e i (x) - i∈N e -λit λ i φ, e i e i (x) = 1 t -∆ -1 φ + i∈N e -λit ∆ -1 φ, e i e i (x) . (4.1)
Consequently, from T 0 G0 |ũ(t, x)| 2 = 0 and (4.1) we find that:

-∆ -1 φ + i∈N e -λit ∆ -1 φ, e i e i (x) = 0 in (0, T ) × G 0 ,
which differentiating in time implies that:

i∈N e -λit φ, e i e i (x) = 0 in (0, T ) × G 0 . (4.2)
Hence, using the analyticity of the solutions of the heat equation we have that (4.2) implies that φ = 0, and thus we have the averaged unique continuation property.

Proof of Theorem 2.3 for density functions which satisfy (2.3)

The proof consists on several steps. First, we show that assuming (2.3) the averaged dynamic are real-analytic and then use this to prove the unique continuation. To begin with, we prove the analyticity:

Proposition 4.1.
Let G be a Lipschitz domain, α any random variable and φ ∈ L 2 (G). Then, the function:

U : t ∈ (0, +∞) → ũ(t, •; φ) ∈ L 2 (G),
is analytic.

Proof. In order to prove Proposition 4.1, we prove that U ∈ C ∞ and that:

∀a 1 , a 2 ∈ (0, +∞) ∃C > 0 : sup t∈[a1,a2] U (i) (t) L 2 (G) ≤ C i i! ∀i ∈ N, (4.3)
which is a characterization of analyticity in R + (see, for instance, [20, Proposition 1.2.12]). Since:

U (t) = +∞ 0 v(αt, •; φ)ρ(α)dα,
for v the solution of (3.5), we can easily see that:

U (i) (t) = +∞ 0 α i ∂ i t v (αt, •; φ)ρ(α)dα, (4.4)
and thus U ∈ C ∞ . Moreover, (4.3) follows from (4.4), the triangular inequality and (3.4).

To continue with, we present the following auxiliary result: Lemma 4.2. Let G be a domain, G 0 ⊂ G be a subdomain and v be an analytic function from (0, +∞) to L 2 (G).

Then, if there is

T * > 0 such that v = 0 in (0, T * ) × G 0 , v = 0 in (0, +∞) × G 0 .
Proof. Lemma 4.2 follows from the fact that if v is analytic from (0, +∞) to L 2 (G) and ψ ∈ L ∞ (G), then ψv is analytic from (0, +∞) to L 2 (G). Indeed, since v is analytic, by definition, for all t ∈ (0, +∞) there

are v i, t ∈ L 2 (G) such that v = ∞ i=0 v i, t(t -t) i in a neighbourhood of t, so ψv = ∞ i=0 (ψv i, t)(t -t) i in
that neighbourhood, and thus ψv is analytic. This implies that v1 G0 is analytic from (0, +∞) to L 2 (G).

Consequently, since v1 G0 = 0 in (0, T * ), by analyticity v1 G0 = 0 in (0, +∞) × G, so v = 0 in (0, +∞) × G 0 . Now we are ready to prove Theorem 2.3:

End of the proof of Theorem 2.3. Let φ ∈ L 2 (G) such that ũ(t, x; φ) = 0 in (0, T ) × G 0 .
-λ r 0 )t r φ L 2 (G) . (4.5)
Thus, by considering (4.5) for large values of t we obtain that φ, e 0 L 2 (G) e 0 = 0 in G 0 , which by Lemma 3.2 implies that φ, e 0 L 2 (G) e 0 = 0, arriving at a contradiction.

To continue with, we can prove in an similar way that if φ is null up to the N -th frequency, then ũ(t, x; φ) = 0 in (0, T ) × G 0 implies that the (N + 1)-th frequency is also null. Consequently, we obtain by induction that ũ(t, x; φ) = 0 in (0, T ) × G 0 implies φ = 0.

Lack of null averaged observability

In this section we prove Theorem 2.4. As for the notation used in this section, C (resp. c) denotes a sufficiently large (resp. small) positive constant that may be different each time it appears and which just depends on G, G 0 , T and ρ. In particular, it does not depend on the index N that we are going to introduce.

In order to prove Theorem 2.4 we construct a sequence φ N satisfying (2.5). For that purpose, we first state and justify the properties which allow to have (2.5):

• The first property is that:

N ∈N supp(φ N ) ⊂⊂ G \ G 0 . (5.1)
This requirement is very natural as G \ G 0 is the part of the domain that cannot be observed when α = 0.

We use it in addition to Lemma 3.4 to obtain that u(t, x; α, φ N ) decays exponentially in {(x, t, α) : x ∈ G 0 , αt < N -1/2 }.

• The second property is that:

φ N ∈ e i ⊥ i∈Λ N . (5.
2)

The benefit of (5.2) is that u(t, x; α, φ N ) decays exponentially in {(x, t, α) : x ∈ G, αt ≥ N -1/2 }, which follows from (3.1).

• The third property is that for C > 0 large enough we have that:

P CN φ N L 2 (G) ≥ √ 3 φ N L 2 (G) /2. (5.3)
This estimate is needed to make sure that ũN (T,

•; φ N ) L 2 (G) / φ N L 2 (G)
does not decay too fast.

Let us construct the sequence φ N . For that, we inspire in [16, Section 6] and consider more or less a linear combination of Dirac masses; that is,

φ N ≈ C √ N i1,...,i N =1 c i1,...,i N ,N δ 0 x i 1 ,...,i N ,N . (5.4)
In fact, the Dirac masses are replaced by

C N d N ς(C N (x -x i1,...,i N ,N
)), for ς a regularizing function and C N to be defined. The property (5.1) is trivial. As for (5.2), we can obtain it by taking the right linear combination.

Indeed, we just have to solve an homogeneous linear system which, for C > 0 large enough, by Weyl's law (see Lemma 3.1) has more unknowns than equations. Finally, we can obtain (5.3) by choosing the right approximation with functions whose support has a diameter proportional to N -1/2 . In particular, we can prove that:

Proposition 5.1. Let G ⊂ R d be a Lipschitz domain and G 0 ⊂ G. Then, there is a sequence (φ N ) N ≥N0 satisfying (5.1), (5.2) and (5.3).

The rigorous proof of Proposition 5.1 is a bit technical, so it is postponed to Appendix A.

Remark 5.2. Since (5.1), (5.2) and (5.3) just depend on G and G 0 , so does the sequence φ N .

Example 5.3. In Figure 1 we illustrate the solutions of the heat equation given by the proof of Proposition 5.1 to get an insight on how they look like. For doing these graphs we have taken G = (0, π), G 0 = (0, π/2), ρ = 1 (0,1) and:

ς(x) = exp -1 10(x -1) 2 (x + 1) 2 1 (-1,1) (x).
(5.5)

We recall that in (0, π) we have that e i (x) = sin(ix) and λ i = i 2 .

Let us now prove rigorously Theorem 2.4:

Proof of Theorem 2.4. We consider φ N given by Proposition 5.1 (for N large enough). We easily find that:

T 0 G0 +∞ 0 u(t, x; α, φ N )ρ(α)dα 2 dxdt ≤ T 0 G0 +∞ 0 |u(t, x; α, φ N )| 2 ρ(α)dαdxdt = T 0 G0 +∞ 0 |u(t, x; α, φ N )| 2 1 αt≤N -1/2 (t, α)ρ(α)dαdxdt + T 0 G0 +∞ 0 |u(t, x; α, φ N )| 2 1 αt>N -1/2 (t, α)ρ(α)dαdxdt ≤ C e -c √ N + e - √ N φ N 2 L 2 (G) . (5.6)
Indeed, for the first inequality of (5.6) we have used that the L 1 norm in a probabilistic space can be estimated by the L 2 norm. As for the second inequality of (5.6), we have used (5.1) and (3.7) for bounding the first integral, whereas we have used (5.2), (3.6) and the Fourier decomposition of the solutions of the heat equation for bounding the second one. (f) The average solution with initial value φ9 zoomed in (0, π/2)

Figure 1: The averaged solutions of the heat equation when G = (0, π), G 0 = (0, π/2), ρ = 1 (0,1) with initial values of the sequence given in the proof of Proposition (see (A.1)) and with ς given by (5.5).

To continue with, using (3.1), (2.4) and (5.3) we obtain that:

ũ(T, •; φ N ) 2 L 2 (G) = i∈N ∞ 0 e -λiαT ρ(α)dα 2 | φ N , e i | 2 ≥ c i∈N e -C(λiT ) r | φ N , e i | 2 ≥ ce -CN r i∈Λ CN | φ N , e i | 2 = ce -CN r P CN φ N 2 L 2 (G) ≥ ce -CN r φ N 2 L 2 (G) . (5.7)
Hence, recalling that r ∈ [0, 1/2) we easily obtain (2.5) from (5.6) and (5.7).

Proof of null averaged observability

In this section we prove Theorem 2.6. As for the notation, C (resp. c) denotes a sufficiently large (resp.

small) positive constant that may be different each time it appears and which just depends on G, G 0 , ρ and r, but which is independent of T ∈ (0, T 0 ), for T 0 (G, G 0 , ρ) small enough.

In order to prove Theorem 2.6 we use the approach introduced in [35, Section 2]. It is not a direct consequence of the results presented in that section because the dynamics of the averaged solution only satisfies a decay property and not a semigroup property.

First, we reformulate [35, Lemma 2.1]: Lemma 6.1. Let G ⊂ R d be a domain, G 0 be a subdomain, T 0 > 0, q ∈ (0, 1) and f be a positive function such that f (t) → 0 as t → 0 + . Suppose that we have for all φ ∈ L 2 (G) and all t 1 , t 2 ∈ (0, T 0 ] satisfying t 1 < t 2 that:

f (t 2 -t 1 ) ũ(t 2 , •; φ) 2 L 2 (G) -f (q(t 2 -t 1 )) ũ(t 1 , •; φ) 2 L 2 (G) ≤ t2 t1 G0 |ũ(τ, x; φ)| 2 dxdτ. (6.1)
Then, we have for all φ ∈ L 2 (G) and T ∈ (0, T 0 ] that:

ũ(T, •; φ) L 2 (G) ≤ f ((1 -q)T ) ũ(•; φ) L 2 ((0,T )×G0) .
The proof of Lemma 6.1 is very similar to that of [35, Lemma 2.1]: a telescopic sum considering t 2,i = T q i and t 1,i = T q i+1 for i ∈ N.

As in [START_REF] Miller | A direct Lebeau-Robbiano strategy for the observability of heat-like semigroups[END_REF], we do not prove (6.1) directly, but we prove a similar version, which is the analogue of [35, Lemma

2.3]:

Lemma 6.2. Let G ⊂ R d be a domain, G 0 be a subdomain, T 0 , β, γ 1 , γ 2 , f 0 , g 0 > 0 satisfying γ 1 < γ 2 . Suppose that we have for all φ ∈ L 2 (G) and all t 1 , t 2 ∈ (0, T 0 ] satisfying t 1 < t 2 the inequality:

f (t 2 -t 1 ) ũ(t 2 , •; φ) 2 L 2 (G) -g(t 2 -t 1 ) ũ(t 1 , •; φ) 2 L 2 (G) ≤ t2 t1 G0 |ũ(τ, x; φ)| 2 dxdτ, (6.2) for f (s) ≥ f 0 exp(-2/(γ 2 s) β ) and g(s) ≤ g 0 exp(-2/(γ 1 s) β ). Then, for any γ ∈ (0, γ 2 -γ 1 ) there is T ∈ (0, T 0 ]
such that for all T ∈ (0, T ] and φ ∈ L 2 (G):

ũ(T, •; φ) L 2 (G) ≤ f -1 0 exp(1/(γT ) β ) ũ(•; φ) L 2 ((0,T )×G0) .
Moreover, if g 0 < f 0 , we can take γ = γ 2 -γ 1 and T = T 0 .

The proof of Lemma 6.2 is the same as [35, Lemma 2.3]: bounding superiorly g(s) f (qs) and using Lemma 6.1.

Now we are ready to prove Theorem 2.6. We do it by following the strategy of [35, Section 2]:

Proof of Theorem 2.6. Let t 1 , t 2 ∈ [0, 1) such that t 1 < t 2 and φ ∈ L 2 (G). We are going to prove:

ce -C((t2-t1) -1/4 + √ λ) ũ(t 2 , •; φ) 2 L 2 (G) -Ce -cλ r (t2-t1) ũ(t 1 , •; φ) 2 L 2 (G) ≤ t2 t1 G0 |ũ(τ, x; φ)| 2 dxdτ, (6.3)
for all λ ≥ λ 0 and then use Lemma 6.2 with the appropriate value of λ (depending this value on t 1 and t 2 ).

First, considering Remark 3.8 and that P λ φ ⊥ P ⊥ λ φ we have that:

ũ(t 2 , •; φ) 2 L 2 (G) ≤ 2 t 2 -t 1 t2 (t1+t2)/2 G (|ũ(τ, x; P λ φ)| 2 + |ũ(τ, x; P ⊥ λ φ)| 2 )dxdτ. (6.4)
From Lemma 3.2 and that P λ φ = φ -P ⊥ λ φ we obtain that:

2 t 2 -t 1 t2 (t1+t2)/2 G |ũ(τ, x; P λ φ)| 2 dxdτ ≤ C e C √ λ t 2 -t 1 t2 (t1+t2)/2 G0 |ũ(τ, x; P λ φ)| 2 dxdτ ≤ Ce C √ λ t 2 -t 1 t2 (t1+t2)/2 G0 |ũ(τ, x; φ)| 2 dxdτ + Ce C √ λ t 2 -t 1 t2 (t1+t2)/2 G |ũ(τ, x; P ⊥ λ φ)| 2 . (6.5)
Moreover, from the decay property of Remark 3.8 and (3.8) we have that:

Ce C √ λ t 2 -t 1 t2 (t1+t2)/2 G |ũ(τ, x; P ⊥ λ φ)| 2 dxdτ ≤ Ce C √ λ ũ (t 2 + t 1 )/2, •; P ⊥ λ φ 2 L 2 (G) ≤ Ce C √ λ-cλ r (t2-t1) ũ(t 1 , •; P ⊥ λ φ) 2 L 2 (G) . (6.6)
Thus, from (6.4)-(6.6), 2 ≤ Ce C √ λ and (t 2 -t 1 ) -1 ≤ Ce C(t2-t1) -(2r-1) -1 (recall that C > 0 is a sufficiently large constant) we obtain:

ũ(t 2 , •; φ) 2 L 2 (G) ≤ Ce C √ λ t 2 -t 1 t2 (t1+t2)/2 G0 |ũ(τ, x; φ)| 2 dxdτ + Ce C √ λ-cλ r (t2-t1) ũ(t 1 , •; P ⊥ λ φ) 2 L 2 (G) ≤ Ce C((t2-t1) -(2r-1) -1 + √ λ) t2 (t1+t2)/2 G0 |ũ(τ, x; φ)| 2 dxdτ + Ce C √ λ-cλ r (t2-t1) ũ(t 1 , •; P ⊥ λ φ) 2 L 2 (G) .
This inequality implies that, for a sufficiently small c > 0, since e -C(t2-t1) -(2r-1) -1 < 1:

ce -C((t2-t1) -(2r-1) -1 + √ λ) ũ(t 2 , •; φ) 2 L 2 (G) ≤ t2 t1 G0 |ũ(τ, x; φ)| 2 dxdτ + Ce -cλ r (t2-t1) ũ(t 1 , •; P ⊥ λ φ) 2 L 2 (G) , (6.7
) which implies (6.3).

We now define:

λ(t 2 , t 1 ) = C(t 2 -t 1 ) -(r-1/2) -1 , (6.8)
for C ≥ λ 0 a positive constant sufficiently large. If we take in (6.3) λ given by (6.8), we obtain (6.2) for the functions:

f (s) = c exp -C s -(2r-1) -1 + C 1/2 s -(2r-1) -1 , g(s) = C exp -cC r s -(2r-1) -1 .
Indeed, we have for all s ∈ (0, 1) that:

f (s) ≥ c exp -CC 1/2 s -(2r-1) -1 .
Since r > 1 2 the functions f and g satisfy the hypothesis of Lemma 6.2 for β = (2r -1) -1 , γ 1 = (cC r ) -1/β and γ 2 = (CC 1/2 ) -1/β by taking C large enough, so we end the proof by using Lemma 6.2.

The controllability problem

In this section we first resume the theoretical study of the controllability problem and then perform some simulations.

A theoretical study

As stated in the introduction, the observability results that we have obtained in this paper have some implications on the controllability of (1.1). Let us consider the controllability problem given by:

       y t -α∆y = f 1 G0 , in (0, T ) × G, y = 0, on (0, T ) × ∂G, y(0, •) = y 0 , on G. (7.1)
In particular, we focus on the following notions of controllability, which are introduced in [41]:

Definition 7.1. System (7.1) is null averaged controllable or null controllable in average if for all T > 0 there is C > 0 such that for any initial value y 0 ∈ L 2 (G) there is f ∈ L 2 ((0, T ) × G 0 ) satisfying:

f L 2 ((0,T )×G0) ≤ C y 0 L 2 (G) ,
and ỹ(T, •; y 0 , f ) = 0. If (7.1) is null averaged controllable, the cost of the null averaged controllability is defined by: K(G, G 0 , ρ, T ) = sup

y 0 ∈L 2 (G)\{0} inf f :ỹ(T,•;y 0 ,f )=0 f L 2 ((0,T )×G0) y 0 L 2 (G) . (7.2)
Definition 7.2. System (7.1) is approximately averaged controllable or approximately controllable in average if for all T > 0, ε > 0 and y 0 , y 1 ∈ L 2 (G), there exists a control f ε such that:

ỹ(T, •; y 0 , f ε ) -y 1 L 2 (G) < ε.
We now recall the duality result between observability and controllability: [START_REF] Lü | Averaged controllability for random evolution Partial Differential Equations[END_REF]). Let G ⊂ R d be a domain and G 0 ⊂ G be a subdomain. System (7.1) is null controllable in average if and only if system (1.4) is null observable in average in G 0 . In that case, K = K (see (2.2) and (7.2)); that is, the cost of the control of null averaged observability equals the cost of null averaged controllability.

Theorem 7.3 ([
Similarly, system (7.1) is approximately averaged controllable if and only if system (1.4) satisfies the unique continuation property in G 0 .

The proof of Theorem 7.3 can be found in [29, Appendix A]. As an immediate consequence we obtain that Theorems 2.4, 2.3 and 2.6 and Remark 2.5 imply the following controllability results for system (7.1):

Corollary 7.4. Let G ⊂ R d be a domain, G 0 ⊂ G be a subdomain and T > 0. Then:

• Under the hypotheses of Theorem 2.3, system (7.1) is approximately controllable in average.

• Under the hypotheses of Theorem 2.4, system (7.1) is not null controllable in average.

• If G 0 = G, system (7.1) is null and approximately controllable in average for any probability distribution ρ.

• Under the hypotheses of Theorem 2.6, system (7.1) is null controllable in average and there are C, T 0 > 0 such that for all T ∈ (0, T 0 ] we have the bound:

K(G, G 0 , ρ, T ) ≤ Ce CT (2r-1) -1 .
As shown in [29, Section A.2], the controls can be obtained (when the system is controlable) by minimizing a quadratic functional, so they can be obtained by running simulations.

Simulations

In this section we illustrate experimentally the controllability results obtained in Corollary 7.4. Our objective is not to develop rigorous numerical methods for obtaining the controls, which goes beyond the objective of this work, but to get an insight on the differences between density functions inside and outside the nullcontrollability regime. For that, we recall that the optimal control for null controllability in average (when (7.1) is null controllable) is given by ϕ(t, x; φ)1 G0 , for ϕ the solution of (1.3) and φ the minimizer of: Figure 4: The state at time t = 1 of the averaged solutions of the heat equation after applying the controls constructed with the minimizer of J in V 20 , V 50 and V 100 with y 0 = 1 2 . In the left figure we have considered ρ = 1 (1,2) and we observe that those functions are controls that take the solutions almost to equilibrium. The computations are obtained by approximating the density functions with Dirac masses, which explains why the solution does not converge exactly to 0. In the right one, ρ = 1 (0,1) and we see that the minimizer of J do not take the solutions to equilibrium, so the functions that we have obtained cannot be considered controls.

J(φ) = 1 2 T 0 G0 +∞ 0 ϕ(t, x; α, φ)ρ(α)dα 2 dxdt + y 0 , +∞ 0 ϕ(0; α, φ)ρ(α)dα .
thus, most often we do not get the same trajectory, as:

t 0 +∞ 0 f α (s)e -α(t-s) ρ(α)dαds = t 0 +∞ 0 f α (s)ρ(α)dα +∞ 0
e -α(t-s) ρ(α)dα ds. A specific counter-example can be given with scalar ODEs. Let us consider a dynamic that behaves half of the times like:

u + u = f, and half of the times like:

v -v = g; that is, P[α = 1] = P[α = -1] = 1/2.
If the initial value is 1 we can take the solutions to rest at time T = 1 with the controls f (t) = -t and g(t) = t -2 (this is done by considering that 1 -t is a valid trajectory). The average of the controls is -1. However, for f (t) = g(t) = -t+t-2

2

= -1, we obtain the trajectories u(t) = -1 + 2e -t and v(t) = 1, whose average is e -t , which is not null at T = 1.

• Random initial data. As in some previous works involving average controllability results (see [START_REF] Zuazua | Averaged control[END_REF], [START_REF] Lü | Averaged controllability for random evolution Partial Differential Equations[END_REF], [START_REF] Coulson | On average controllability of random heat equations with arbitrarily distributed diffusivity[END_REF] and [START_REF] Lü | A Concise Introduction to Control Theory for Stochastic Partial Differential Equations[END_REF]) we fix the initial value. However, we can prove by linearity that averaged approximate controllability with fixed initial data implies averaged approximate controllability with random initial data

y 0 α if +∞ 0 y 0 α 2 L 2 (G) ρ(α)dα < ∞.
Nonetheless, in the case of null controllability in average to argue with linearity we would need to prove that the range of the average with random initial data is the same as the range with fixed initial value, which is an open problem. Alternatively, applying duality as in Theorem 7.3 could be a possibility. Indeed, we have to prove that If the initial value is independent of α, then +∞ 0

ϕ T → +∞ 0 y 0 α , ϕ(0, •; α) L 2 (G) ρ(α)dα
y α 0 , ϕ(0, •; α) L 2 (G) ρ(α)dα = y 0 , +∞ 0 ϕ(0, •; α)ρ(α)dα L 2 (G)
, so the equivalent observability inequality is

φ(0, •) L 2 (G) ≤ C φ 2 L 2 ((0,T )×G0) ,
but if the initial value depends on α, we need to prove that:

+∞ 0 ϕ(0, x; α, ϕ T ) 2 L 2 (G) ρ(α)dα ≤ C φ 2 L 2 ((0,T )×G0) . (8.1)
The main obstacle to prove (8.1) by replicating the proof of Theorem 2.6 is to prove an analogous result for (6.5). In recent papers where random initial values are considered the difficulty of a random initial value is bypassed by using exact averaged controllability (see, for instance, [START_REF] Lazar | Control of parameter dependent systems[END_REF]), which is satisfied by finite dimensional or hyperbolic systems but not by parabolic ones, or by assuming that there is no randomness in the dynamics, just on the initial value (see, for instance, [START_REF] Privat | Optimal shape and location of sensors for parabolic equations with random initial data[END_REF]). We highlight that understanding the controllability properties of (7.1) when y 0 α is a random initial value is an open problem whose resolution would help to have a more complete picture.

• Acting on the boundary. We may prove analogous controllability results to Theorems 2.3, 2. • More general random variables. Even if all the results in this paper have been stated for random variables with a density function, they are true for any random variable whose law satisfies the analogous inequalities of (2.3) and (2.4). Indeed, the proofs can be replicated step by step.

• More regular norms. Even if we have obtained all the results in this paper for the final state in L 2 (G) and we have made the observation in L 2 ((0, T ) × G 0 ), analogous results are valid for final states in H s1 (G) and the observation in H s2 ((0, T ) × G 0 ) (for any s 1 , s 2 ∈ R + ) for a domain G sufficiently regular. Indeed, the proofs are very similar with the only difference of some polynomial factors of N or λ. We recall that Lemma 3.2 can be adapted to observe a higher norm with the L 2 -norm. In fact, for any function φ = i∈Λ λ a i e i we have:

φ H 1 0 (G) = i∈Λ λ a 2 i λ i 1/2 ≤ √ λCe C √ λ i∈Λ λ a i e i L 2 (G0) ≤ Ce (C+1) √ λ φ L 2 (G0) .
• Analyticity of the space variable. If (α, G) satisfies the hypotheses of Theorem 2.6, we can easily prove as in [13, Theorem 1] that the free averaged solutions of the heat equation preserve the analyticity with respect to the spatial variable in the interior of the domain.

• Diffusion with negative values. Regarding the cases where the diffusion takes strictly negative values we do not have null averaged controllability. Indeed, under that hypothesis we can easily prove that:

lim n→∞ ũ(T, •; e n ) L 2 (G) ũ(•; e n ) L 2 ((0,T )×G)
= +∞.

• Some density functions in which the problem remains open. There are some density functions which satisfy neither (2.3) for some r > 1/2 nor (2.4) for some r < 1/2. For those density functions their (non-)observability properties are still unproved, for instance, those satisfying (1.6). It is an infinite dimensional class since it contains all functions provided by (2.7) for k = 1.

• Unique continuation. It would be interesting to have a proof of the unique continuation property for the averaged dynamics of any random variable α, even when it takes negative values. Indeed, there are some random variables whose density functions do not satisfy (2.3) for any r > 0 (for instance, ρ(α) = 2α1 (0,1) (α)), so their unique continuation is still unproved. In particular, we wonder if the unique continuation is preserved when ρ is too irregular, as a counter-example would probably be of such type.

• Measurable control domains. The observability properties proved in Theorem 2.6 can be extended to sets of the type E × G 0 , for E a measurable set. Indeed, we can use the approach of [START_REF] Phung | An observability estimate for parabolic equations from a measurable set in time and its applications[END_REF][START_REF] Apraiz | Observability inequalities and measurable sets[END_REF][START_REF] Lü | Averaged controllability for random evolution Partial Differential Equations[END_REF][START_REF] Coulson | On average controllability of random heat equations with arbitrarily distributed diffusivity[END_REF], which complement the ideas of [START_REF] Miller | A direct Lebeau-Robbiano strategy for the observability of heat-like semigroups[END_REF] with some results from Measure Theory.

• More general heat equation. An interesting problem that remains open is the study of the averaged observability properties of the random heat equation when the lower terms are also random terms, as:

y t -div(σ(x, α)∇y) + A(x, α) • ∇y + a(x, α)y = 0.
In particular, this is interesting when the averaged convection operator and the averaged diffusion operator do not commute. Unfortunately, the techniques presented in this paper do not help in that direction since they rely on being the eigenfunctions associated to the elliptic operator independent of α to ensure that the averages of the respective eigenfunctions remain orthogonal. Consequently, they can only be applied to equations of the type y t -αLy, for L an elliptic self-adjoint operator of compact resolvent. Thus, a theoretical or an in-depth numerical study would be of high interest for those operators.

• Numerics and simulations. It is an interesting question to develop rigorously the numerics for (7.1), including the speed of convergence. Determining precise methods that reproduce accurately the optimal control and the equation is also an open problem. Moreover, it would be interesting to illustrate with high accuracy that when ρ = 1 (0,1) the norms of the optimal controls for taking the average to a distance of at most ε of the origin explodes. In addition, it would be interesting to perform similar simulations with ρ = e -α 1 (0,+∞) and ρ = √ 21 (0,1) with the objective of determining numerically if with those density functions the approximate controllability in average holds.

• Other random equations. There are many other interesting questions involving random PDEs such as Schrödinger, wave, beam or Stokes equations:

-The Schrödinger equations with random diffusions satisfying the uniform, exponential, Laplace, normal, Chi-squared and Cauchy distributions were studied in [START_REF] Lü | Averaged controllability for random evolution Partial Differential Equations[END_REF]. There, the authors show that the averaged dynamics may be conservative or diffusive depending on the probability density, which leads to averaged controllability properties of very different kind. They consider the uniform distribution in any segment of R and the exponential distribution in [1, +∞), though their proof is valid in any segment of the type [K, +∞) for any K ∈ R. However, the problem of determining the dynamics and controllability properties of the averaged Schrödinger equations with arbitrary distributions is still open.

-The wave equation with random discrete diffusion was studied in [START_REF] Lazar | Averaged control and observation of parameter-depending wave equations[END_REF] and more abstractly in [START_REF] Lohéac | Averaged controllability of parameter dependent conservative semigroups[END_REF],

whereas understanding the general case is still an open challenge.

-Another interesting equation for which we can consider randomness in the higher order term is the Beam equation. In fact, in [START_REF] Marín | Robust averaged control of vibrations for the Bernoulli-Euler beam equation[END_REF] an optimization problem involving the cost of the control and the avarage of the square of the mass at a final time T is studied, but it is an open question to determine its exact (or null) controllability.

-The Stokes equation with random diffusion has not been studied in the literature. However, we can get analogous versions of Theorems 2. Let us fix q = (q 1 , . . . , q d ) and > 0 such that:

K := [q 1 , q 1 + ] × • • • × [q d , q d + ] ⊂⊂ G \ G 0 .
We also fix a positive non-trivial function ς ∈ D(B R d (0, 1)). We define for (γ 1 , . . . , γ d ) ∈ [0, 1] d : p(γ 1 , . . . , γ d ) := q + (γ 1 , . . . , γ d ),

which is a parametrization of K. With this in mind, we define the functions: • In order to have (5.2) we just need to find a non-trivial solution of the system: { φ N , e i L 2 (G) = 0, ∀i ∈ Λ N .

(A.3)

We remark that the system (A.3) is a linear homogeneous system with C√ N d unknowns (the constants (c i,N ) i ) and |Λ N | equations, so from Weyl's law (see Lemma 3.1) and by taking C large enough we obtain that there are more unknowns than equations, which implies that (A.3) has a non-trivial solution. In particular, we can fix (c i,N ) i a non-null tuple such that φ N is a solution of (A.3).

• In order to prove (5.3) it suffices to prove that for C > 0 large enough and all N ∈ N we have that:

∆φ N L 2 (G) ≤ CN 2 φ N L 2 (G) . (A.4)
Indeed, from (A.4) we obtain that:

φ N L 2 (G) ≥ 2 ∆φ N L 2 (G) CN ≥ 2 P ⊥ CN ∆φ N L 2 (G) CN ≥ 2 P ⊥ CN φ N L 2 (G) ,
so we find that: Consequently, the sequence φ N satisfies (A.4) for C large enough depending on C, and hence it also satisfies (5.3).

P CN φ N 2 L 2 (G) = φ N 2 L 2 (G) -P ⊥ CN φ N 2 L 2 (G) ≥ 3

0 e

 0 ) that: ũ(t, •; φ) = +∞ 0 e -αλ0t ρ(α)dα φ, e 0 L 2 (G) e 0 + i∈N * +∞ 0 e -αλit ρ(α)dα φ, e i L 2 (G) e i = +∞ -αλ0t ρ(α)dα φ, e 0 L 2 (G) e 0 + O e -(λ r 1

( a )

 a The control constructed with the minimizer of J in V20 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0, 1]. (b) The control constructed with the minimizer of J in V20 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0, 0.5]. (c) The control constructed with the minimizer of J in V20 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0.9, 1]. (d) The control constructed with the minimizer of J in V50 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0, 1]. (e) The control constructed with the minimizer of J in V50 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0, 0.5]. (f) The control constructed with the minimizer of J in V50 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0.9, 1]. (g) The control constructed with the minimizer of J in V100 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0, 1]. (h) The control constructed with the minimizer of J in V100 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0, 0.5].(i) The control constructed with the minimizer of J in V100 for y 0 = 1 2 , ρ = 1 (1,2) and t ∈ [0.9, 1].

Figure 2 :

 2 Figure 2: Graphs of the controls for ρ = 1 (1,2) and y 0 = 1 2 constructed with the minimizer of the functional J in V 20 , V 50 and V 100 . In the left column we illustrate the whole graphs, in the middle column we illustrate the graphs with the time variable zoomed in [0, 1/2], and in the left column zoomed in [0.9, 1].

4 and 2. 6 when•

 6 the control acts on the boundary. This can be done by repeating the proofs almost step by step, with the only difference of using [2, Theorem 9] instead of Lemma 3.2. We also remark that we cannot use an extension reduction technique as in [3, Section 3.3] since the trace on the boundary would depend on the parameter, and by the same argument of the first remark in Section 8, averaging the traces does not suffice. Neumann boundary conditions. We have analogous results of Theorems 2.3 and 2.6 for the controllability of the averaged solutions of the heat equation with random diffusion and Neumann boundary conditions. Indeed, we can repeat the proof step by step of those theorems since (3.3) is also true for Neumann boundary conditions (see [27, Theorem 2]). However, whether the analogous of Theorem 2.4 is true remains an open question since we do not have an analogous result of Lemma 3.4 for Neumann boundary conditions.

  3 and 2.6 for the Stokes equation with random diffusion as of the heat equation by considering [9, Theorem 3.1]. Nonetheless, determining if the analogous of Theorem 2.4 is true remains an open problem because a lack of a comparison theorem prevent from using analogous arguments.A Proof of Proposition 5.1As in Section 5, C denotes a sufficiently large positive constant that may be different each time it appears and which just depends on G, G 0 and ρ. In particular, it does not depend on the index N . Similarly, C is a constant sufficiently large that just depends on G, G 0 and ρ and C is a sufficiently large constant depending on those parameters and C. Finally, • denotes the floor function of a real number.

φ• 2 )

 2 N (x) := C√ N i1,...,i d =0 c i,N ς i,N (x), for ς i,N (x) := ς i := (i 1 , . . . , i d ) and c i,N and C a large constant to be defined later on (see Figure 5 for an illustration of how K and the support of φ 1 may look like). Let us check that for some C and c i,N the sequence φ N given by (A.1) satisfies (5.1)-(5.3): We have that: supp(φ N ) ⊂ x : d(x, Since the right-hand side of (A.2) is a decreasing sequence of sets and since K ⊂⊂ G \ G 0 we can easily prove (5.1) for C large enough.

Figure 5 :

 5 Figure 5: An illustration of the support of φ 1 in a domain belonging to R 2 .

4 φ N 2 L 2 ( 2 L 2 (= CN 2 φ N 2 L 2

 222222 G) , which is (5.3) squared. So, let us prove (A.4). We clearly have for all i, ĩ ∈ 0, . . . , C√ N d satisfying i = ĩ that supp(ς i,N ) ∩ supp(ς ĩ,N ) = ∅. Thus, we have that:∆φ N (G) . (A.5)

  By Proposition 4.1 and Lemma 4.2 we have that ũ(t, x; φ) = 0 in (0, +∞) × G 0 . Let us show that the first frequency of φ is null by contradiction. If the first frequency is not null, we obtain from (3.1) and (3.17

Acknowledgements:

This work was started in a scientific visit of the first author to the Friedrich-Alexander-Universität Erlangen-Nürnberg, funded by the Research Grants -Short-Term Grants, 2020 (57507441) of the German Academic Exchange Service. This project has also received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement NO: 694126-DyCon). The work of the first author was also supported by grants from Région Ile-de-France. The work of the second author is also partially supported by the Air Force Office of Scientific Research (AFOSR) under Award NO: FA9550-18-1-0242, by the Grant MTM2017-92996-C2-1-R COSNET of MINECO (Spain), by the Alexander von Humboldt-Professorship program, the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No.765579-ConFlex, and the Transregio 154 Project "Mathematical Modelling, Simulation and Optimization Using the Example of Gas Networks" of the German DFG.

M = 20, M = 50 and M = 100. This is motivated by the results presented in [START_REF] Abdelli | Numerical approximation of the averaged controllability for the wave equation with unknown velocity of propagation[END_REF] for the random wave equation.

In fact, we may prove as in [1, Section 5], by using variational and duality arguments, that if the system is null controllable in average the controls constructed with the minimizer of J in V M converge to the optimal control in L 2 ((0, T ) × G 0 ). We say that a control is constructed with the minimizer of J in V M if it is the restriction to (0, T ) × G 0 of the solution of (1.3) with ϕ T replaced by the minimizer. Since V M is a finite dimensional space, computing the minimizer of J is equivalent to solving numerically a linear system, which can be easily done by using any numerical computing environment (in our case MATLAB). We have obtained the following illustrations:

• We illustrate in Figure 2 (resp. in Figure 3) the controls constructed with the minimizer of J restricted to the spaces V M for ρ = 1 (1,2) (resp. for ρ = 1 (0,1) ). For ρ = 1 (1,2) we obtain a sequence of functions which seems to converge with respect to M to some function, which is something that can be seen in an even more clear way when t ∈ [0, 1/2]. Of course, the closer the time is to 1, the more slowly the punctual values of the control converge pointwise with M (and in t = 1 it diverges), but this is a well-known behavior when controlling a parabolic dynamics (see, for instance, [START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems[END_REF][START_REF] Münch | Numerical approximation of null controls for the heat equation: ill-posedness and remedies[END_REF][START_REF] Fernández-Cara | Strong convergent approximations of null controls for the 1D heat equation[END_REF]). However, for ρ = 1 (0,1) the sequence of functions does not seem to converge at all, which is something that we can appreciate in a more detailed way when t ∈ [0, 1/2]. What explains the difference is the controllability properties of each density.

• We illustrate in Figure 4 the state at t = 1 of the respective solutions of the averaged heat equation with the previously obtained functions. For ρ = 1 (1,2) we see that the state is taken slowly to 0, whereas for ρ = 1 (0,1) the picture shows some oscillations. This matches the theoretical results obtained in Corollary 7.4 as for ρ = 1 (1,2) the system is known to be averaged null controllable, whereas for ρ = 1 (0,1) the system is known to be not null averaged controllable. The computations are obtained by approximating the density functions with Dirac masses, which explains why the solution in the left figure does not converge exactly to 0.

Further comments and open problems

In this section we underline some extensions of our results to analogous situations and comment some interesting open problems:

• Average of the controls. A naive and incorrect way of computing a control that takes the average to rest is to compute the average of the controls that take each of the instances to rest. However, we do not get the same trajectory if we consider some source terms f α and then average on the solutions or if we compute the solutions with source term ∞ 0 f α ρ(α)dα and then average. In fact, the solutions of the equation:

u + αu = f α , are given by:


Davide Pierangeli

Mushegh Rafayelyan

Claudio Conti

Sylvain Gigan

Scalable spin-glass optical simulator

Many developments in science and engineering depend on tackling complex optimizations on large scales. The challenge motivates intense search for specific computing hardware that takes advantage from quantum features, nonlinear dynamics, or photonics. A paradigmatic optimization problem is finding low-energy states in classical spin systems with fully-random interactions. To date no alternative computing platform can address such spin-glass problems on a large scale. Here we propose and realize an optical scalable spin-glass simulator based on spatial light modulation and multiple light scattering. By tailoring optical transmission through a disordered medium, we optically accelerate the computation of the ground state of large spin networks with all-to-all random couplings. Scaling of the operation time with the problem size demonstrates optical advantage over conventional computing. Our results point out optical vector-matrix multiplication as a tool for spin-glass problems and provide a general route towards large-scale computing that exploits speed, parallelism and coherence of light.

INTRODUCTION

Non-deterministic polynomial-time (NP) problems are crucial from biochemistry to quantum physics. Their solution using polynomial resources requires nondeterministic Turing machines, which are unconventional computing models where a defined state can result in different outcomes [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. Alternative computing architectures exploit quantum annealing [START_REF] Johnson | Quantum annealing with manufactured spins[END_REF][START_REF] Boixo | Evidence for quantum annealing with more than one hundred qubits[END_REF], stochastic elements [START_REF] Borders | Integer factorization using stochastic magnetic tunnel junctions[END_REF], nonlinear dynamics with gain and losses [START_REF] Marandi | Network of time-multiplexed optical parametric oscillators as coherent Ising machine[END_REF][START_REF] Mcmahon | A fully-programmable 100spin coherent Ising machine with all-to-all connections[END_REF][START_REF] Inagaki | A coherent Ising machine for 2000node optimization problems[END_REF][START_REF] Böhm | A poor man's coherent Ising machine based on opto-electronic feedback systems for solving optimization problems[END_REF][START_REF] Inagaki | Large-scale Ising spin network based on degenerate optical parametric oscillators[END_REF][START_REF] Takesue | Simulating Ising Spins in External Magnetic Fields with a Network of Degenerate Optical Parametric Oscillators[END_REF], in-memory operations [START_REF] Traversa | Memcomputing NP-complete problems in polynomial time using polynomial resources and collective states[END_REF][START_REF] Cai | Harnessing Intrinsic Noise in Memristor Hopfield Neural Networks for Combinatorial Optimization[END_REF], or photon's speed and coherence [START_REF] Pierangeli | Large-Scale Photonic Ising Machine by Spatial Light Modulation[END_REF][START_REF] Prabhu | Accelerating recurrent Ising machines in photonic integrated circuits[END_REF][START_REF] Wu | An optical fiber network oracle for NPcomplete problems[END_REF][START_REF] Xu | A scalable photonic computer solving the subset sum problem[END_REF][START_REF] Farhat | Optical implementation of the Hopfield model[END_REF][START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF][START_REF] Estakhri | Inverse-designed metastructures that solve equations[END_REF][START_REF] Bueno | Reinforcement learning in a large-scale photonic recurrent neural network[END_REF][START_REF] Hamerly | Large-Scale Optical Neural Networks Based on Photoelectric Multiplication[END_REF][START_REF] Marcucci | Theory of neuromorphic computing by waves: machine learning by rogue waves, dispersive shocks, and solitons[END_REF][START_REF] Rafayelyan | Large-Scale Optical Reservoir Computing for Spatiotemporal Chaotic Systems Prediction[END_REF]. Among them, Ising machines are special-purpose processors designed for finding a ground state of a Ising spin model. They are currently attracting broad attention, since tasks such as partitioning, routing, and encrypting can be mapped on Ising Hamiltonians [START_REF] Lucas | Ising formulations of many NP problems[END_REF]. Devices based on various physical mechanisms have been recently realized using superconducting networks [START_REF] Harris | Phase transitions in a programmable quantum spin glass simulator[END_REF][START_REF] Hamerly | Experimental investigation of performance differences between coherent Ising machines and a quantum annealer[END_REF], optical parametric oscillators [START_REF] Mcmahon | A fully-programmable 100spin coherent Ising machine with all-to-all connections[END_REF][START_REF] Inagaki | A coherent Ising machine for 2000node optimization problems[END_REF][START_REF] Böhm | A poor man's coherent Ising machine based on opto-electronic feedback systems for solving optimization problems[END_REF], polariton condensates [START_REF] Berloff | Realizing the classical XY Hamiltonian in polariton sim-ulators[END_REF][START_REF] Kalinin | Polaritonic XY-Ising Machine[END_REF], coupled laser cavities [START_REF] Parto | Realizing spin Hamiltonians in nanoscale active photonic lattices[END_REF][START_REF] Pal | Rapid fair sampling of the XY spin Hamiltonian with a laser simulator[END_REF], nanophotonic circuits [START_REF] Prabhu | Accelerating recurrent Ising machines in photonic integrated circuits[END_REF][START_REF] Okawachi | Demonstration of chipbased coupled degenerate optical parametric oscillators for realizing a nanophotonic spin-glass[END_REF] and spatial light modulators (SLM) [START_REF] Pierangeli | Large-Scale Photonic Ising Machine by Spatial Light Modulation[END_REF][START_REF] Pierangeli | Noise-enhanced spatial-photonic Ising machine[END_REF][START_REF] Kumar | Large-scale Ising emulation with four body interaction and all-to-all connections[END_REF][START_REF] Fang | Experimental Realization of Gauge Transformation and Observation of Phase Transition in Spatial Optical Spin Systems[END_REF][START_REF] Pierangeli | Adiabatic evolution on a spatial-photonic Ising machine[END_REF]. Scalability with respect to the problem size is the main factor hampering their nearterm application. In fact, several Ising machines, such as the D-Wave quantum annealer [START_REF] Hamerly | Experimental investigation of performance differences between coherent Ising machines and a quantum annealer[END_REF], rely on local interactions between their elementary units, a fact that strongly limits long-range connectivity and imposes redundant schemes difficult to scale in practice [START_REF] Lechner | A quantum annealing architecture with all-to-all connectivity from local interactions[END_REF]. Other platforms, such as coherent Ising machines (CIMs) [START_REF] Mcmahon | A fully-programmable 100spin coherent Ising machine with all-to-all connections[END_REF][START_REF] Inagaki | A coherent Ising machine for 2000node optimization problems[END_REF][START_REF] Böhm | A poor man's coherent Ising machine based on opto-electronic feedback systems for solving optimization problems[END_REF][START_REF] Inagaki | Large-scale Ising spin network based on degenerate optical parametric oscillators[END_REF], provide all-to-all connectivity and can host dense spin networks made of thousands of elements, but with couplings that are not fully programmable and assume only a few possible values. For these reasons, the relevant NP problems that can be implemented and solved heuristi-cally on Ising machines on large scale are still not exhaustive.

In this Article, we report a pivotal step toward "Ising computing" by realizing a scalable photonic device that can simulate large-scale spin problems with continuous random couplings. We demonstrate use of optical random vector-matrix multiplications to implement the energy function of a spin-glass system. Since the optical setting enables simultaneous processing of all spin interactions in parallel, our approach exhibits an optical advantage at large scale over digital computing. The photonic hardware accelerates the solution of the spinglass problem independently of the used algorithm, which suggests that our setup may potentially speed-up any minimization approach. We apply the optical simulator to the number-partitioning problem, thus proving it can be useful for a vast class of practical combinatorial optimization tasks. Although we program only the coupling distribution, our scheme may develop into a fullyprogrammable special-purpose optical processor using reconfigurable transmissive elements [START_REF] Leedumrongwatthanakun | Programmable linear quantum networks with a multimode fibre[END_REF][START_REF] Matthès | Optical complex media as universal reconfigurable linear operators[END_REF][START_REF] Spall | Fully reconfigurable coherent optical vector-matrix multiplication[END_REF].

Finding the minimum energy configuration of a spin glass (SG) is a benchmark NP-hard problem [START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF][START_REF] Barahona | On the computational complexity of Ising spin glass models[END_REF], and its computer intractability continuously inspires novel heuristic algoritms [START_REF] Goto | Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems[END_REF][START_REF] Roques-Carmes | Heuristic recurrent algorithms for photonic Ising machines[END_REF]. The system can be illustrated as in Fig. 1(a), where a set of N unitary Ising spins σ i ∈ {+1, -1} occupies the sites of a disordered lattice. Due to strong lattice distortions, the effective interaction J ij between the i-th and j-th spin takes a broad spectrum of values. The quadratic SG Hamiltonian has the form

H(σ) = - 1 2 N i,j=1 J ij σ i σ j (1)
where the J ij elements come from a Gaussian distribution function P (J ij ). The model is a cornerstone of statistical mechanics also known as Sherrington-Kirkpatrick (SK) model [START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF]. Each problem instance corresponds to a graph of N all-to-all connected nodes with a set of randomly weighted links [Fig. 1(b)].

The operating principle of our optical SG simulator is shown in Fig. 1(c). The basic idea is to encode the spins on a coherent wavefront by spatial light modulation [START_REF] Pierangeli | Large-Scale Photonic Ising Machine by Spatial Light Modulation[END_REF] and their interaction on the optical transmission matrix (TM) of a disordered medium [START_REF] Matthès | Optical complex media as universal reconfigurable linear operators[END_REF][START_REF] Popoff | Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media[END_REF]. A similar approach has been recently investigated as an instruments to access spin-glass dynamics and its complexity [START_REF] Leonetti | Optical computation of a spin glass dynamics with tunable complexity[END_REF]. Specifically, we consider the optical field transmitted via multiple scattering E m = i t m i E i , where 1 < i < N and t m i is the complex TM element connecting the i-th input mode (spin) generated by an SLM to the m-th output mode detected by a camera [START_REF] Popoff | Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media[END_REF]. The total intensity transmitted over M output modes is thus

I T = m |E m | 2 = m i,j tm i t m j Ēi E j .
Defining the spin variables via the optical phase delays φ i ∈ {0, π}, so that σ i = exp(iφ i ) = E i , we obtain (see Appendix A)

I T = -H(σ); J ij = M m=1 J m ij = M m=1 Re( tm i t m j ). (2)
Equation ( 2) establishes a direct relation between the scattered intensity and the SG energy. Since the TM of a disordered medium is a random full-rank matrix [START_REF] Popoff | Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media[END_REF], when M = N we find that J ij has uncorrelated random elements with Gaussian distribution P (J ij ), as for the SG model in Eq. [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. Shaping the binary input phase distribution to maximize the transmitted intensity corresponds to looking for the SG ground state. Energy minimization can be performed with any iterative method while the spin system is optically emulated.

EXPERIMENTAL SETUP AND METHODS

The optical SG is numerically simulated by forming N pixel blocks from a square mesh (SLM plane). The initial optical field E I has constant amplitude and its phase is a random configuration of N binary phases, φ i = 0, π. A unitary TM matrix W with random complex numbers is generated. At each iteration, a single spin (phase value φ i ) is randomly selected and flipped; the optical field linearly propagates, E T = W • E I , and the input phase is updated only if the output total intensity I T increases. Numerical evaluation of I T corresponds to a measurement with a 64-bit sensitivity detector in a noiseless system (Fig. 2). In general, within this scheme, ∼ 10N iterations are sufficient for a good convergence. We normalize the transmitted intensity to the initial transmission, which allows us to compare the result with experiments at any input optical power.

The experimental device follows the setup illustrated in Fig. 3(a). A continuous-wave laser beam at λ = 532 nm is expanded, polarization controlled, and impinges on a reflective liquid-crystal SLM (Meadowlark Optics HSP192-532, 1920 × 1152 pixels) performing phase-only light modulation. The SLM area is divided into N addressable optical spins by grouping several pixels. Binary modulated light is projected on the objective backfocal plane (OBJ1, 10x, NA = 0.1) and it is focused on a strongly scattering medium (a thick diffuser made of teflon, DIFF) with 0.5mm thickness. Scattered light is collected by a second objective (OBJ2, 20x, NA = 0.4) and the transmitted intensity speckle pattern is detected by a CCD camera (Basler acA2040-55um, 2048 × 1536 pixels) with 8-bit (256 gray-levels) intensity sensitivity on each pixel. Each camera pixel has a size comparable with the spatial extent of a speckle grain and thus corresponds to an output spatial mode. The SLM and CCD have communication bandwidths of 4000MB/s and 600MB/s, respectively.

The ground-state search is conducted sequentially by means of the digital recurrent feedback. Computation starts from a random configuration of N binary phase blocks (spins) on the SLM. The measured intensity distribution determines the feedback signal. The SNR is approximatively 50, with the main noise sources that are associated to flickering effects of the phase modulator, power laser fluctuations, and tiny mechanical effects on the scattering medium. At each machine cycle a batch of spins is randomly selected and flipped; the intensity transmitted on M camera pixels is detected and the spin state is updated if the change increases I T . The batch size is selected as 2.5% of the spin number, which ensures that a single change on the SLM is detected over camera noise. This values determines the maximum achievable accuracy for the experimental ground state.

RESULTS

The scheme of the optical SG simulator is numerically validated in Fig. 2, where we model a large-scale device with N = M =1024. Linear optical propagation through the scattering medium is simulated by the randomly generated TM in Fig. 2(a). According to Eq. ( 2), the TM gives a J ij set following a Gaussian probability density with zero mean and deviation J = 1/4N (Fig. 2

Supplementary Material [45]

). In Fig. 2(c) we show the total intensity, normalized to the initial transmission, that is transmitted on the output modes during the optimization procedure. While I T increases and saturates to a final speckle distribution [Fig. 2(e)], the binary phases on the SLM converge towards a state minimizing the corresponding spin energy [Fig. 2(d)]. The final intensity corresponds to the ground-state energy [Eq. ( 2)], apart from a constant factor. To demonstrate that the method solves the spin problem, we benchmark the ground-state energies with simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] (SA) at zero temperature on the same random graph. Results on 100 independent runs are in Fig. 2(f) and indicate that our model operates with an accuracy comparable with a standard robust optimization algorithm. States with lower energy can be found by refining the iterative method by introducing effective temperature variations both in SA and the optical model (Supplementary Fig. 2 [45]). Simulations on various instances (Supplementary Fig. 1 [45]) further indicate the effectiveness of our approach in finding heuristic SG solutions. We find that the low-energy states obtained for different J ij realizations overlaps with the solutions found varying only the initial condition. Therefore, the SG low-energy space can be sampled by replicating either over the initial spin configuration or realization of the random couplings.

Experimental implementation

We realize the optical SG simulator according to the experimental setup in Fig. 3(a). Ising spins are encoded on a laser wavefront by a phase-only SLM, a volumetric diffuser provides multiple scattering, and camera pixels are the output modes. The optical device works in a measurement and feedback scheme. At each machine iteration, we measure the intensity I T on M = N fixed camera pixels and update the spins in order to maximize the transmission (see details in Appendix B). The setup can make use of any optimization scheme in this operation, i.e., it is algorithm-agnostic. After a few thousands of iterations, we get a transmission enhancement (normalized transmitted intensity) close to the expected value, with variations depending on the the random input condition [inset in Fig. 3(a), N =256]. In analogy with the numerical findings, from the measured final intensity we get the ground-state energy for each realization.

Using the optical setting, we performed sets of compu- tations for random SG problems of different sizes, up to more than 10 4 spins and 10 8 connections. We quantify the solutions found by analyzing their SG energy in comparison with numerical models. Main results are summarized in Fig. 3(b-c). For N =256, the Ising machine finds an approximate solution to the NP-hard problem with an accuracy comparable with the optical SG model and SA. Optical computing was successful also for largescale graphs with 1024 nodes, although with lower performance. For a fixed number of machine iterations (Appendix B), we find the ground-state accuracy decreases as a function of the system size [Fig. 3(d)]. On the contrary, numerical results indicate that the optical SG simulator as modelled in Fig. 2 is able to perform independently of the system size. The observed behavior is thus a direct consequence of the experimental conditions, i.e., practical non-idealities of the device. It is important to remark that a similar scaling effect occurs in any other quantum and classical optimizers built in practice, even at much smaller sizes [START_REF] Hamerly | Experimental investigation of performance differences between coherent Ising machines and a quantum annealer[END_REF].

To understand how to improve the computational ability of our proof-of-principle device when the number of spins increases, we consider the effect of various experimental factors on the optical SG model. Specifically, we analyze the dependence of the ground-state accuracy on the detection noise level and the impact of the finiteprecision of the camera at various sizes. Figure 3(e) shows the results when varying the signal to noise ratio (SNR) for a N =1024. A rapid decrease of the solution accuracy is observed as the noise level exceeds 10 -2 , which indicates that optical noise partially explains the experimental performance. A key role is also played by the finite-precision of the camera. The accuracy improves considerably as we increase the detector bit precision, even in presence of noise [Fig. 3(e)]. The effect of the components precision on the computation is even more crucial as the problem size increases (see Supplementary Fig. 4 [45]). Additional evidences in Supplementary Material [45] further indicate that the performance scaling in the realized simulator [Fig. 3(d)] is successfully modelled. The overall analysis suggest that, by concomitantly reducing noise and increasing the precision of the optical readout, the optical simulator can be effective also on large scales.

From the statistical physics point of view, the operation of the optical SG simulator is limited by a finite effective temperature T * . We estimate it numerically by using an inverse numerical approach. Exploiting SA, we anneal the SG up to the measured ground-state energy, we let it equilibrate and extract its temperature. For N=4096 we obtain β * = 1/T * ≈ 0.9. Cooling to lower temperatures can be achieved improving the device construction to make it more sensitive and noise tolerant.

A crucial parameter of our spin-glass simulator is the number of output modes. In fact, according to Eq. ( 2), the transmitted intensity, and thus spin Hamiltonian, depends on the output mode number M . It determines the rank of the interaction matrix and, consequently, the correlations and the distribution of its values. As reported in Supplementary Fig. 3 [45], by increasing M the couplings J ij evolves from sharp peaked distribution to the Gaussian probability density P (J ij ). Ising problems in which M = N , although differ from the considered SK model, can also be particularly interesting for applications. For instance, the M = 1 case directly maps to the number-partitioning problem [START_REF] Mertens | Phase transition in the number partitioning problem[END_REF], which is the combinatorial optimization problem in which a set of real number must be divided into two subsets differing as little as possible in their weight. The problem is NP-complete and represents a typical task encountered in resource allocation [START_REF] Lucas | Ising formulations of many NP problems[END_REF]. For example, in cases where we have to divide a set of assets fairly between two people. Application of the experimental SG simulator to number partitioning is demonstrated in Supplementary Fig. 5 [45]. The partition solution is optically found with good accuracy. Remarkably, in this case the performance does not degrade with the problem size [START_REF] Fang | Experimental Realization of Gauge Transformation and Observation of Phase Transition in Spatial Optical Spin Systems[END_REF], and the efficiency is maintained for sets exceeding 10 4 random numbers. This underlines that our optical device can be directly applied to specific computing tasks, and can be beneficial in a broad range of applications. Among these, we mention finding cliques in networks, which is central for understanding social dynamics.

Optical advantage

The key advantage of our optical SG simulator is its possible scalability to sizes intractable with conventional hardware. In fact, common algorithms require to evaluate Eq. ( 1) at each iteration, an operation which time and memory consumption grows quadratically with the spin number. The optical part of our scheme executes such matrix multiplication fully in parallel, independently of the problem size and feedback algorithm. The scaling advantage is demonstrated in Fig. 4 by measuring the iteration time versus the problem size. In contrast to the quadratic scaling of the SG model on a conventional computer, the optical computation time scales only linearly, with a mild slope depending only on the limited communication bandwidth of the electronic feedback. Therefore, independently of the machine operation frequency, scaling laws ensure the existence of an optical advantage region at large scales. The sensitivity of light modulators and detectors, and more generally, optical noise, rules the maximum size that can be efficiently solved on optical platforms.

CONCLUSION

In conclusion, we have reported a scalable optical device able to solve random spin problems. Exploiting spatial light modulation and coherent optical propagation of light, our scheme allows parallel information processing for arbitrary problem sizes and without any fabrication constraints. Our setup can be exploited as an optical accelerator for the solution of spin glasses with any optimization algorithm [START_REF] Pierangeli | Adiabatic evolution on a spatial-photonic Ising machine[END_REF]. Given the 8-bit precision of the intensity detector and the effect of experimental noise, the analog spin simulator finds ground states with energies higher than those obtained with simulated annealing on a 64-bit digital processors. However, in principle, the same accuracy can be reached by tuning further the optical setup and the optimization procedure.

The use of a physical medium to encode spin interactions also opens interesting perspectives for programming arbitrary Ising problems, which could be done by selecting various subset of input and output modes [START_REF] Leedumrongwatthanakun | Programmable linear quantum networks with a multimode fibre[END_REF][START_REF] Matthès | Optical complex media as universal reconfigurable linear operators[END_REF], or by directly tailoring the transmission matrix using either microfabrication, or a second spatial light modulator [START_REF] Spall | Fully reconfigurable coherent optical vector-matrix multiplication[END_REF]. Our approach points out a parameters region where optical computing presents a favourable scaling with respect to electronic hardware. Developments in photonic technology would allow to optically address many NPhard combinatorial optimizations deep into this region, where neuromorphic computing can also find its natural application [START_REF] Bueno | Reinforcement learning in a large-scale photonic recurrent neural network[END_REF][START_REF] Hamerly | Large-Scale Optical Neural Networks Based on Photoelectric Multiplication[END_REF][START_REF] Rafayelyan | Large-Scale Optical Reservoir Computing for Spatiotemporal Chaotic Systems Prediction[END_REF]. The TM models monochromatic transmission through a linear optical system at the mesoscopic level. Its complex coefficients t m i connect amplitude and phase of the optical field between the m-th output mode and the ith input element, E m = N i t m i E i , where we adopt superscripted indices only for clarity. As shown in Ref. [START_REF] Popoff | Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media[END_REF], the elements t m i are uncorrelated random complex numbers when a thick disordered medium is placed between a SLM and a camera, and they can be measured experimentally. We first consider the intensity on a single output mode, which reads as

I m = | i t m i E i | 2 = tm j t m i Ēi E j .
Defining the spins via the binary phase delays φ i ∈ {0, π}, so that E i = exp(iφ i ) = σ i up to a global phase factor, we get the Ising Hamiltonian

I m = -H m = ij J m ij σ i σ j with J m ij = Re( tm j t m i )
, apart from constant factors. Pairs of spins with positive (negative) interaction correspond to points of the optical field resulting in constructive (destructive) interference. In this case, the couplings are correlated, i.e., rank(J m ij ) = 1 (the interaction is specified by only N degrees of freedom). This case corresponds to a class of Ising problems, known as Mattis SG, that have an exact groundstate solution [START_REF] Mattis | Solvable spin systems with random interactions[END_REF]. In optics such solution corresponds to the optimal wavefront shaping for focusing on a single output mode, which, in conditions of negligible noise, gives an enhanced trasmission proportional to N [START_REF] Rotter | Light fields in complex media: mesoscopic scattering meets wave control[END_REF]. In combinatorial optimization, such configuration directly maps to number partitioning problem [START_REF] Fang | Experimental Realization of Gauge Transformation and Observation of Phase Transition in Spatial Optical Spin Systems[END_REF]. Application of the optical SG simulator for finding specific partitions in a set of random numbers is detailed in Supplementary Material [45].

The interaction matrix and its probability distribution varies considerably when increasing the number of output channels (Supplementary Fig. 3 [45]). For M modes, we have

I T = m I m = m ( j tm j σ j )( i t m i σ i ) = m i,j tm i t m j E i E j
, which gives the equivalence in Eq. ( 2):

I T = -ij J ij σ i σ j with J ij = M m=1 Re( tm i t m j
). The coupling matrix rank is now rank(J ij ) = M . When M = N , we get a full-rank matrix (N 2 variables specify the couplings) describing random uncorrelated spin interactions. The J ij are distributed with a Gaussian density P (J ij ) (zero mean and deviation 1/4N ), as verified in Supplementay Material [45] for both numerical and experimental data (Supplementary Fig. 1). Simultaneous maximization of I T over N output modes is equivalent to minimizing the energy of a SG, with interactions encoded in the TM.

APPENDIX B: NUMERICAL AND EXPERIMENTAL DETAILS

During the optimization, the SG energy is evaluated by applying Eq. ( 2) on the optical phase distribution. At the ground state, such energy is related to the optimized transmitted intensity by a constant factor. This factor depends only on the spin number. We exploit this property in experiments, where the factor extracted from a set of measured TMs at a given size is used also on optical computing runs in which the TM is varied in each realization. The ground-state accuracy in Fig. 3

is defined as 1 -[(G -G min )/((G + G min ))],
where G and G min are the mean SG energy measured on the device (numerical implementation) and its computer model (SA), respectively for purple (green) dots. Error bars indicate one standard deviation over 20 realizations.

To characterize the device in experimental conditions, we introduce various ingredients in the model. Finiteprecision of the camera is obtained by discretizing the output intensity in 2 n levels, n being the number of bits. We introduce optical noise at the readout, by adding uncorrelated Gaussian fluctuations to the intensity I T , with the SNR that determines the fluctuations amplitude (noise level). All codes are implemented in MATLAB on an Intel processor with 6 cores running at 4.1 GHz and supported by 16 GB ram. In Fig. 4, iteration times for standard computing refer to this specific CPU. We note that high-performance computing on dedicated systems can substantially reduce these iteration times; however, their quadratic scaling will remain unaltered. As for simulated annealing (SA), a custom optimized version has been implemented following Ref. [START_REF] Isakov | Optimised simulated annealing for Ising spin glasses[END_REF]. The code exploits various methods including sequential updating, forward energy computation and fast pre-computed random numbers. It has been benchmarked on standard graphs, including K2000, with results analogous to Ref. [START_REF] Goto | Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems[END_REF] in terms of ground-state energies.

In our scheme, each SG graph corresponds to a measured TM with size N 2 [START_REF] Popoff | Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media[END_REF]. The TM is experimentally reconstructed using Non-negative Matrix Factorization and a phase retrieval algorithm [START_REF] Boniface | Non-invasive focusing and imaging in scattering media with a fluorescencebased transmission matrix[END_REF]. Slight translations/rotations of the disordered medium result in a different TM. The optical stability of the scattering medium (approximatively one hour) fixes the physical time for which the interaction matrix remains unaltered. This factor limits the optimization effectivness over long times. We thus kept the number of iterations in each run constant to 16000 and collect tens of computation varying only the initial condition. However, faster optical elements can considerably lower the total computation time of our optical SG simulator. The optical setup operates at 150Hz and, according to the employed SLM technology, the iteration time can be reduced up to 1.4 milliseconds, maintaining its linear dependence on the problem size.

FIG. 1 .

 1 FIG. 1. Scheme of the optical spin-glass simulator. (a) Sketch of Ising spins on a disordered lattice. (b) Graph representation of a SG problem for various sizes. Each spin is a node of a fully-connected network where the coupling matrix Jij is represented by color-coded links. (c) Optical scheme mapping the SG model. A spatial light modulator (SLM) inscribes Ising spins in separated spatial points of the optical field Ei. The spin network is encoded in a disordered medium that mixes all the incoming modes according to its scattering matrix. Any i-th and j-th two spins contribute to the m-th output mode by a coupling coefficient J m ij . The computation works by optimizing the total intensity transmitted on a set of output modes.

FIG. 2 .

 2 FIG. 2. Design and validation of the optical spin-glass scheme. (a) TM modeling multiple light scattering between N = 1024 optical spins and M = N detector modes. (b) Coupling's probability distribution with Gaussian fit (line). (c) Total transmitted intensity during the computation for different initial conditions and (d) energy of the corresponding SG. Insets in (c-d) show a ground-state phase configuration and its problem graph. (e) Final intensity with transmission maximized on the dashed region. (f) Energy histograms of ground states found using the optical SG model and zero-temperature SA.

  (b) and

FIG. 3 .

 3 FIG. 3. SG computing device. (a) Experimental setup, where recurrent feedback from the measured intensity updates the spin configuration on the SLM. Insets show a measured TM and corresponding transmitted intensity during computation for several runs. (b-c) SG ground-state energy histogram for N =256 and N =1024. Inset table indicates values for solutions computed optically and with computer algorithms. Vertical lines are average solutions from numerical simulations of the optical SG model. (d) Accuracy of the ground-state energy varying the spin number. Experimental results are compared to the scaling behavior of the ideal device. (e) Modelling of the experimental device: ground-state accuracy as a function of optical noise at the detection (SNR, signal to noise ratio) for intensity detectors with different finite precision. The shaded area indicates the parameter region in which the realized device (8-bit precision) is expected to be operating.

FIG. 4 .

 4 FIG.4. Optical advantage. Time for updating the SG configuration (iteration time) versus its size. Values refer to our scheme implemented on our optical setup (green dots) and on a standard digital computer (blue dots). Lines are fit functions showing linear and quadratic scaling of optical and electronic computing. Inset highlights the setup-dependent crossing point (N * = 1460) between the scaling behaviors, which delineate an optical advantage region for large scales (green-shaded area). The vertical dotted line, referred as noise limit, is an hypothetical border to indicate the maximum size solvable with an optical SG simulator given a finite optical noise level.

  APPENDIX A: SPIN-GLASS HAMILTONIAN IN THE TM FRAMEWORK

ACKNOWLEDGEMENTS D.P. and C.C. acknowledge funding from SAPIExcellence 2019 (SPIM project), QuantERA ERA-NET Cofund (Grant No. 731473, project QUOMPLEX), PRIN PELM 2017 and H2020 PhoQus project (Grant No. 820392). M.R. and S.G. acknowledge support from H2020 European Research Council (ERC) (Grant No. 724473), DARPA (HR00111890042) and Institut Universitaire de France. We thank J. Dong and A. Boniface for useful discussions and technical support in the laboratory. Author contributions. D.P. and M.R. carried out experiments, numerical simulations and data analysis. D.P wrote the paper with contributions from all the authors. C.C. and S.G. co-supervised the project.