Quantum simulation of thermodynamics in an integrated quantum photonic processor

One of the core questions of quantum physics is how to reconcile the unitary evolution of quantum states, which is information-preserving and time-reversible, with evolution following the second law of thermodynamics, which, in general, is neither. The resolution to this paradox is to recognize that global unitary evolution of a multi-partite quantum state causes the state of local subsystems to evolve towards maximum-entropy states. In this work, we experimentally demonstrate this effect in linear quantum optics by simultaneously showing the convergence of local quantum states to a generalized Gibbs ensemble constituting a maximum-entropy state under precisely controlled conditions, while introducing an efficient certification method to demonstrate that the state retains global purity. Our quantum states are manipulated by a programmable integrated quantum photonic processor, which simulates arbitrary non-interacting Hamiltonians, demonstrating the universality of this phenomenon. Our results show the potential of photonic devices for quantum simulations involving non-Gaussian states.


INTRODUCTION
One of the long-standing puzzles of theoretical physics is how notions of statistical physics and of basic quantum mechanics fit together in closed systems [1]. Statistical mechanics is concerned with probabilistic, stationary ensembles that maximize entropy under external constraints. Elementary quantum mechanics, in contrast, describes the deterministic evolution of quantum states of closed systems under a specified Hamiltonian. It has become clear [2][3][4][5] that these seemingly contradictory premises can be resolved by making the distinction between global unitary dynamics and local relaxation (see Fig. 1). The physical mechanism is that local expectation values converge to those of statistical ensembles, while the entire closed quantum system undergoes unitary dynamics. Large-scale, closed quantum systems therefore appear locally thermal without the need to postulate an external heat bath. Crucially, this local equilibration behaviour is believed to be ubiquitous, in the sense that one has to fine-tune the Hamiltonian in order to not observe it [5,6].
The mechanism of local equilibration is particularly clearcut under non-interacting quadratic bosonic Hamiltonians, such as describe linear quantum optics. If the initial state is non-Gaussian, it is expected to 'Gaussify' in time, i.e., to locally converge to Gaussian states that maximize the entropy given all second moments of the state [7][8][9][10][11]. In this case, < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 S z I g H V U G 1 k c E 6 8 2 A 6 q z w g l Z 8 Y w = " > A A A B + X i c d Z D L S s N A F I Y n 9 V b r L e r S z W A R 6 s K S t P a y E Y p u u q x g 2 k I b y 2 Q 6 a Y d O J m F m U i i h b + L G h S J u f R N 3 v o 2 T t o K K / j D w 8 Z 9 z O G d + L 2 J U K s v 6 M D J r 6 x u b W 9 n t 3 M 7 u 3 v 6 B e X j U l m E s M H F w y E L R 9 Z A k j H L i K K o Y 6 U a C o M B j p O N N b t J 6 Z 0 q E p C G / U 7 O I u A E a c e p T j J S 2 B q b p F N r n 8 A q S + + S C N q G a D 8 y 8 V b T s U r 1 S h i l Y l W o l B a t U r d n Q 1 p A q D 1 Z q D c z 3 / j D E c U C 4 w g x J 2 b O t S L k J E o p i R u a 5 f i x J h P A E j U h P I 0 c B k W 6 y u H w O z 7 Q z h H 4 o 9 O M K L t z v E w k K p J w F n u 4 M k B r L 3 7 X U / K v W i 5 V f d x P K o 1 g R j p e L / J h B F c I 0 B j i k g m D F Z h o Q F l T f C v E Y C Y S V D i u n Q / j 6 K f w f 2 q W i X S 2 W b y / z j e t V H F l w A k 5 B A d i g B h q g C V r A A R h M w Q N 4 A s 9 G Y j w a L 8 b r s j V j r G a O w Q 8 Z b 5 9 L U J I t < / l a t e x i t > U (V ) = e iHt < l a t e x i t s h a 1 _ b a s e 6 4 = " B f / L c T k 7 + A u D K y v s R D k G I / 8 q 6 L 0 = " > A A A B 8 H i c d Z D L S g M x F I Y z X m u 9 V V 2 6 C R b B V Z l p 7 W V Z d O O y g r 1 I O 5 R M e t q G J p k h y Q h l 7 F O 4 c a G I W x / H n W 9 j p q 2 g o j 8 E P v 5 z D j n n D y L O t H H d D 2 d l d W 1 9 Y z O z l d 3 e 2 d 3 b z x 0 c t n Q Y K w p N G v J Q d Q K i g T M J T c M M h 0 6 k g I i A Q z u Y X K b 1 9 h 0 o z U J 5 Y 6 Y R + I K M J B s y S o y 1 b u + 9 n i J y x K G f y 7 s F 1 y v W y i W c g l u u l F N w i 5 W q h z 0 L q f J o q U Y / 9 9 4 b h D Q W I A 3 l R O u u 5 0 b G T 4 g y j H K Y Z X u x h o j Q C R l B 1 6 I k A r S f z B e e 4 V P r D P A w V P Z J g + f u 9 4 m E C K 2 n I r C d g p i x / l 1 L z b 9 q 3 d g M a 3 7 C Z B Q b k H T x 0 T D m 2 I Q 4 v R 4 P m A J q + N Q C o Y r Z X T E d E 0 W o s R l l b Q h f l + L / o V U s e J V C 6 f o 8 X 7 9 Y x p F B x + g E n S E P V V E d X a E G a i K K B H p A T + j Z U c 6 j 8 + K 8 L l p X n O X M E f o h 5 + 0 T B I a Q k g = = < / l a t e x i t > |1i < l a t e x i t s h a 1 _ b a s e 6 4 = " c / X F 3 O W N t 3 I 3 g t X V g 1 l 4 I E E j 6 o o = " > A A A B 8 H i c d Z D L S g M x F I b P e K 3 1 V n X p J l g E V 2 W m t Z d l 0 Y 3 L C v Y i 7 V A y a d q G J p k h y Q h l 7 F O 4 c a G I W x / H n W 9 j p q 2 g o j 8 E P v 5 z D j n n D y L O t H H d D 2 d l d W 1 9 Y z O z l d 3 e 2 d 3 b z x 0 c t n Q Y K 0 K b J O S h 6 g R Y U 8

k b R p m O O 1 E i m I R c N o O J p d p v X 1 H l W a h v D H T i P o C j y Q b M o K N
t W 7 v 3 Z 7 C c s R p P 5 d 3 C 6 5 X r J V L K A W 3 X C m n 4 B Y r V Q 9 5 F l L l Y a l G P / f e G 4 Q k F l Q a w r H W X c + N j J 9 g Z R j h d J b t x Z p G m E z w i H Y t S i y o 9 p P 5 w j N 0 a p 0 B G o b K P m n Q 3 P 0 + k W C h 9 V Q E t l N g M 9 a / a 6 n 5 V 6 0 b m 2 H N T 5 i M Y k M l W X w 0 j D k y I U q v R w O m K D F 8 a g E T x e y u i I y x w s T Y j L I 2 h K 9 L 0 f / Q K h a 8 S q F 0 f Z 6 v X y z j y M A x n M A Z e F C F O l x B A 5 p A Q M A D P M G z o 5 x H 5 8 V 5 X b S u O M u Z I / g h 5 + 0 T A v u Q k Q = = < / l a t e x i t > |0i < l a t e x i t s h a 1 _ b a s e 6 4 = " c / X F 3 O W N t 3 I 3 g t X V g 1 l 4 I E E j 6 o o = " > A A A B 8 H i c d Z D L S g M x F I b P e K 3 1 V n X p J l g E V 2 W m t Z d l 0 Y 3 L C v Y i 7 V A y a d q G J p k h y Q h l 7 F O 4 c a G I W x / H n W 9 j p q 2 g o j 8 E P v 5 z D j n n D y L O t H H d D 2 d l d W 1 9 Y z O z l d 3 e 2 d 3 b z x 0 c t n Q Y K 0 K b J O S h 6 g R Y U 8 4 k b R p m O O 1 E i m I R c N o O J p d p v X 1 H l W a h v D H T i P o C j y Q b M o K N t W 7 v 3 Z 7 C c s R p P 5 d 3 C 6 5 X r J V L K A W 3 X C m n 4 B Y r V Q 9 5 F l L l Y a l G P / f e G 4 Q k F l Q a w r H W X c + N j J 9 g Z R j h d J b t x Z p G m E z w i H Y t S i y o 9 p P 5 w j N 0 a p 0 B G o b K P m n Q 3 P 0 + k W C h 9 V Q E t l N g M 9 a / a 6 n 5 V 6 0 b m 2 H N T 5 i M Y k M l W X w 0 j D k y I U q v R w O m K D F 8 a g E T x e y u i I y x w s T Y j L I 2 h K 9 L 0 f / Q K h a 8 S q F 0 f Z 6 v X y z j y M A x n M A Z e F C F O l x B A 5 p A Q M A D P M G z o 5 x H 5 8 V 5 X b S u O M u Z I / g h 5 + 0 T A v u Q k Q = = < / l a t e x i t > |0i < l a t e x i t s h a 1 _ b a s e 6 4 = " m y 1 / r m w N H e 0 E o S T e R q S b 1 V m + F B 4 = " > A A A B 7 X i c d Z D L S s N A F I Z P 6 q 3 W W 9 W l m 8 E i u C p J t Z d l 0 Y 3 L C v Y C b S i T y a Q d O 8 m E m U m h h L 6 D G x e K u P V 9 3 P k 2 T t o K K v r D w M d / z m H O + b 2 Y M 6 V t + 8 P K r a 1 v b G 7 l t w s 7 u 3 v 7 B 8 X D o 4 4 S i S S 0 T Q Q X s u d h R T m L a F s z z W k v l h S H H q d d b 3 K d 1 b t T K h U T 0 Z 2 e x d Q N 8 S h i A S N Y G 6 s z m P p C q 2 G x Z J d t p 9 K o X q A M 7 G q t m o F d q d U d 5 B j I V I K V W s P i + 8 A X J A l p p A n H S v U d O 9 Z u i q V m h N N 5 Y Z A o G m M y w S P a N x j h k C o 3 X W w 7 R 2 f G 8 V E g p H m R R g v 3 + 0 S K Q 6 V m o W c 6 Q 6 z H 6 n c t M / + q 9 R M d N N y U R X G i a U S W H w U J R 1 q g 7 H T k M 0 m J 5 j M D m E h m d k V k j C U m 2 g R U M C F 8 X Y r + h 0 6 l 7 N T K F 7 e X p e b V K o 4 8 n M A p n I M D d W j C D b S g D Q T u 4 Q G e 4 N k S 1 q P 1 Y r 0 u W 3 P W a u Y Y f s h 6 + w Q Z Q Y 9 + < / l a t e x i t > . . .

Local equilibra-on
Global purity < l a t e x i t s h a 1 _ b a s e 6 4 = " c / X F 3 O W N t 3 I 3 g t X V g 1 l 4 I E E j 6 o o = " > A A A B 8 H i c d Z D L S g M x F I b P e K 3 1 V n X p J l g E V 2 W m t Z d l 0 Y 3 L C v Y i 7 V A y a d q G J p k h y Q h l 7 F O 4 c a G I W x / H n W 9 j p q 2 g o j 8 E P v 5 z D j n n D y c s R p P 5 d 3 C 6 5 X r J V L K A W 3 X C m n 4 B Y r V Q 9 5 F l L l Y a l G P / f e G 4 Q k F l Q a w r H W X c + N j J 9 g Z R j h d J b t x Z p G m E z w i H Y t S i y o 9 p P 5 w j N 0 a p 0 B G o b K P m n Q 3 P 0 + k W C h 9 V Q E t l N g M 9 a / a 6 n 5 V 6 0 b m 2 H N T 5 i M Y k M l W X w 0 j D k y I U q v R w O m K D F 8 a g E T x e y u i I y x w s T Y j L I 2 h K 9 L 0 f / Q K h a 8 S q F 0 f Z 6 v X y z j y M A x n M A Z e F C F O l x B A 5 p A Q M A D P M G z o 5 x H 5 8 V 5 X b S u O M u Z I / g h 5 + 0 T A v u Q k Q = = < / l a t e x i t > |0i < l a t e x i t s h a 1 _ b a s e 6 4 = " u W Z W E V s C / 0 / n W U W e g T 8 X 5 Z S 9 f R s = " < l a t e x i t s h a 1 _ b a s e 6 4 = " m y 1 / r m w N H e 0 E o S T e R q S b 1 V m + F B 4 = " > A A A B 7 X i c d Z D L S s N A F I Z P 6 q 3 W W 9 W l m 8 o W c 6 Q 6 z H 6 n c t M / + q 9 R M d N N y U R X G i a U S W H w U J R 1 q g 7 H T k M 0 m J 5 j M D m E h m d k V k j C U m 2 g R U M C F 8 X Y r + h 0 6 l 7 N T K F 7 e X p e b V K o 4 8 n M A p n I M D d W j C D b S g D Q T u 4 Q G e 4 N k S 1 q P 1 Y r 0 u W 3 P W a u Y Y f s h 6 + w Q Z Q Y 9 + < / l a t e x i t > . . . < l a t e x i t s h a 1 _ b a s e 6 4 = " B f / L c T k 7 + A u D K y v s R D k G I / 8 q 6 L 0 = " > A A A B 8 H i c d Z D L S g M x F I Y z X m u 9 V V 2 6 C R b B V Z l p 7 W V Z d O O y g r 1 I O 5 R M e t q G J p k h y Q h l 7 F O 4 c a G I W x / H n W 9 j p q 2 g o j 8 E P v 5 z D j n n D y L O t H H d D 2 d l d W 1 9 Y z O z l d 3 e 2 d 3 b z x 0 c t n Q Y K w p N G v J Q d Q K i g T M J T c M M h 0 6 k g I i A Q z u Y X K b 1 9 h 0 o z U J 5 Y 6 Y R + I K M J B s y S o y 1 b u + 9 n i J y x K G f y 7 s F 1 y v W y i W c g l u u l F N w i 5 W q h z 0 L q f J o q U Y / 9 9 4 b h D Q W I A 3 l R O u u 5 0 b G T 4 g y j H K Y Z X u x h o j Q C R l B 1 6 I k A r S f z B e e 4 V P r D P A w V P Z J g + f u 9 4 m E C K 2 n I r C d g p i x / l 1 L z b 9 q 3 d g M a 3 7 C Z B Q b k H T x 0 T D m 2 I Q 4 v R 4 P m A J q + N Q C o Y r Z X T E d E 0 W o s R l l b Q h f l + L / o V U s e J V C 6 f o 8 X 7 9 Y x p F B x + g E n S E P V V E d X a E G a i K K B H p A T + j Z U c 6 j 8 + K 8 L l p X n O X M E f o h 5 + 0 T B I a Q k g = = < / l a t e x i t > |1i < l a t e x i t s h a 1 _ b a s e 6 4 = " m y 1 / r m w N H e 0 E o S T e R q S b 1 V m + F B 4 = " > A A A B 7 X i c d Z D L S s N A F I Z P 6 q 3 W W 9 W l m 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " M o j C v d 9 t I C a d z k a I o L 9 R y A r s Photonic simulation of quantum equilibration. A closed, many-body quantum system, initialized in a product state and undergoing unitary evolution generated by a Hamiltonian, necessarily remains in a pure state. However, local observables may exhibit a generalized thermalization. Entanglement builds up between sub-systems until after some time teq, each sub-system appears to have approximately relaxed into a maximum entropy state. The paradigmatic case of a non-Gaussian bosonic state evolving under a quadratic Hamiltonian can be probed via a photonic simulation platform. A fully programmable linear optical chip can provide 'snapshots' of the local and global system dynamics for arbitrary times and interaction ranges by implementing the appropriate unitary U (V ) = e −iĤt with V ∈ U(m) for m modes. tum systems, both measured in the number of optical modes and in the number of photons. The fact that photonic quantum interference without explicit photon-photon interactions carries computational hardness, as demonstrated by the hardness of boson sampling [33][34][35] shows that even non-universal photonic processors can perform operations beyond the capabilities of classical devices [29,[36][37][38]. The technological contribution of this work is to go a substantial step further and investigate to what extent the newly found levels of control and system size can be exploited for photonic quantum simulation of systems of interest, contributing to placing integrated optical devices in the realm of quantum technological devices [26,[39][40][41][42] for quantum simulation.

Local equilibration
In any setting governed by closed-system Hamiltonian dynamics, equilibration can only happen locally for local observables, since the global entropy must be preserved in time. In the setting considered, the global system is a multimode linear-optical system initially prepared in a highly non-Gaussian state ρ on m bosonic degrees of freedom, namely |ψ⟩⟨ψ| with |ψ⟩ = |1, . . . , 1, 0, . . . , 0⟩ of n = 3 single photons in m = 4 optical modes. The bosonic modes are associated with annihilation operatorsb 1 , . . . ,b m . The subse-quent integrated linear optical circuit is given by a unitary V ∈ U(m) that linearly transforms the bosonic modes. Any unitary from the group U(m) of m × m unitary matrices can be realized by a suitably designed linearly optical circuit. In state space, such linear optical circuits are reflected by ρ → σ := U (V )ρU (V ) † , where U (V ) is the physical implementation of the passive mode transformation V that linearly transforms a set of bosonic operators to a new set as The representation of the mode transformation in Hilbert space V → U (V ) is commonly referred to the metaplectic representation in technical terms. Finally, the output distribution is measured in the Fock basis using quasi-photon-number-resolving detectors, giving measurements of the form µ → P (µ) with where µ = (n 1 , . . . , n m ) is a given pattern of detection events. For our purpose of showing local equilibration, we interpret the evolution U (V ) = e −iĤt as the evolution under a Hamil-tonianĤ for time t > 0, which distributes information. In the linear optical system at hand, we will implement two Hamiltonians, a quadratic bosonic translationally invariant 'hopping' Hamiltonian, resembling the non-interacting limit of a Bose-Hubbard Hamiltonian, and a Haar random transformation V ∈ U(m) corresponding to a Hamiltonian with random long-range interactions. In a fixed-size optical system, we can simulate the evolution at various times by tuning the strength of the evolution, interpreting t as scaling the strength rather than the duration of the interaction.
As the time t gets larger, increasingly longer-ranged entanglement builds up. This means that the expected moments of the local photon numbern j :=b † jb j of each of the output modes labeled j = 1, . . . , m of the state σ will increasingly, in the depth of the circuit, equilibrate and lead to a distribution that resembles that of a (generalized) Gibbs ensemble. In other words, as seen in Fig. 1, one encounters local equilibration where the reduced quantum states of a subset of the modes, or individual modes, equilibrate and take thermal-like values. Equivalently, we can say that the state will locally thermalize, in the sense that it results in the same expectation values for local observables as if the entire system had relaxed to a thermal equilibrium state.
Strictly speaking, here we observe a generalized thermalization in the following sense. The Gibbs or canonical state reflecting thermal equilibrium is given by ξ := e −βĤ /tr(e −βĤ ) for a suitable inverse temperature β > 0 that is set by the energy density. For non-interacting bosonic systems, local equilibration for subsystems consisting of several modes is instead expected to converge to a generalized Gibbs state. To be specific, here, the initial state is a product state (and hence has obviously short-ranged correlations) -albeit not being translationally invariant -and the bosonic quadratic Hamiltonian will on the one hand be translationally invariant before it undergoes a time evolution generated by U (V ) = e −iĤt (or the Haar-random V ∈ U(m)). The situation is particularly transparent whereĤ is a hopping Hamiltonian which is transla-tionally invariant. Defining the momentum space occupation numbers asN e 2πik(y−x)/mb † xby (2) one finds that the generalized Gibbs ensemble is then given by the maximum-entropy state ω given by ω := argmax{S(η) : tr(ηN k ) = ⟨ψ|N k |ψ⟩ for all k}, (3) associated with an inverse temperature per momentum mode, where S(η) = −tr(η log η) is the von Neumann entropy. For an infinite system, convergence to such a state is guaranteed [7][8][9]11], in the sense that the global pure state will remain pure, but again, all reduced states (and for that matter, all expectation values of local observables) will for most times take the values of this generalized Gibbs ensemble. For finite systems, it has been rigorously settled in what sense the state is locally approximated by such a generalized Gibbs ensemble [7,[12][13][14][15] before recurrences set in. We discuss the specifics of this mechanism in more detail in Supplementary Note 2. For the Haar-random unitaries, we still find Gaussification in expectation, creating an interesting state of affairs, as here the theoretical underpinning is less clear.
For subsystems consisting of a single bosonic mode only, canonical or Gibbs states as well as generalized Gibbs ensembles both give rise to identical photon number distributions reflecting Gaussian states: The state 'Gaussifies' in time. The situation at hand is particularly simple in the situation where the expectation value of the photon number is the same for each of the m output modes. Then for a Gaussian state, the probability of observing k photons reduces to where D := n/m is the photon density per mode, which acts as an effective temperature. Interestingly, generalized Gibbs ensembles are still not quite thermal or canonical Gibbs states, which would be maximum-entropy states given the expectation value of the energy, but a generalization of that state, due to the noninteracting nature of the Hamiltonian. For example, in full non-equilibrium dynamics under large-scale interacting Bose-Hubbard Hamiltonians (as can be probed with cold atoms in optical lattices [21]) one expects an apparent relaxation to a Gibbs state. In contrast, a generalized Gibbs ensemble maximizes the von-Neumann entropy under the constraint of the energy expectation and the momentum space occupation numbers which are preserved under the non-interacting translationally invariant evolution t → e −iĤt . Therefore, one can say that each of the momentum modes is then associated with its own temperature, as sketched in Fig. 1, and the system 'thermalizes' up to the constraints of the momentum space occupation numbers being preserved.
Such generalized Gibbs ensembles are also interesting from the perspective of quantum thermodynamics [19,43,44]. The presence of the additional conserved charges indeed alters the thermodynamic properties and comes in as a further constraint. It is also found that the minimum-work principle can break down in the presence of a large number of conserved quantities [43]. Resource theories for thermodynamic exchanges of non-commuting and hence non-Abelian observables are also strongly altered for generalized Gibbs ensembles compared to their thermal counterparts [44].

Certification
In this section, we lay out the certification tools that we have developed to verify that the experiment has worked close to its anticipated functioning. Crucially, time evolution preserves the purity of a quantum system; the system only appears to be equilibrated when considering the local dynamics. Therefore, in the ideal case, it should be possible to undo the time evolution after applying U . This leads to the evolution U † U = I, meaning that a revival of the initial, non-Gaussian state is observed. In a noiseless experiment, this operation would function perfectly, and all entanglement will be formed between the photons as opposed to between the photons and the environment. This latter form of entanglement corresponds to decoherence and cannot be time-reversed by acting only on the photons. Therefore, the extent to which one observes a revival of the initial state serves as a measure of the degree of photon-photon entanglement versus the degree of decoherence.
We further formalize this idea in the form of a fidelity witness [45,46] that certifies the fidelity F (σ, |ψ t ⟩) = | ⟨ψ t | σ |ψ t ⟩ | between the experimentally prepared state σ and a pure target state described by a state vector |ψ t ⟩ := e −iĤt |ψ⟩. The procedure requires a well-calibrated, programmable measurement unitary and number resolving (but not spectral-mode resolving) detectors. It consists of two settings for the measurement unitary: the inverse of the target unitary and the inverse followed by a Fourier transform U F . The constant number of measurement settings and polynomial classical computation resources required mean the procedure is efficiently scalable to arbitrary system sizes. Here, we consider the specific case of witnessing against the specific target state of our experiment leaving the generalization to the supplementary material.
For the first measurement, we measure the state U † σU in the photon number basis. More specifically, we measure the fraction p 1 of detection events which correspond to our input state (i.e., exactly one photon in the first three input modes and no photon in the fourth mode). If our photodetectors would perfectly resolve the temporal and spectral degrees of freedom of the photons, this measurement in itself would be sufficient for certification [45]. However, in our system, the detectors only resolve the spatial mode. Neglecting this and naively carrying out the above procedure could result in certifying a large fidelity even with photons in distinct temporal modes, i.e., distinguishable states.
To rule this out, we employ an additional measurement setting, as part of a two-step certification process: we implement U † followed by a Fourier transformation and count photons. From the first setting, we upper bound the probability p 1 of seeing one photon in each of the first three spatial modes and no photon in the fourth. From the second, we upper bound p 2 , the overlap probability of σ with the distinguishable subspace. This is done by monitoring the fraction of observed interference patterns that would be forbidden for truly indistinguishable photons following a Fourier transform [47][48][49].
In this way, we arrive at a fidelity bound of the form where ϵ > 0 is the probability that the bound is correct, and δ is the corresponding statistical penalty, which arise from the observed photon counting statistics on p 1 and p 2 . This bound is derived from Chebyshev's inequality and holds with very few assumptions on the underlying distribution (for a full derivation, see the supplemental material). If one is merely interested in establishing the presence of entanglement in the system, one can derive a simple entanglement witness W from the estimated fidelity. We use the following definition of an entanglement witness [50] where λ 2 max is the maximal Schmidt coefficient in the decomposition of |ψ t ⟩ over a given partition, whose classical computation is not scalable, but feasible in our case. It follows then that F > λ 2 max is a witness of entanglement.

Integrated photonic platform
We use an integrated quantum photonics architecture as our experimental platform (see Fig. 2). Integrated quantum photonics constitutes a platform for non-universal quantum simulation based on bosonic interaction between indistinguishable photons [27,[51][52][53][54][55]. In integrated quantum photonics, quantum states of light are fed into a large-scale tuneable interferometer and measured by single-photon-sensitive detectors.
Our interferometer is realized in silicon nitride waveguides [56,57], and has an overall size of n = 12 modes and an optical transmission of 2.2 − 2.7 dB, i.e., 54% − 60% depending on the input channel. Reconfigurability of the interferometer is achieved by a suitable arrangement of unit cells consisting of pairwise mode interactions realized as tuneable Mach-Zehnder interferometers [58]. Each unit cell of the interferometer is tuneable by the thermo-optic effect. For a full 12 mode transformation, the average amplitude fidelity where U set and U get are the intended and achieved unitary transformations in the processor, respectively. The processor preserves the secondorder coherence of the photons [57].
We implement a quantum simulation of thermalization and a verification experiment in two separate sections of the interferometer. These two sections are indicated in blue and yellow, respectively in Fig. 2; the area below the dotted line in Fig. 2 is not used. These two sections both form individual universal interferometers on the restricted space of four optical modes, allowing us to apply two arbitrary optical transformations U 1 and U 2 in sequence.
We use the first section to simulate time evolution of our input state. We select two families of Hamiltonians to simulate: A hopping HamiltonianĤ NN = γ kb † kb k+1 + h.c. which consists of equal-strength nearest-neighbour interactions between all modes, which simulates the superfluid, non-interacting limit of the Bose-Hubbard model, and a set of 20 randomly chosen long-range HamiltoniansĤ LR = i,j γ i,jb † jb i + h. c., which we generate by applying the matrix logarithm to a set of Haar-random unitary matrices [59].
The second section of the interferometer, indicated in Fig. 2 in yellow, is used for certification. When we wish to directly measure the quantum state generated by the first section, we set this area to the identity, leaving the state after U 1 untouched. However, we can also use this second section to make measurements in an arbitrary basis on the quantum state generated by U 1 , which allows us to certify the closeness of our produced quantum state to the ideal case.
Our photon source is a pair of periodically poled potassium titanyl phosphate (ppKTP) crystals operated in a Type-II degenerate configuration, converting light from 775 nm to 1550 nm [60], with an output bandwidth of ∆λ ≈ 20 nm. By using a single external herald detector and conditioning on the detection of three photons after the chip, we post-select on the state vector |ψ⟩ = |1, 1, 1, 0⟩ [52]. By tuning the relative arrival times of our photons, we can continuously tune the degree of distinguishability between our photons. Onchip measurements via the Hong-Ou-Mandel (HOM) effect [61] lower bound the wave function overlap between photons x = | ⟨ψ i |ψ j ⟩ |, according to V = x 2 , where |ψ i ⟩ is the wave function of photon i, and V is the visibility of the HOM dip. We measure visibility's of 89% and 92% for photons of different sources, and 94% for photons of the same source. Photon detection is achieved with a bank of 13 superconducting single-photon detectors [62,63], which are read out with standard correlation electronics. For each of our four modes of interest, we multiplex three detectors to achieve quasi-photon number resolution [64] with the thirteenth detector used as the herald. By means of adjusting the time delay between the photons we can adjust their degree of mutual distinguishability. We can switch between indistinguishable particles, which produce an overall entangled state (i.e., exhibiting both modal and particle entanglement), which will exhibit thermalization, and distinguishable particles, in which each photon traverses the experiment unaffected by the others, corresponding to a product state of the single-photon wave functions, which does not exhibit local thermalization. are spontaneously split in two red photon pairs. One of these four photons is used as a herald and the other three are injected in the first three modes of our 12 × 12 integrated photonic programmable processor. The processor output is sent to small fibre-beam-splitter networks and superconducting nanowire single photon detectors (SNSPDs), which act as pseudo-number counting detectors. In the processor, we program the unitary U1 used to simulate the temporal dynamics (the blue block). In addition, we can program a second unitary U2 for the verification process (the yellow block). The zoom-in shows a Mach-Zehnder interferometer that implements one of the programmable beam-splitters. The inset shows a photograph of a fibre-connected integrated optical chip nominally identical to the one used in the experiment. Photo credit for the inset photo: Gijs van Ouwerkerk (PHIX Photonics Assembly).
tion time indicated at the head of the column, and the rows indicate different measurement settings, i.e., either the experiment itself or the corresponding certification measurements. The data in these figures was acquired over 20 minutes for the photon number distribution, 320 minutes per certification measurement for the hopping Hamiltonian and 220 minutes for each certification measurement of the long-range Hamiltonian, with four-photon events (three photos in the processor plus herald) occurring at a rate of 4 Hz.
The first row of the two sub-figures displays the singlemode photon-number statistics k → p(k) as generated after the application of U in the first section of the processor. The output statistics were measured for the first output mode. The experiment was carried out for both distinguishable (blue points) and indistinguishable (red points) photons. The grey bars show the expected distribution at full equilibration given by Eq. (4). For both Hamiltonians, initially, the input state is still clearly present, as indicated by the high probability to observe exactly one photon in the observed output mode. However, entanglement builds up as time evolves, since the photons increasingly equilibrate. Consequently, for the indistinguishable photons, the initial input state evolves to a thermal-like state at t = 1. For both Hamiltonians, the distinguishable photons (whose output statistics correspond to those of classical particles) do not approach the canonical thermal state, demonstrating the intrinsic link between entanglement and thermalization.
For the hopping Hamiltonian, at later times (t = 2, t = 5), the finite size of our Hamiltonian gives rise to a recurrence, i.e., the state moves away from equilibrium again and evolves back towards the initial input state [7,8]. For the long-range Hamiltonian, in contrast, the long-range interactions mean that recurrences are pushed away to later time not included in the simulation. These results suggest the presence of longrange order (as opposed to structured, nearest-neighbour interactions) tends to increase the time for which a system will continue to exhibit local relaxation. Whilst this picture is intuitive, a rigorous understanding of these effects is an exciting open problem for theory and future experiments. The general agreement across a large range of randomly chosen Hamiltonians also represents strong experimental evidence for the ubiquity of these effects [5,6].
The second row of the two sub-figures shows the full output-state distribution µ → p(µ) after only the application of U , measured with indistinguishable photons. The bars in the background correspond to the expected distributions. For the long-range Hamiltonian, a single representative example of our 20 Hamiltonians is plotted. From this data, it can be clearly seen that at the point of thermalization, the photons are spread over many possible output configurations, whereas a recurrence manifests as a transition back to fewer possible output configurations.
The third and fourth rows show the output-state distributions after the first and second certification measurement, respectively. In these rows, the output configurations which contribute positively to the fidelity witness are indicated in green, and those which contribute negatively are indicated in red. The first certification measurement undoes the entanglement generated by U and ideally only results in state vectors of the form |ψ out ⟩ = |1, 1, 1, 0⟩. The second certification Theoretical predictions (Th) are represented by bars and the experimental results (Exp) are represented by circles. The green-coloured data corresponds with outcomes that benefit the certification protocol, whereas the red data is forbidden, i.e., ideally should not occur. b) Long-range Haar-random model: In panel I, the time evolution of photon-number probability distribution in spatial output mode 1 for 20 different random Hamiltonians are plotted. The black points (squares) show the theoretical prediction for indistinguishable (distinguishable) particles, while coloured points correspond to experimental data. Panels II-IV show the observed output distributions for the first long-range Hamiltonian. These rows correspond to the output distributions of the first long-range Hamiltonian (panel II), the first certification measurement U −1 (panel III) and the second certification measurement UF U −1 (panel IV). Theoretical predictions (Th) are represented by bars and the experimental results (Exp) are represented by circles. The green-coloured data corresponds with outcomes that benefit the certification protocol, whereas the red data is forbidden, i.e., ideally should not occur. measurements also applies a three-mode Fourier to the generated states. Ideally, this results in only four allowed output configurations. These certification measurements show good agreement with the ideal allowed states, demonstrating the high degree of control over the experiment. For the second certification measurement, most of the deviations from the expected distribution can be attributed to the known photon indistinguishabilities. From the data presented in the third and fourth row, we extract the values of p 1 and p 2 , respectively, which are used in the fidelity witness as laid out in Eq. (5). Fig. 4a) and 4b) show the certified fidelities for both the hopping Hamiltonian and the first random long-range Hamiltonian, respectively. The three horizontal ticks on each data point correspond to confidence values of ϵ = 0.7, ϵ = 0.8 and ϵ = 0.9. The line shows the entanglement witness, corresponding to a bi-partition between mode 1 and the remaining modes. The relatively constant fidelity to the target global state contrasts against the conversion of the local, single-mode statistics to thermal statistics, as seen in Fig. 3. Fig. 4a) shows that entanglement is certified for t = 1 in the hopping Hamiltonian system. The observed fidelity F = 0.359 is above the threshold of the entanglement witness. Similarly, Fig. 4b) shows a unambiguous certification for the first long-range Hamiltonian at t = 2. The fidelity F = 0.360 is well above the certification threshold. Both of these entanglement certifications hold with a confidence of at least 90%.
The certification fidelities are limited by imperfect control over the processor. This follows from the certification fidelity at t = 0.2 for the long-range Hamiltonian. This fidelity F = 0.462 is significantly higher than others. Closer inspection shows a near optimal value for p 2 , which is now only limited by the partial distinguishability of the generated photons. This implies that the certification at other time steps is limited by imperfect chip control, i.e., a limited fidelity at which any measurement can be implemented. A second factor limiting the certification is detector blinding, which affects the obtained values of p 1 (see the supplemental material for more details on detector blinding and the convergence of the certification statistics).

DISCUSSION
In conclusion, we have experimentally shown that a pure quantum state in a closed environment can locally behave like a thermal state because of entanglement with the other modes. To this end, we simulated both the non-interacting limit of a Bose-Hubbard hopping Hamiltonian and 20 random long-range Hamiltonians on a programmable 12-mode photonic processor. Previous experiments in this direction have not been able to show this kind of reversibility since creating a sufficiently isolated quantum system and controllable evolution is notoriously difficult. However, our experiment is fully time-reversible, just like quantum mechanics itself. This reversibility has allowed us to certify that equilibration and thermalization are due to entanglement between the quantum particles rather than with the environment. These results also provide experimental evidence for the universality of these phenomena and shed new light on the role of long-range interactions on relaxation dynamics. From the point of view of the development of quantum technologies, these experiments showcase the degree of control, low decoherence and rapidly growing size of integrated quantum photonic processors as instances of a near-term quantum computational platform.

Photon source and input state preparation
Distinguishable and indistinguishable photonic quantum state vectors of the form |ψ⟩ = |1, 1, 1, 0⟩ are generated by a multi-photon source consisting of two free-space Type-II SPDC sources. Two non-linear 2 mm length ppKTP crystals (Raicol Crystals) are pumped by a Ti:Sa mode-locked laser (Tsunami, Spectra Physics) at 775 nm with a spectral bandwidth of 5.4 nm FWHM. Pulses are generated with a repetition frequency of 80 MHz and 150 fs pulse duration. Each crystal is pumped by approximately 10 mW pump power, generating degenerate signal-idler pairs at 1550 nm with generation probability < 1% per pulse. Typical heralding efficiencies for individual crystals are around 40-45%, while typical two-photon event rates are ∼ 0.20 MHz coincidence counts at 40 mW pump power. While the source is designed to produce as pure photons as possible, residual energy and momentum conservation result in spectral signal-idler correlations. These correlations are attributable to the periodically-poled structure of the non-linear crystals. We suppress these correlations by using a spectral bandpass filter of ∆λ = 25 nm. Halfwave plates are used to remove the distinguishability in photon polarization and to match the TE mode supported by our quantum photonic processor. Three motorized linear stages (SLC-2475, Smaract GmbH) are used to control relative photon arrival times, used to switch the distinguishability of the photons.

Quantum photonic processor
Our quantum photonic processor consists of a photonic chip, the control electronics which actuate this chip, and peripheral systems such as cooling. The photonic chip implements arbitrary linear optical transformations on 12 waveguides. The waveguides are implemented as stoichiometric silicon nitride (Si 3 N 4 ) asymmetric double-stripe (ADS) waveguides with the TriPleX technology [56]. The waveguides are optimized for light of a wavelength of 1550 nm, and have propagation loss of < 0.1 dB/cm. The waveguides have a minimal bending radius of 100 um. Coupling on and off the chip is achieved by adiabatic mode converters, which are implemented by removing the top layers of the ADS stack. These converters have coupling losses down to 0.9 dB / facet. The overall measured loss budget of the processor is 2.5 ± 0.2 dB, with roughly 1.8 dB attributable to the two adiabatic couplers and 0.7 dB to propagation losses on chip.
Universality of the optical transformation is achieved by a network of beam splitters in a checkerboard geometry. Each tuneable beam splitter is implemented as a Mach-Zehnder interferometer (MZI), with two static 50/50 directional couplers. To tune the MZI, two thermo-optical phase shifters are used, one inside the MZI which enables shifting of light amplitude between adjacent optical modes, and one external to the MZI which allows for a relative phase shift between the two modes. The thermo-optic phase shifters are implemented as 1 mm long platinum heaters, and have V π = 10 V, and dissipate roughly 400 mW of power each. This power is carried off the chip through a Peltier element which is itself actively cooled with water cooling. A bank of 132 digital-to-analog converters converts signals from a control computer to voltages over the heaters. A dedicated software package is used for communication, and to compute the required voltages. Control over the processor to the precision required in this experiment requires understanding of the crosstalk between these control channels, which is achieved in a dedicated software package.

Photon detection system
A suite of 13 superconducting nanowire single-photon detectors (SNSPDs) is used for photo-detection. These detectors are biased close to their critical current (8 to 22 µA range), operating at quantum efficiencies of around 90% for 1550 nm photons with typical 200 Hz dark counts. Fourfold coincidence rates within a 750 ps window are monitored by a time tagger device (Timetagger Ultra, Swabian). From the combination of photon generation rates and dark count rates, we estimate that less than one in a million measured four-fold coincidence events are expected to be triggered by a dark count. Polarization maintaining (PM) fibres are used in combination with polarization controllers to optimize and stabilize output counts in each channel. Pseudo-number resolution detection is realized by multiplexing detectors in a 1-to-3 quasi-photon number resolving detector (q3PNRD) configuration by fibre beam splitters on the four optical modes of interest, with the thirteenth detector used as a herald.

Photon detection calibration
In order to sample from µ → P (µ) in an unbiased way, as required in this work, it is important to characterize the relative output losses from the different detectors. The SNSPDs have variation in their detection efficiency, and the same holds for the output coupling of the various optical modes of the photonic processor. Non-uniformity in the overall detection efficiency of our experiment biases the sampling of P (µ), since it will suppress some outcomes while relatively enhancing others. Note that this does not hold for any inhomogeneities in the in-coupling, due to post-selection. Furthermore, we assume that on-chip losses are reasonably uniform, which is evidenced by the high matrix amplitude fidelities. Furthermore, note that an absolute detection calibration (a notoriously difficult problem at the single-photon level) is not necessary, only a relative one between the 12 detectors of interest.
Non-uniform detection channel losses are characterized by directly transmitting heralded single photons from input mode 1 to all four output modes consecutively; these optical transformations can be performed with high fidelity. In each of these four consecutive experiments, the heralded singles count rate of each detector in the q3PNRD behind the output mode of interest is measured. All measured heralded singles count rates S i originate from the same on-chip uniform heralded single photon rate R 1 , therefore, it is convenient to pool all other losses such as out-coupling efficiencies, detection efficiencies and splitting ratios for each detection channel i in a lumped factor p i , to get Since we are only interested in relative efficiencies, we introduce relative weight factors for each detection channel, which are then normalized with respect to the maximum measured heralded singles rate and defined by In our experiments, we achieved excellent weight factor stability. Typically, we observed less than 1% relative fluctuations over more than 15h time span. Similar to nonuniform detection efficiency, the fact that each qPNRD is effectively less efficient when detecting multiple photons as opposed to a single photon biases the output distribution and must be corrected for. Experimentally, we measure heralded threefold coincidence rates CC p,q,r , which denote the rate at which detectors p, q and r and the herald detector fire simultaneously, normalized to the overall frequency of successful experiments. The challenge is then to convert these probabilities into an unbiased estimate of P (µ).
To compensate for q-PNR effects, we enumerate all combinations of threefold detection event which would give rise to a particular output pattern µ. For probabilistic multi-photon detection, the probability of measuring j photons behind mode i when k photons are injected is denoted P i (j|k). We note that for P i (1|1) and P i (2|2) there are three possible permutations, while for P i (3|3) there is just one permutation. More explicitly, we find where w are the weight factors determined above and w pi + w qi + w ri ≤ 1 due to incorporated losses. Since all P i (j|k) are independent probability events, we find for an estimate for P (µ) P (µ) = (p,q,r)∈µ CC p,q,r P 1 (n 1 |n 1 )P 2 (n 2 |n 2 )P 3 (n 3 |n 3 )P 4 (n 4 |n 4 ) , (13) where µ = (p, q, r) denotes all combinations of detection events contributing to the same µ and n i is the number of photons detected in a mode i for a given µ. These results are used to correct raw measurement data.

Data availability
All experimental and simulated data used in this study are available in the 4TU.ResearchData database [65].

Code availability
All data post-processing and simulation code used in this study are available in the 4TU.ResearchData database [65].

SUPPLEMENTARY NOTE 1 -DERIVATION OF FIDELITY WITNESS
A fidelity witness provides guarantee that the fidelity of some target state with an experimental output is at least a certain threshold value with at least a certain probability [45]. Here, we present a derivation of such a witness, including finite-size statistics, that is efficient in terms of experimental effort and classical computation. In the first place we can use this threshold as evidence that our global system retains approximately the same fidelity with a target pure state whilst the local systems exhibit apparent entropy increase. A natural further question that arises is, is a particularly meaningful fidelity threshold? In this experiment, where the key feature of interest is the role of entanglement in producing local entropy production, we will use previously established relationships between fidelity and entanglement (see, e.g., Supplementary Ref. [50] to establish useful benchmarks). The idea is that the fidelity between a separable state and an entangled target state vector |ψ t ⟩ cannot exceed a certain threshold, which is set by the largest Schmidt coefficient. If the fidelity exceeds that threshold, i.e., if F > λ 2 max , entanglement must be present. Because the size of the largest Schmidt coefficient decreases with the amount of entanglement, for more entangled states lower fidelities are sufficient to witness the presence of entanglement.
Ideal case: Fully-mode-resolving detectors. The fidelity between a quantum state σ and a target state σ t = |ψ t ⟩ ⟨ψ t | is defined as In our case, the target state vector is an initial state vector Thus, the fidelity can then be written as where we have suppressed the argument V for brevity, and is lower bounded in terms of photon number operators by [45] F (n) = (n + 1 −n) wheren j = ∞ nj =0 n j |n j ⟩ ⟨n j | =b † jb j are the photon number operators, whose eigenvalues are the number of photons in mode j and ⟨ jn j ⟩ = n is the global photon number. When one post-selects for a constant global photon number in each run (which we do in our experimental setup, to n = 3), then the bound simplifies to Note that the only calculation necessary is to compute the Hermitian conjugate of the given matrix U .
Real case: Spatial-mode-resolving detectors. In our experimental setup we do not have fully-mode-resolving detectors, meaning that we have no physical equivalent ofn j . In particular, our detectors can only resolve spatial modes and no temporal ones. This leaves an uncertainty regarding the exact mode of the photon after the measurement. Instead of projecting onto a unique mode (a single pure quantum state M k = |k⟩ ⟨k|), our detectors project onto a set of states, which are spread out over the temporal degrees of freedom and, without further work at least, cannot be distinguished. This uncertainty severely limits the fidelity that can be established. For three temporal modes, the measurements at each of the four spatial modes correspond to the following set of operators, where M 0 detects the absence of photons and M 1 , M 2 and M 3 measure one, two and three photons, respectively, in a given spatial mode. Since we set up our experiment such that the ideal initial and the target state consists of one photon per spatial mode, our measurement operator of interest will be M 1 and our certification scheme will be based on the measurement of where A, B, C and D label the four spatial modes. We will want to verify whether the initial state is being recovered after implementing the unitary and its inverse, i.e., that the three photons are in the same temporal mode and each in a different one of the first three spatial modes. The problem with using these measurements and naively applying the bound in Supplementary Eq. (18) is that, considering Supplementary Eq. (19), it is clear these measurements can produce the ideal click pattern even if the photons were completely distinguishable and no quantum interference or multi-photon entanglement has been present.
At this point we introduce some alternative notation that will come in handy later: The numbers in the ket-vector label the occupied temporal mode (there are three temporal modes, so the numbers go from 1 to 3) and the subscript 1 indicates that each spatial mode is occupied by one photon only (which is the case for M 1 ). For example, all photons being in the first temporal mode reads |1, 1, 1⟩ 1 = |1, 0, 0; 1, 0, 0; 1, 0, 0; 0, 0, 0⟩. The first photon in the first, second photon in the second, and third photon in the third temporal modes reads |1, 2, 3⟩ 1 = |1, 0, 0; 0, 1, 0; 0, 0, 1; 0, 0, 0⟩. Now that we have pointed out the ambiguity of just counting photons in spatial modes, and equipped with useful notation, we next turn to the question of how to overcome the uncertainty regarding the temporal or spectral degree of freedom. An answer lies in the observation that certain interference patterns can be clearly associated with non-synchronous, i.e., distinguishable, photon states (similar to a HOM dip [61]) -we call those interference patterns forbidden patterns. We make use of this effect in practice by implementing a Fourier transform U F after the unitary U and its inverse and use the overlap between the distinguishable subspace of states and the image of the forbidden patterns under the Fourier transform U F to sharpen the lower bound on the fidelity. It is worth mentioning that this overlap is not 1:1 and some ambiguity will remain. It does, however, reduce the ambiguity significantly and thereby increases the estimated fidelity in a useful way.
We now give a full derivation of the fidelity bound. The bound, as in Main Eq. (5), has two components. 1) A lower bound p 1 on seeing one photon per spatial mode and 2) an upper bound p 2 on the overlap of σ with the distinguishable subspace. The two components correspond to two different measurement settings: 1. Implement the unitary and its inverse and count photons.

Implement the unitary and its inverse, implement a Fourier transform (an interference experiment) and then count photons.
First measurement setup. For simplicity and without loss of generality we fix M A 1 to the first temporal mode, i.e., M A 1 = |1, 0, 0⟩ ⟨1, 0, 0|. The first measurement setup can be expressed as the operator product |1, 0, 0⟩ ⟨1, 0, , the probability of seeing one photon in each of the first three spatial modes (regardless of the temporal modes), can be estimated experimentally with accuracy ϵ 1 . This estimation takes the form of a lower bound p 1 -the result of the first round of measurements.
which expanded yields The first term (black) reflects the fidelity given by F = tr U † σU (|1, 1, 1⟩ 1 ⟨1, 1, 1| 1 ) . The other terms correspond to the overlap of U † σU with those states where one (blue) or two (red) photons are distinguishable. We summarize those states asP 1 andP 2 , respectively, and callP =P 1 +P 2 the distinguishable subspace. With that shorthand notation Supplementary Eq. (21), simplifies to We can calculate p 1 by counting the relative number of instances of M 1 in our first measurement setup.
Second measurement setup. Next, in order to improve the fidelity bound in Supplementary Eq. (22), we need to upper bound its second term tr(U † σUP ), which enters the fidelity bound with a minus sign. This term contains the overlap of σ with the FIG. 5. Forbidden measurement probability. This figure shows how distinguishable and indistinguishable states translate to forbidden and allowed interference patterns after the Fourier transform. Forbidden patterns are strictly suppressed for indistinguishable states (in fact they are asserted via a suppression rule), hence indistinguishable states result in allowed patterns with unit probability. Distinguishable states result in forbidden patterns with probability λi, which we compute for all n-photon distinguishable states in all distinguishable sub-spaces Pi to establish a bound on tr(U † σUP ), the overlap between the state U † σU and the distinguishable subspaceP . distinguishable subspaceP . We upper bound it by implementing a Fourier interference experiment on the first three modes. This is described by the unitary U F = U Four is the appropriate Hilbert space operator corresponding to the physical implementation of the mode transformation Ideally, the first three modes should each be occupied by one perfectly indistinguishable photon (the first term in Supplementary Eq. (21)). For a state of this form, some counting patterns corresponding to projections onto certain output states are impossible [47][48][49]. By contrast, as we will now show, situations involving anything other than the ideal case will result in some forbidden measurement outcomes with a certain probability ( Supplementary Fig. 5). Working backwards from this, we can use the observed frequency of the forbidden states to get a worst case upper bound on tr(U † σUP ). This kind of discrimination is something we could not do using only the first setup. In this case, the forbidden patterns of photon numbers are all those which correspond to state vectors not included in the set {|1, 1, 1, 0⟩, |3, 0, 0, 0⟩, |0, 3, 0, 0⟩, |0, 0, 3, 0⟩} [47][48][49]. Defining the projection onto all forbidden states as M f , the quantity obtained with our second measurement setup can be written as where p 2 is the the probability of observing a forbidden state. Using the cyclicity of the trace in the middle term we see we can also interpret p 2 as the overlap of U † σU with the image of the forbidden states based on the photon counting measurements. As mentioned before, some states within that image clearly correspond to distinguishable states, that is, states inP and other do not. To be more clear, we can express the image operator in terms of the distinguishable sub-spaces P i ∈P (i ∈ {1, 2}) as where is the probability that a state in the distinguishable subspace P i ∈P results in a forbidden pattern andP ⊥ projects onto the complementary space toP . Since tr(U † σUP ⊥ ) ≥ 0, we have To calculate these probabilities what really matters is the number of mutually distinguishable photons. For example, to calculate the probability of the forbidden state vector |2, 0, 1, 0⟩ from states in P 1 with one photon distinguishable from two other identical photons), one first computes the probability of the indistinguishable photons transforming through the Fourier transform into a state that could be transformed into a forbidden pattern by the third photon. This can be computed via Supplementary Eq. (37) derived in Supplementary Ref. [48]. In this case, that would be permutations of the outputs |1, 0, 1, 0⟩ or |2, 0, 0, 0⟩. Then one computes the probability that the final, distinguishable photon, would fall in the correct mode to produce the forbidden output. For example, if the indistinguishable photons evolve the state vector |1, 0, 1, 0⟩, then the distinguishable boson evolving to either |1, 0, 0, 0⟩ or |0, 0, 1, 0⟩ would result in a forbidden output pattern. The transition probabilities for distinguishable and indistinguishable photons through a Fourier transform are highly symmetric. In the above example, it turns out to makes no difference in which mode the one distinguishable boson initially resides. Working through each of the states in P 1 and P 2 we find they all result in the same forbidden state probabilities of λ 1 = 4 9 and λ 2 = 2 3 respectively. Substituting into Supplementary Eq. (26) we have where we have divided through by 4 9 . We can use Supplementary Eq. (28) as an upper bound on tr U † σUP = tr U † σU (P 1 + P 2 ) , which is what we originally set out to do. We find that and thus arrive at the following final expression for the fidelity bound in Supplementary Eq. (22), to get Note that this step renders the lower bound loose in general, and will tend to provide a pessimistic estimate of the fidelity. Finally, we turn to the question of finite-size statistics. Many tools have been developed for this situation, and here we will make use of a result from Supplementary Ref. [45] based on Chebyshev's inequality which states that given k independent samples and an observed fraction p 1 we can say that the 'true' probability of that outcomep 1 must satisfy where 0 ≤ Σ < ∞ is the variance of the distribution and k is the number of measurements. This can be put together to obtain Main Eq. (5) which holds with probability ϵ = ϵ 1 ϵ 2 , and δ(ϵ) = δ(ϵ 1 ) + δ(ϵ 2 ) arising from applying Supplementary Eq. (31) to the experimental observations of p 1 and p 2 .
Generalization to larger systems. This certification method can be extended to arbitrarily many modes and photons. Whilst a detailed investigation of the robustness and performance of this method is beyond the scope of this work, we briefly explain how the protocol generalizes and make some comments. The scheme can be used to certify the creation the fidelity multi-partite entangled state vectors |ψ t ⟩ = U LO (V ) |ψ⟩ created by acting an m-mode linear optical unitary on an initial state vector |ψ⟩, where the first n ≤ m modes are populated with indistinguishable photons as there exist forbidden states for arbitrary n. In fact, the technique can be slightly generalized further to any initial state where the photons are arranged in a periodic pattern. The generalised expression for the first measurement setting would read, leading to a bound where the P i are all the different sub-spaces corresponding to the existence of different numbers bosons partitioned into different distinguishable 'species'. There can be up to n of species (i.e., one distinguishable, two distinguishable, . . . , n distinguishable -if the number of species is equal to the number of photons n, then all photons are mutually distinguishable). The second measurement setting is already described for arbitrary P i and hence n in Supplementary Eq. (26) and, recalling that allows us to obtain the bound in the general case. This scheme is manifestly efficient in the number of measurement settings (two) and also scales well in terms of the the total sample size for each probability estimate. However, to evaluate the bound we naturally need to know the value of the λ i and also the set of forbidden states to determine p 2 . These calculations are a one-off cost in the sense that it need only be performed once ahead of time for any value of n and can then be used to certify all states in the corresponding class. In this sense, it is not counted in the scaling cost of the protocol, nevertheless it is a non-trivial overhead and we discuss the calculation in some more detail.
To explain things further we briefly recall some notation and results from Supplementary Refs. [47][48][49]. Let r = (r 1 , r 2 , . . . , r m ) and s = (s 1 , s 2 , . . . , s m ) be the input and output mode occupation list with i s i = i r i = n. In our case, we have m = 4 modes and n = 3 photons. Our input mode occupation has been r = (1, 1, 1, 0). A useful alternative notation for the mode occupation list is the mode assignment list d(q), which is structured in terms of photons rather than modes. Its entries represent the photons and the numerical value indicates the mode that is being occupied by that photon (the list has as many entries as there are photons as opposed to as many entries as there are modes). For example, the mode occupation list r = (2, 0, 0, 1) becomes d(r) = (1, 1, 4) (the first and second photon being in mode one and the third photon in mode four). The general expression for the mode assignment list given a mode occupation list q reads For bosons, the transition probabilities through a Fourier transform are proportional to the permanent of an n × n sub-matrix M of the m × m Fourier matrix V Four n . With this notation of mode assignment lists, we can neatly express the transition probabilities as for the case of indistinguishable photons and for distinguishable ones. The n × n matrix constructed from V Four n , referred to as M , is defined as where d j (r) is the j th element of the mode assignment list d(r) and the elements of V Four n are given in Supplementary Eq. (23). The forbidden patterns can efficiently be calculated as the strictly suppressed output states of the indistinguishable case with respect to the chosen Fourier transform V Four n . More precisely, in Supplementary Ref. [47] it has been shown that, for a given (potentially p-periodic) initial state r, final states s are suppressed through quantum interference when the criterion holds, i.e., if the above criterion holds, then the transition probability P (r, s, V Four n ) in Supplementary Eq. (37) vanishes. Having found the suppressed (i.e., forbidden) patterns, one can then go ahead and compute the λ i (probability that states in the various P i would result in a forbidden state). As explained above, for a given P i one needs to consider the probabilities that the populations of the distinguishable species can combine to result in a forbidden state. Strictly speaking, to evaluate Supplementary Eq. (35) we only need the value of the smallest λ i . Based on preliminary investigations we conjecture that the case of n − 1 indistinguishable bosons and 1 distinguishable boson is the minimal case. If it were necessary to check all of the λ i , it is not trivial to determine how many calculations this would entail as it corresponds to the problem of placing n indistinguishable objects in k indistinguishable boxes (the bosons in each species are of course mutually distinguishable, here we are using indistinguishable in the sense that the situation that the arrangement with, say, 3 bosons in species 1 and 2 bosons in species 2 is, for our purposes, equivalent to 3 bosons in species 2 and 2 bosons in species 1) for which there is no compact form. A crude upper bound for a given n and an number of species would be to count the number of ways of placing n objects in k distinguishable boxes which would upper bound the number of λ i to be calculated via n j=2 n−1 j−1 . Even if our conjecture is true, calculating a single λ i would still involve evaluating the transition probability (and hence matrix permanent) for n − 1 bosons, which is classically hard in general. Nevertheless, to our knowledge the hardness for the specific case of a Fourier transform remains open, which leaves the total complexity of this calculation unclear for the present.

SUPPLEMENTARY NOTE 2 -GAUSSIFICATION
For completeness, in this subsection we recall some previous work on Gaussification and present some numerical results illustrating approximate Gaussification for the systems considered in this work. For non-interacting quadratic bosonic Hamiltonians, such as describe linear quantum optics experiments of the kind considered here, the mechanisms for equilibration have been well studied [7][8][9][10][11]. In particular, it has been rigorously shown that systems will tend to 'Gaussify', meaning that after a sufficiently long enough time has elapsed, any subsystem (or even a block of subsystems) will converge to Gaussian, maximum entropy states and remain there. For finite number of modes and bosons we will never find the subsystems in perfectly Gaussian state, nor will they remain in such a state indefinitely. Instead the system will approximately Gaussify with the closeness of the approximation depending upon the system size.
More formally, consider a state vector |ψ⟩ ∈ H U = H S ⊗ H E of the entire 'universe' of our experiment which comprises m modes/sites which we can think of as a system S and environment E with given second moments in the creation annihilation operators ⟨b † ib j ⟩ > (i.e., the photon occupancy). Define a reduced state of a subsystem ϱ S = tr E {|ψ⟩ ⟨ψ|}. In our work, we have been focusing on the state of this subsystem and considering just a single mode, but one could also think of a larger subsystem. We are then interested in the dynamics as a function of time and system size. Taking the thermodynamic limit will involve fixing ratio of photons to modes (n/m) and then considering the limit m → ∞. In that case we want to know if ϱ S will eventually equilibrate to ϱ mc S , the micro-canonical state on the subsystem given the constraints on the second moments. In this limit, for fixed second moments, the maximum entropy state is a Gaussian state ϱ me S = ϱ G . Approximate Gaussification can then be expressed as the condition that, for any S and any ε > 0 there exists an m and relaxation and recurrence times t Rec and t Relax , such that Even for the modest system sizes at play in this work, it is still possible to observe substantial Gaussification. In Supplementary  Fig. 6 we numerically calculate the Wigner functions for each of the 4 modes for the initial input states at t = 0 and an evolved state at t = 1 for both the hopping Hamiltonian and one of the Haar-random, long-range Hamiltonians. Whilst the initial state exhibits substantial non-Gaussianity and Wigner negativity for the modes initially occupied with a single photon, after evolution all modes appear as approximately Gaussian with additional modulation caused by finite-size effects and all Wigner negativity has vanished. It is interesting to note that mode 4, which has initially been in a perfectly Gaussian vacuum state, is technically less Gaussian after evolution. Nevertheless, the system still exemplifies the phenomena described in Supplementary Eq. (41), as Gaussification is a claim that must hold for all subsystems and not just one. After evolution the condition that all modes can be well-approximated by a Gaussian is satisfied, whereas it is radically violated by the initial input state.
FIG. 6. Approximate Gaussification. Theoretical Wigner functions of each mode (columns) for the initial input state at t = 0 (top row) and the evolved state (t = 1) for the hopping Hamiltonian (middle row) and a long range Hamiltonian (bottom row). The Wigner quasiprobability W is plotted as a function of dimensionless continuous momentum p and position q eigenvalues.

SUPPLEMENTARY NOTE 4 -QUANTUM PHOTONIC PROCESSOR OPERATION FIDELITY
To test the fidelity of the implemented optical transformations, we perform a calibration experiment. For this, classical CW light from a 1550 nm super luminescent diode (Thorlabs S5FC1005P) is injected into the input modes and an array of calibrated photodiodes (Thorlabs FGa01FC) are used to detect the output signal. The amplitude fidelities are defined as F := 1 n Tr(|U † set ||U get |), where U get denotes observed transfer matrix, U set is the target transfer matrix and the absolute signs indicate the element-wise absolute value of the matrix elements, and n = 12 modes is the size of the transfer matrices. For a set of 150 random permutation matrices, a value of F = 0.992 ± 0.002 is found, whereas for a set of 100 Haar-random matrices we find F = 0.979 ± 0.01. The full histograms of these measurements are shown in Supplementary Fig. 8.

SUPPLEMENTARY NOTE 5 -PHOTON DETECTOR BLINDING
Although post-selection on heralded three-photon events allows for extracting events based on the input state vectors |ψ⟩ = |1, 1, 1, 0⟩, other, unwanted states are frequently produced because of the probabilistic nature of SPDC sources, e.g., when one source produces a photon pair but the other one does not. Unwanted by-product states cause detector blinding, which in combination with imperfect matrix fidelities biases observed photon statistics. This effect is illustrated in Supplementary  Fig. 9 for the identity matrix transformation. Therefore, the observed photon statistics is dependent on used pump power levels. Especially photon statistics for unitary matrix transformations close to the identity matrix transformation are affected as can be seen in Main Fig. 3, as detector blinding is more likely due to most of the light being directed to a limited set of detectors.

SUPPLEMENTARY NOTE 6 -DETAILED EXPERIMENTAL RESULTS
In this section, we provide additional information on the results of our experiments. To better illustrate the convergence to Main Eq. (4), the total variation distance between the experimental data and the canonical probability density function (from the first panel of Main Fig. 3a and 3b) is shown in Supplementary Fig. 10. Furthermore, we consider the achievable fidelity bound as a function of measurement time. The certification fidelity only becomes meaningful when the probability of error 1 − ϵ is sufficiently small. This requires prolonged measurement times to accumulate sufficient statistics. Moreover, detector blinding (see Supplementary Fig. 9) limits the average pump power for the fidelity certification measurements to only 5 mW per crystal. Such low pump power results in a fourfold coincidences rate of around 4 Hz. To prevent a bias caused by long term drift, all certified time steps are measured 'interleaved', i.e., each certification measurement is repeated throughout multiple times for short run time.
Supplementary Figs. 11 and 12 depict how the certification fidelities converge when the total measurement time is increased. Each data point is the result of a 20 minute measurement per certification step, i.e., 40 minutes per data point. The horizontal axis is linearized as 1/ √ T , where T is the measurement time in hours. Empirically, we find that the certified fidelity decreases linearly on this scale (i.e., increases linearly with −1/ √ T ). This is consistent with the independent nature of the separate experimental runs. The red horizontal dashed lines are the bi-partition (spatial mode 1 vs spatial modes 2,3,4) fidelities required to certify entanglement. Finally, the blue solid line is a linear fit through the data points. The fit is extrapolated to 100 hours of measurement time. The number at the end (left) of the fit is the corresponding maximum fidelity expected based on this extrapolation. Supplementary Fig. 11 shows the convergence of the certification fidelity for the non-interacting Bose-Hubbard, or superfluid, Hamiltonian. Each panel corresponds with one of the six simulated time steps. There are a total of 16 batches for each time step, which is almost sufficient to certify τ = 1.0 against the bi-partition with ϵ = 0.9. Therefore, we included 10 more batches measured under similar conditions to increase the fidelity for τ = 1.0 from F = 0.335 to F = 0.359. Furthermore, the extrapolations indicate that longer measurement times are not going to certify the remaining simulated time steps. Similarly, the first long-range Haar-random system's certification converges as shown in Supplementary Fig. 12. There are a total of 11 batches for each time step. Here, the simulated time step of τ = 2.0 is clearly certified against the bi-partition. Unfortunately, the other time steps will not be able to reach the required certification fidelities when the measurement time is increased.