Many-body interference in bosonic dynamics

We develop a framework to systematically investigate the influence of many-particle interference on the dynamics of generic $-$ possibly interacting $-$ bosonic systems. We consider mixtures of bosons which belong to several distinguishable species, allowing us to tune the level of many-particle interference, and identify the corresponding signatures in the time-dependent expectation values of observables. Interference contributions to these expectation values can be classified based on the number of interfering particles. Interactions are shown to generate a series of additional, higher-order interference contributions. Finally, based on a decomposition of the Hilbert space of partially distinguishable bosons into irreducible representations of the unitary group, we determine some spectral characteristics of (in)distinguishability.

Whether or not quantum objects can be told apart has profound consequences on their phenomenology, which can be traced back to the mathematical structures needed to appropriately describe states of indistinguishable particles. The equilibrium states of indistinguishable bosons and fermions * a.buchleitner@physik.uni-freiburg.de are long known to obey different statistics from those of distinguishable particles, but indistinguishability also affects the dynamics of many-particle systems. The effect of indistinguishability on the dynamics stems from the coherent superposition of alternative many-particle pathways leading to the same final state. The resulting many-particle interference was first revealed by Hong, Ou and Mandel (HOM) in 1987 [1] with a pair of identical photons impinging on a fair beam splitter. The perfect bunching of the photons at the output differs drastically from the classically expected distribution, which is however recovered if there exists a property (polarization, arrival time, . . . ) by which the photons can be distinguished from each other.
Experiments involving ever larger numbers of photons in multimode interferometers [2][3][4][5][6][7][8][9] have been spurred by the fact that the output distribution of such devices cannot be efficiently sampled with classical means provided the photons are fully indistinguishable [10,25]. In realistic situations, however, perfect indistinguishability is difficult to obtain since it requires total control over all degrees of freedom that might be used to distinguish the particles. As a consequence, the theory of many-body interference for particles with distinguishing degrees of freedom has been developed [11][12][13][14][15][16][17][18], and signatures of partial distinguishability have been identified in the output correlations of multi-mode interferometers [19][20][21][22][23], as well as in the occurrence of bunched [24] and suppressed [25][26][27][28][29] events. More generally, the effects of manipulating distinguishing degrees of freedom on many-particle interference have been investigated both theoretically and experimentally [30][31][32][33][34][35][36] These studies have mostly considered photonic setups, and were therefore restricted to non-interacting, scattering scenarios. However, it is clear that interference also occurs in continuously evolving systems and in the presence of interactions between the particles. Actually, quantum many-body physics typically deals with systems of identical -and therefore interfering -particles, but interactions, rather than interference, are generally seen as the main source of complexity. Developing a formalism to systematically investigate the unfolding of many-particle interference phenomena over time in general, interacting many-body systems is therefore not only of fundamental interest, but we think that it will provide new insights into long-standing problems in many-body physics such as dynamical equilibration after a quench [37][38][39], formation and spreading of correlations [40,41], or transport in interacting systems [42,43]. So far, only small systems with two modes or two particles were studied from this perspective: e.g. HOM-like interference [44][45][46][47][48][49], the dynamics of a bosonic Josephson junction [50,51] or two-particle quantum walks [52][53][54][55][56].
Hence, it is the purpose of this work to systematically explore the impact of particle (in)distinguishability on the time evolution of generic interacting many-body systems. We consider bosons which occupy a discrete set of coupled modes and whose mutual (in)distinguishability is controlled by an additional "internal" degree of freedom, as we put forward in Ref. [60]. In Sec. II, we introduce our framework by considering the simple case of two (in)distinguishable bosons in a double-well potential. We relate their dynamics in the absence of interaction to the HOM effect, and explain the additional features that appear when the particles interact with each other. Following this introductory section, we proceed in Sec. III with a systematic classification of the indistinguishability-induced contributions to the expectation values of observables. This formalism is then applied to the study of time-dependent expectation values in noninteracting (Sec. IV A) and interacting systems (Sec. IV B). In both cases, we consider general multi-mode quantum systems and use the Bose-Hubbard model to illustrate our results. Finally, in Sec. V, we propose a complementary approach based on a decomposition of the Hilbert space of partially distinguishable particles into irreducible representations of the unitary group. This allows us to identify spectral signatures of (in)distinguishability both in the presence or absence of interactions. In Sec. VI, we close with a summary of our results and provide a short outlook.

II. HONG-OU-MANDEL INTERFERENCE IN A DOUBLE-WELL
To introduce the ideas that we develop in the rest of this paper, we take the well-known Hong-Ou-Mandel (HOM) experiment as a starting point [1]. Two photons impinge on a beam splitter, which performs the following transformation of the input modes a out 1,α = cos(θ)a in 1,α − i sin(θ)a in 2,α , a out 2,α = cos(θ)a in 2,α − i sin(θ)a in 1,α .
Here, the real parameter θ defines the reflection probability R = cos 2 (θ) of the beam splitter and a in m,α (resp. a out m,α ) is the annihilation operator associated with input (resp. output) mode m = 1, 2 for a photon whose other degrees of freedom (polarization, arrival time, . . . ) are labeled by α. In contrast to the external modes m = 1, 2 which are mixed by the beam-splitter, the degrees of freedom described by α are fixed and we will refer to them as internal degrees of freedom, or species. Particles are said to be indistinguishable if they share the same internal state (i.e. belong to the same species), and distinguishable if they are in orthogonal internal states (i.e. belong to different species). As input states for the HOM experiment, we consider two photons, one in each mode, that are either indistinguishable or distinguishable from one another: Here |∅ is the vacuum state and α and β refer to orthogonal internal states. The two-particle HOM interference is revealed by a correlated measurement at the output of the device. Specifically, we consider the observable n 1 n 2 , where n m = α (a out m,α ) † a out m,α is the total number of particles in output mode m = 1, 2. The sum over internal states α ensures that this observable does not resolve the internal state of the particles, we say it is speciesblind. In the HOM setup, the expectation value (EV) of n 1 n 2 coincides with the coincidence probability, i.e. the probability to detect one particle in each output mode. Applying the beam splitter transformation (1) together with the bosonic commutation relations a m,α , a n,β = a † m,α , a † n,β = 0 and a m,α , a † n,β = δ mn δ αβ , we find For distinguishable particles, the coincidence probability is the sum of the probabilities for both photons to be reflected and for both to be transmitted. This classical reasoning fails in the case of indistinguishable particles, where the transition amplitudes associated with these two two-particle paths interfere destructively. In particular, for a fair beam splitter (θ = π/4, R = 1/2), we observe the complete suppression of coincidence events. However, we note that the EV of the particle number in only one output is insufficient to detect interference, e.g. n 1 D = n 1 I = 1. On the basis of this observation, we define k-particle observables in Sec. III A and discuss their sensitivity to many-particle interference processes. Let us introduce an analogue of the HOM experiment based on the dynamics of particles trapped in a double-well potential (and, in the rest of this paper, in larger multi-mode systems). This allows us to very naturally introduce interactions between the particles and investigate its effect on many-particle interference. We therefore consider the time-continuous dynamics of two bosons evolving with the two-mode Bose-Hubbard Hamiltonian The external modes m = 1, 2 now refer to the spatial modes localized in each potential well and we again consider an internal degree of freedom α through which the particles can be 1. (a) Two-mode photonic scattering setup. A photon (blue dot) initially in mode m = 1 is reflected from the beam splitter (gray line), i.e. remains in this mode, with probability R = cos 2 θ and it is transmitted to mode m = 2 with probability 1 − R = sin 2 θ. (b) Dynamics of a single particle in two tunnel-coupled potential wells. The particle remains in the initial well with probability cos 2 (Jt/ ) and tunnels to the other well with probability sin 2 (Jt/ ). The evolution is the same for both systems after the identification θ = Jt/ .
distinguished. J is the tunnel coupling between the two wells, U is the on-site interaction strength and n m = α a † m,α a m,α counts the total number of particles in mode m. By choosing the tunnelling and interaction parameters to be independent of the particles' species, we rendered the Hamiltonian H DW species-blind. As a consequence, the unitary evolution operator U(t) = e −iHDWt/ does not modify the internal degrees of freedom and acts on particles independently of their species. The dynamics is thus only affected by the (in)distinguishability of the particles in the initial state, which controls the occurrence of many-particle interference. This system has for example been realized experimentally with rubidium atoms trapped in a double-well potential generated by two optical tweezers and which can be distinguished by their magnetic hyperfine state [45].
In the absence of interactions (U = 0), each particle independently performs Rabi oscillations between the two modes and the double-well setup is equivalent to the HOM one. Indeed, the evolution of the annihilation operators of the two modes in the Heisenberg picture then reads a 1,α [t] = cos(Jt/ )a 1,α − i sin(Jt/ )a 2,α , (8a) a 2,α [t] = cos(Jt/ )a 2,α − i sin(Jt/ )a 1,α , with O[t] = U † (t)OU(t) (we set the initial time to t 0 = 0). By comparison with Eqs. (1), we see that after an evolution time t, the non-interacting double-well implements a beam splitter transformation with reflectivity parameter θ = Jt/ . The relation between these two processes is illustrated in Fig. 1. In Fig. 2, the non-interacting evolution of the EV n 1 [t]n 2 [t] is plotted with full lines for both states |Ψ I (black) and |Ψ D (red). In particular, we observe the cancellation of n 1 [t]n 2 [t] at t = π/4J for the state of indistinguishable particles.
In presence of interaction, the simple linear relation (8) between the annihilation operators at time zero and at a later time t breaks down. In Sec. IV B, we show that, as a consequence, EVs of observables pick up new contributions which are sensitive to higher-order interference effects. While in the non-interacting case a two-particle observable n 1 n 2 was necessary to detect (two-particle) interference, in the pres- ence of interactions this is already possible with a singleparticle observable such as the tunnelling operator H tun = J α a † 1,α a 2,α + a † 2,α a 1,α . As an example, we show in Fig. 3 the EV of H tun for both states |Ψ I and |Ψ D and for various values of the interaction strength U . While both states yield the same EV (identically equal to zero) for U = 0, the introduction of interactions allows to tell them apart.
To advance in the study of interacting systems, we develop, in Sec. V, an approach based on the structures imprinted upon the many-particle Hilbert space by distinguishability. To give a flavor of this method, let us write the Hamiltonian matrix for two interacting particles in the double-well.
If these particles are distinguishable, of species α and β, it is instructive to use the following basis: Indeed, the Hamiltonian H DW [Eq. (7)] then takes a blockdiagonal form which reflects the symmetry of the states under the exchange of α and β: one 3 × 3 block operating on the subspace spanned by the symmetric states |Ψ 11 , |Ψ + and |Ψ 22 , and one 1 × 1 block acting on the antisymmetric state For two indistinguishable particles, there are only three basis states and a short calculation shows that the Hamiltonian matrix in that basis is identical to the top-left block of Eq. (13). The state of distinguishable particles |Ψ D is a superposition of the above symmetric and antisymmetric states (10) and (12), with |Ψ + following the same evolution as the state of indistinguishable particles |Ψ I , while |Ψ − , being alone in the antisymmetric subspace, remains invariant. The EV of a speciesblind observable O also decomposes into symmetric and antisymmetric parts One can verify that the EVs of n 1 n 2 and H tun with respect to |Ψ D and |Ψ I , shown in Figs. 2 and 3, obey this relation whatever the value of the interaction strength U .
In larger systems (more particles and/or modes), the Hamiltonian describing mixtures of particles from different species also decomposes into blocks corresponding to specific symmetries under exchange of the species, as we discuss in Sec. V A. Although it always contains a block which coincides with the Hamiltonian of indistinguishable bosons, we will see that distinguishability can give rise to distinct new frequencies in the dynamics. We explore such consequences of distinguishability for the many-body spectrum in Secs. V B and V C.

III. INTERFERENCE IN THE EXPECTATION VALUES OF OBSERVABLES
In this section, we systematize and generalize the concepts introduced in the previous section, beyond the simple example of HOM interference. In particular, we define a hierarchy of operators based on the number of particles that they act upon and analyse the many-particle interference contributions to their EVs.

A. Order of operators
We consider bosons which can populate a set of external modes m = 1, . . . M and belong to one of several species α = 1, . . . S. The operator a † m,α creates a particle of species α in mode m and obeys the bosonic commutation relations (4). To systematically study the influence of (in)distinguishability on the EVs of observables it is useful to classify operators according to the number of particles that they act upon (which we also refer to as the order of the operator). Moreover, we want to consider operators which act solely on the external degrees of freedom of the particles, i.e. that are species-blind in the sense that they act on particles regardless of their species and without modifying it.
A species-blind one-particle operator (1PO) can be built as a linear combination of operators of the form S α=1 a † m,α a n,α , which move one particle from mode n to m, irrespective of its species. The general structure of a 1PO is thus given by the expression where the o n m are complex matrix elements. Analogously, we define a species-blind two-particle operator (2PO) as a sum of terms with two creation and two annihilation operators: a † m,α a † m ,α a n,α a n ,α .
Here and in the following, we omit the bounds of the sums when there is no ambiguity. Note that we use normal ordering, i.e. annihilation operators are written to the right of creation operators. To ensure species-blindness of the 2PO in Eq. (19), we sum over species indices α and α , each of which appears as subscript of both a creation and an annihilation operator. This construction can be generalized to species-blind k-particle operators (kPO), where the sum is over mode vectors m = (m 1 , . . . m k ) and n = (n 1 , . . . n k ), with m i , n j ∈ {1, . . . M }, and species vectors α = (α 1 , . . . α k ), with α i ∈ {1, . . . S}. We represent such operators diagrammatically by a box with k legs on the left and k legs on the right, corresponding, respectively, to the creation and annihilation operators, as shown in Fig. 4(a). If the coefficients obey o n m = (o m n ) * then the operator is an observable. For simplicity, we will use the abbreviation kPO both for k-particle operators and k-particle observables, as the distinction should be obvious from the context. Examples of species-blind 1POs are the particle density in mode l, or the tunnelling Hamiltonian on a one-dimensional (1D) chain, with tunnelling strength J, On the other hand, an on-site interaction of the form is a species-blind 2PO. Finally, the probability of measuring a given N -particle configuration N = (N 1 , . . . N M ), with N m particles on site m, independently of their species -a quantity often considered in the many-photon interference context -can be expressed as the EV of the species-blind N -particle observable Higher-order operators are naturally obtained as products of lower-order operators. Because of the normal-ordering requirement, the product of a kPO and an lPO gives rise to operators of orders max(k, l) through k + l, i.e.
as can be seen by application of Wick's theorem [57]. Diagrammatically, the contraction of an annihilation operator of O k with a creation operator of P l corresponds to joining the corresponding legs. The diagrams resulting from the product of two 2POs are shown in Fig. 4(b). The terms of order k + l in the products O k P l and P l O k are identical since they are obtained by normally ordering all creation and annihilation operators appearing in O k and P l . They correspond to diagrams such as the last one shown in Fig. 4(b), where no legs are contracted. As a consequence, the commutator [O k , P l ] only contains terms of orders max(k + l) through k + l − 1.

B. Expectation values of operators and order of interference
The hierarchy of operators defined above is relevant for the study of (in)distinguishability since a k-particle observable can only be sensitive to interference processes involving at most k particles. In particular, the EVs of 1POs are independent of the particles' mutual (in)distinguishability, while second-and higher-order operators are affected [58]. In the following, we derive this result within our framework.
We consider EVs in many-body Fock states with N m,α particles of species α in mode m: We denote N = ( α N m,α ) m=1,...M the corresponding populations of the modes and S = ( m N m,α ) α=1,...S the species distribution.
The EV of a kPO [Eq. (20)] in state |Ψ only receives contributions from terms associated with mode vectors m and n, and species vectors α, all of length k, such that there exists a permutation that sends m to n without changing α. This condition ensures that, after successive action of the annihilation and creation operators on the Fock state |Ψ , one recovers the same state, multiplied by the number N Ψ m,α of ways of picking k particles, with mode and species labels This leads us to the following expression for the EV of an arbitrary species-blind kPO: where is the set of unique mode vectors n generated by a permutation of m, under the restriction that this permutation does not change the species vector α. In group-theoretic language (see e.g. [59]), S α (m) is the orbit of m generated by the stabilizer of α in the group S k of permutations of k objects. Diagrammatically, one assigns one particle of |Ψ to each left leg of the observable (thereby defining vectors m and α). The right legs are then associated to the same species vector α but to a potentially different mode vector n. If there exists a way of connecting left and right legs, pertaining to the same mode and species, one-to-one, then the corresponding matrix element o n m contributes to the EV. Indeed, by connecting the legs one defines a permutation π ∈ S k such that m i = n π(i) and α i = α π(i) . All particles belonging to a given cycle [59] of π must belong to the same species. Note that there may be more than one such permutation, and that this does not affect the EV. The order of the corresponding interference process can thus be defined as the minimum, over such permutations π, of the longest cycle length of π.
In In diagrams (a) and (c), the mode vectors m and n are identical, such that the identity permutation can always be used to connect the legs, independently of the species of the particles (represented here by the coloured lines). These ladder diagrams therefore correspond to classical contributions which involve no many-particle interference. In panels (b) and (d), a pair of modes is exchanged between m and n. If these modes are populated by particles of the same species, one can draw crossed diagrams, which correspond to a two-particle (HOMlike) interference process. In diagram (e), m and n are related by a cyclic permutation, corresponding to a three-particle interference process.
Evaluating Eq. (28) for an arbitrary species-blind 1PO [Eq. (18)] yields where N m = α N m,α is the total number of particles in mode m. Hence, as expected, the EV of a 1PO is completely independent of the species of the particles. In contrast, the EV of an arbitrary species-blind 2PO does depend on the species distribution. It reads In the first term, we have gathered the diagonal matrix elements o mm mm , associated with ladder diagrams as represented in Fig. 5 (a). This term depends only on the total populations N m and is therefore independent of the distinguishability of the state. The second term contains elements o m m mm where the order of the lower and upper indices is related by a transposition π = (12) (using cycle notation [59]). These contributions are associated with crossed diagrams as shown in Fig. 5 (b) and only occur if mutually indistinguishable particles are present in modes m and m .
For the EV of a species-blind 3PO, we can again distinguish contributions associated to mode vectors m and n that can be related by i) the identity, ii) a transposition, iii) a cyclic permutation, yielding (we use 1, 2, 3 as a shorthand for m 1 , m 2 , m 3 to lighten the notation) These contributions can thus be seen as resulting from i) no interference, ii) two-particle interference, iii) three-particle interference. Diagrams representing these processes are drawn in panels (c), (d) and (e) of Fig. 5, respectively. From this reasoning, it is clear that the EV of a kPO will in general receive contributions from interference processes involving up to k particles, in addition to the classical ladder terms, which are the only ones to survive in a state of distinguishable particles. Before closing this section, we remind the reader that we have only considered EVs in Fock states of the form (26). Although this is a rather tight restriction at first sight, note that we have total freedom in the choice of the single-particle mode basis (which does not change the generic form of species-blind kPOs). Moreover, for superpositions of Fock states with different numbers of particles per species, the EV of species-blind observables is additive; i.e. for |Ψ = i c i |Ψ i , with |Ψ i of form (26), we have Ψ|O|Ψ = i |c i | 2 Ψ i |O|Ψ i as long as each |Ψ i has a distinct species distribution S. The absence of cross terms follows from the fact that species-blind observables do not modify the internal degrees of freedom of the |Ψ i , which are and remain in orthogonal states. This allows to easily treat the case of particles with non-orthogonal internal states. For example, the input states |Ψ I and |Ψ D of the HOM experiment [Eq. (2) and (3)] can be combined into a state of partially distinguishable particles where the particle in the second mode is in a superposition of internal states α and β. By changing the value of η, one simply interpolates between the distinguishable and indistin- guishable cases: C. The degree of indistinguishability The considerations of the previous section provide a convenient way of quantifying the indistinguishability of a given Fock state |Ψ , namely by comparing the contributions to the EV of an observable which result from many-particle interference to those which do not. In Ref. [60], we thus defined a measure of the degree of indistinguishability (DOI) based on the relative numbers of crossed and ladder contributions to the EV of a 2PO [see Eq. (31)]: A state of indistinguishable particles, i.e. a Fock state with all particles belonging to the same species, has a DOI value of one. On the other extreme, a Fock state where each particle belongs to a distinct species has a DOI value of zero. Actually, it suffices that no two particles of the same species occupy distinct modes for the DOI to vanish. Generic Fock states (26) have a DOI value between zero and one, depending on the specific distribution of particles among modes and species. For example, in Fig. 6, we depict three Fock states with identical repartition of the particles over the modes but different distributions of the species. State |Φ 1 is fully indistinguishable and, accordingly, has a DOI value of one. State |Φ 3 gives rise to no many-particle interference, having a DOI value of zero. The partially distinguishable state |Φ 2 has I = 3/11. The proposed DOI measure depends only on the Fock state |Ψ , and not on the system's dynamics nor on the measured observable. It evaluates the magnitude of interference contributions to EVs of 2POs but says nothing of the constructive or destructive character of the interference, which depends on the observable through the complex coefficients  (36) is based on the EVs of 2POs, which are not affected by interference processes involving more than two particles. In contrast, the EVs of higher-order observables also receive contributions from higher-order interference processes, with a different dependence on the populations N m,α [see e.g. Eq. (32)] and whose importance might not well be captured by the DOI measure. Note, however, that for the states (26) considered here, two-particle interference is a prerequisite to observe these higher-order processes, since these require more than two indistinguishable particles populating distinct modes. Moreover, low-order observables, i.e. fewbody correlation measurements, are experimentally easier to access. We therefore believe that our measure is well suited to quantify indistinguishability in actual physical setups.
In the following section, we discuss, inter alia, the connection of the DOI measure I with the time-dependent EVs of various low-order observables in non-interacting [Sec IV A] and interacting [Sec. IV B] many-particle systems.

IV. INFLUENCE OF INDISTINGUISHABILITY ON TIME-DEPENDENT EXPECTATION VALUES
We are interested in the impact of the mutual (in)distinguishability of particles on their dynamics. To study this influence, we compare the evolution of Fock states |Ψ [Eq. (26)] that have the same total density distribution N but differ by the number of particles per species and/or how the different species are distributed over the modes. An example of such states is given in Fig. 6. Moreover, we consider evolution and measurements in the external degrees of freedom only, i.e. both the Hamiltonian and the measured observables are species-blind, such that any observed differences between the evolutions of the various states can be attributed to the sought-after interference effects.

A. Non-interacting systems
We consider a system of non-interacting particles, i.e. particles evolving with a 1PO Hamiltonian [compare Eq. (18)] The time dependence of the annihilation operators in the Heisenberg picture (taking t 0 = 0) is then given by a n,α [t] = U † 0 (t)a n,α U 0 (t) = n c nn (t)a n ,α , where U 0 (t) = e −iH0t/ is the non-interacting manyparticle evolution operator while the coefficients c nn (t) = (e −iht/ ) nn are the matrix elements of the single-particle evolution operator -obtained by exponentiating the singleparticle Hamiltonian matrix (h n m ) mn -and are independent of the species α. Note the explicit distinction between the full time dependence of an operator in the Heisenberg picture, denoted with square brackets, and the time dependence of the operators' matrix elements, denoted with round brackets.
As a consequence of the linear relation between annihilation operators at times zero and t (and the analogous relation for the creation operators), the order of an operator, as defined in (20), is preserved under non-interacting evolution: Hence, the time-dependent EV in a Fock state (26) will exhibit the same form as the static EV given in Eq. (28), except for the explicit time dependence of the matrix elements, as given in (41). All the conclusions that we have drawn in Sec. III therefore apply directly to the time-dependent EVs of noninteracting systems. For illustration, we consider the example of bosons on a M -site 1D lattice with nearest-neighbour tunnelling and hardwall boundary conditions, as described by the Hamiltonian (22). The upper plot of Fig. 7 shows the time-dependent variance of the particle density on the first site, ∆n 2 , in a system of M = 4 modes for the three initial Fock states illustrated in Fig. 6. As discussed in Sec. III C, the two-particle crossed terms which contribute to n 2 1 [t] all come with real and positive coefficients, so that the curves tend to order according to the initial state's DOI value I. In contrast, for the covariance cov(n , shown in the lower plot of Fig. 7, the sign of the crossed contributions is not fixed, and this ordering is lost. Larger contrast of the oscillations for states with a higher DOI however still prevails. The top panel of Fig. 7 suggests that, in a non-interacting system, one can estimate the DOI of an unknown state by comparing the variance of an on-site density n l [t] (or, equivalently, n 2 l [t] ) to the values obtained with completely distin-guishable and indistinguishable particles. To reduce the dependence of this estimation on the choice of site l and eliminate the time t, we perform a time average, which we denote by an overline. We then define which measures the excess fluctuations of the mode's population due to indistinguishability. Here |Ψ I and |Ψ D respectively denote an indistinguishable (all particles of the same species) and a distinguishable state (no particles of the same species in different modes) with the same total density distribution N as the state |Ψ , generalizing the states (2) and (3) introduced in the HOM setting. With Eqs. (41) and (31), we find Comparison with Eq. (36) shows that F(Ψ) I(Ψ) for a sufficiently narrow distribution of the coefficients |c l,m (t)c l,m (t)| 2 . To be more precise, the difference between F and I can be approximately bounded (in the sense that the bound holds for most, but not all, states |Ψ ) by [60] |F − I| σ cc µ cc min(I, 1 − I), where σ cc and µ cc are, respectively, the standard deviation and mean of the set of M (M − 1) coefficients |c l,m (t)c l,m (t)| 2 with m = m . In particular, one expects a stronger correlation between F and I in translation-invariant systems, where the transition probabilities between any two sites should be comparable for all pairs of sites, on average. For the state |Φ 2 shown in Fig. 7, the indistinguishable and distinguishable benchmark states are |Ψ I = |Φ 1 and |Ψ D = |Φ 3 , and we find F(Ψ 2 ) ≈ 0.275, which is indeed in very good agreement with the state's DOI value I = 3/11 ≈ 0.273. To confirm that this correlation also exists in larger systems and for arbitrary initial Fock states, we show in the upper (resp. lower) panel of Fig. 8 the relation between F and I in a one-dimensional lattice with M = 8 (resp. M = 12) modes and N = 2M particles. For systems with S = 2, 3 and 4 different species, we sample 10 5 initial Fock states uniformly by randomly drawing an integer between one and the total Hilbert space dimension and taking the Fock state with the corresponding basis index. The DOI I and the excess fluctuation F of the squared particle density n 2 l (for M even, results are independent of the chosen site l [61]) are calculated for each Fock state using Eqs. (36) and (42), and the resulting density histograms are shown in green. The yellow lines give the approximate bound (44).
We confirm that, for the great majority of states, we indeed have F ≈ I within the bounds set by Eq. (44). In this system, measuring the excess density fluctuations therefore provides a good estimate of the initial state's DOI. Comparison of the upper to the lower panel shows that the correlation between F and I becomes stricter as the system size M and the particle number N are increased (keeping the ratio N/M = 2 constant). In addition to the total density histograms in green, we also show the projections of the density histogram for each fixed number of species S in black (S = 2), red (S = 3) and blue (S = 4). The maxima of the distributions are roughly located at one half, one third and one fourth, respectively, which implies that a randomly sampled Fock state most likely has a DOI value of approximately 1/S. This is in agreement with the intuition that, on average, the presence of a larger number of distinct species leads to lower levels of indistinguishability (as quantified by I) and, accordingly, a weaker contribution of many-particle interferences (quantified by F) to the observed signal.
In summary, we have shown in this section that the influence of indistinguishability on time-dependent EVs in noninteracting systems can be directly understood with the operator hierarchy introduced in Sec. III A, together with the structure of their EVs discussed in Sec. III B. Indeed, the noninteracting Heisenberg time-evolution preserves the order of operators and simply introduces a time dependence in their matrix elements. The time-averaged EVs of squared 1POs, where the crossed contributions always induce constructive interference, therefore correlate very well with the DOI measure proposed in Sec. III C, providing a convenient way of measuring this value. In the following section, we discuss how the presence of particle-particle interactions alters this picture.

B. Interacting systems
To study the influence of interactions, we add a manyparticle term V to the non-interacting Hamiltonian H 0 [Eq. (37)]: Here V is a 2PO for two-body interactions, a 3PO for threebody interactions, etc., and assumed to be species-blind. The relative strength of the two terms is controlled by the real parameter . For instance, the species-blind Bose-Hubbard (BH) Hamiltonian with on-site two-body interaction U is of the above form: As noted in the previous section, the non-interacting evolution U 0 (t) = e −iH0t/ does not alter the order of operators, as defined in Eq. (20). We will now see that this is no longer true for the full interacting evolution U(t) = e −iHt/ . It is therefore convenient to introduce the interaction picture operators, identified by the subscript I, whose order remains fixed in time The Heisenberg time evolution of an arbitrary observable O then reads with the interaction picture evolution operator where the time-ordering operator → T rearranges operators on its right-hand-side in temporal order, with time arguments decreasing from left to right.
Expression (48) can be expanded in orders of , with the nth order term O (n) (t) consisting of the nested commutator of O I (t) with n time-averaged interaction operators V I (τ )dτ : 2PO contribution to O If O is a kPO and V a v-particle interaction, i.e. a vPO, then O (1) (t) contains terms of order up to k + v − 1, as explained in Sec. III A. Analogously, is obtained by applying a commutator with a vPO twice to a kPO and therefore contains terms of order up to k + 2(v − 1). Pursuing this reasoning, the highest order term in O (n) (t) is a k+n(v−1)PO. For instance, if V describes two-body interaction, then the first-order correction O  1 (t) contains 2PO and 3PO terms, etc. Diagrammatic representations of some of these terms are shown in Fig. 9.
The fact that an observable develops higher-order contributions through interaction has an interesting implication for 1POs: in the presence of interaction, their time-dependent EVs become sensitive to the mutual (in)distinguishability of the particles. In particular, to first order in , O 1 [t] Ψ receives two-particle ladder and crossed contributions through the term O For illustration, the upper panel of Fig. 10 shows the numerically exact time-dependent EV of the single-mode density n 1 [t] in a BH system of M = 4 modes. In order to break the bipartite symmetry of the lattice, which leads to the cancellation of contributions proportional to an odd power of U [60], we add a (1PO) tilt F M m=1 m n m to the Hamiltonian (46), with F = 3J. The EVs are taken in initial states with a fixed total density distribution N = (4, 2, 0, 0), but variable numbers and distributions of particles per species. The DOI value [Eq. (36)] of each state is color-coded as indicated by the color bar. For short evolution times, we see that the EVs are indeed ordered according to each state's DOI value. Similar results are obtained for other initial density distributions, as shown in the lower two panels of Fig. 10. In the inset, we plot the values of n 1 [t = /J] as a function of U for the same set of states and compare them to the analytical predictions to first order in U (dashed lines). This confirms that the splitting of the EVs for small U t/ is dominated by two-particle crossed contributions. For longer evolution times and/or larger interaction strengths, higher-order contributions to the EV become increasingly relevant and the correlation with the DOI value is degraded.
Analogously, the time-dependent EV of an arbitrary 2PO develops an infinite hierarchy of interaction-induced contributions. In Fig. 11, we show the particle density variance ∆n 1 [t] in the system which we already considered in the noninteracting case (Fig. 7), but now for a finite on-site interaction strength U = 0.25J. For short evolution times, the influence of interactions on the time-dependent EVs is limited, and the variance evolves similarly as in Fig. 7. For longer evolution times, higher-order contributions (both in the interaction strength and in the number of particles) build up and superpose, leading to distinctly different dynamics, with reduced oscillation contrast as compared to the non-interacting case.
Surprisingly, despite the strong influence of interaction on the dynamics, we can observe that the curves are still ordered according to the states' DOI value most of the time. This is a consequence of the positivity of the matrix elements    of the interaction strength and for different system sizes in the right panel of Fig. 12. The low value of µ(|F − I|) for weak interaction strengths U indicates good correlation between F and I. For stronger interactions, µ(|F − I|) displays an erratic behavior that remains to be understood.

C. Counting particles of a given species in a mode
Inspection of specific initial states leads to the interesting observation that states where all but one particles are localized on a single site have F = I. This means that the lone particle can be used as a probe to count the number of particles with identical species label in the other occupied mode. We show in the following that this is in fact true for any such initial configuration, any observable, any interaction strength and at all times. Let us denote bym the mode that contains a single particle and byα that particle's species. The remaining N −1 particles are in mode m. From our discussion of the EV of kPOs [see Eq. (28)], we conclude that, in this case, only two-particle interference can occur, and that the corresponding contributions to the EV of an arbitrary observable O[t] come with the multiplicity Nm ,α N m,α = N m,α . It follows that the EV must be of the form where the functions f and g depend on the specific observable, Hamiltonian and initial state. Hence, for all times t, the EVs of the states are ordered according to the number of particles N m,α of speciesα in mode m. In Fig. 13 we illustrate this relation for the EVs of n 3 [t] (upper panel) and n 2 3 [t] (lower panel) in a three-site BH system with initial populations N = (5, 1, 0).
Renormalizing the above EV with respect to the corresponding distinguishable state |Ψ D , where no particle of speciesα is initially in mode m, and the fully indistinguishable state |Ψ I , where all particles are of speciesα, one indeed obtains the DOI value, which in this case is simply the fraction of particles of speciesα in mode m: To sum up this section, we have learned that interactions induce a series of higher-order many-particle interference contributions in the EVs of kPOs. In particular, 1POs become sensitive to the DOI of the initial state [recall Fig 10]. Despite these additional contributions, we observed that timeaveraged EVs of squared 1POs remain well correlated with the two-particle indistinguishability measure I of the initial state for weak interaction strengths, as shown in Fig. 12. Furthermore, Fig. 13 illustrates that there are particular initial configurations where only two-particle interference can take place, independently of the observable and interaction strength.
However, as far as perturbative arguments were employed above, these are no longer appropriate for strong interactions and/or long times, where higher-order terms (both in interaction strength and number of particles) become important. In the following section, we change our perspective on the problem and identify signatures of (in)distinguishability in the time-dependence of observables based on the symmetry of the initial state under permutations of the particles. We will see that this new approach rejoins the preceding analysis in some respects, but that it often offers a complementary picture.

V. SPECTRAL PROPERTIES OF PARTIALLY DISTINGUISHABLE BOSONS
With the help of group-theoretical arguments, we first show in Sec. V A that the Hilbert space decomposes into subspaces that are invariant under the action of a species-blind Hamiltonian. We then study energy degeneracies between these subspaces in Sec. V B and explore the spectral signatures of (in)distinguishability in Sec. V C.

A. Hilbert space structure
As observed at the end of Sec. II, the Hamiltonian of two interacting particles in a double-well potential can be written in block-diagonal form provided the chosen basis vectors have a definite symmetry under the exchange of species. We now show how this decomposition can be generalized to systems with more particles and modes. This hinges on an important result from representation theory, Schur-Weyl duality [62], which relates the action of the symmetric group (which permutes particles) to that of the unitary group (which transforms single-particle states). For bosons with both internal and external degrees of freedom, the requirement that the complete (internal and external) wavefunction be symmetric under particle exchange has consequences for the structure of unitaries acting only on the external -or only on the internal -degrees of freedom. This result is known as unitary-unitary duality [62][63][64][65] and it can be stated as follows: For N bosons with external states (modes) m ∈ {1, 2, . . . M } and internal states (species) α ∈ {1, 2, . . . , S}, there exists a basis {|λ, µ, ν } of the Hilbert space -we elaborate on the meaning of the indices λ, µ and ν below -such that for all species-blind Similarly, for all "mode-blind" operators Σ (defined in analogy to species-blind operators, exchanging the roles of internal and external degrees of freedom, see Sec. III A), Here o (λ) and s (λ) are matrices representing O and Σ on the subspaces labelled by λ. Let us now be more precise about the meaning of these formulas, and in particular of the symbols λ, µ and ν.  The vector-valued label λ refers, simultaneously, to irreducible representations (irreps) of the unitary groups (and associated algebras) U (M ) and U (S) -which respectively act on the external and internal degrees of freedom-and of the symmetric group S N of permutations of N objects. The fact that a single label can be used to identify the irreps of different groups is known as "duality". These irreps can be labelled by integer partitions of N , i.e. sequences of integers λ = (λ 1 , λ 2 , . . . λ k ) with λ 1 ≥ λ 2 · · · ≥ λ k > 0 and   ing only those blocks where states with a given species distribution can have non-zero weight.
For a Fock states with N m,α particles of species α in mode m [see Eq. (26)], the known (classical) algorithms to calculate the expansion coefficients c λ,µ,ν are not efficient [66,67]. However, given the species populations S = ( m N m,α ) α=1,...S , we can try to determine for which λ and ν the coefficients c λ,µ,ν can take non-zero values.
First of all, since all symbols in the same column of a semistandard tableau must be different, the diagram λ cannot have more rows than there are modes or species. Let us now consider a state where all bosons belong to the same species. The only diagram λ that can be filled with N identical species labels according to the rules specified above consists of a single row of length N , i.e. λ = (N ) (see also Tab. I). Since for any set of species labels there is exactly one possible filling of this diagram, the corresponding subspace, of dimension dim((N ), M ) = N +M −1 N , appears once in the Hilbert space of any mixture of species. In this subspace, the Hamiltonian is identical to that of N indistinguishable bosons.
For a completely distinguishable state of N particles in N different species, all diagrams λ with at most M rows can appear (possibly several times). Provided M ≥ N , this includes the diagram consisting of a single column. The dynamics in the corresponding subspace(s), of dimension dim ((1, 1, . . . , 1), M ) = M N , is identical to that of N indistinguishable fermions.
Other Young diagrams are associated with mixed exchange symmetries, those of immanons [68], and appear in the description of partially distinguishable bosons or fermions. For a generic Young diagram λ, the number of fillings ν compatible with a given species distribution S (which can also be seen as a partition of N ) is given by the Kostka number K(λ, S) [69]. Closed expressions for the Kostka numbers are, however, not known in general. Nonetheless, the existence of this block-diagonal structure already allows to make qualitative statements about the many-body energy spectrum, as we explain in the following sections.

B. Degeneracies in systems of partially distinguishable particles
In this section, we analyse degeneracies in the spectra of species-blind Hamiltonians generating the dynamics of a mixture of different species of bosons. Although we focus on energy degeneracies, which have consequences for the dynamics of the system, the following arguments apply for any speciesblind observable.
As we have learned from the previous section, degeneracies occur independently of the form of the (species-blind) Hamiltonian when a given Young diagram λ can be filled in more than one way with the set of species labels determined by the species distribution S, i.e. when the Kostka number K(λ, S) > 1. In this case, Hamiltonian blocks belonging to different ν but the same λ have identical spectra. This can happen as soon as there are particles of at least three different species in the system. For example, in Fig. 14, the block associated with λ = (2, 1) (in red) appears twice if S = (1, 1, 1) and the corresponding energy spectra are identical.
These degeneracies are directly derived from the symmetry properties of species-blind observables and occur independently of the choice of a particular species-blind Hamiltonian. However, additional degeneracies can occur for specific Hamiltonians. Let us start by considering the case of noninteracting particles. We choose our mode basis to coincide with the single-particle energy eigenbasis, i.e. a † m,α creates a particle in a single-particle energy eigenstate with energy e m , m = 1, . . . M , such that the Hamiltonian reads In this basis, the Fock states |Ψ = |{N m,α } defined in Eq. (26) are many-body eigenstates with energy m N m e m . From this expression it is clear that these eigenstates are typically degenerate since the energy does not depend on the individual N m,α but only on the total populations N m = α N m,α . We can also look at these degeneracies of non-interacting systems in the basis {|λ, µ, ν } introduced in the previous section. If µ is a Young tableau of shape λ which contains N m times the index m associated with the single-particle energy e m , then |λ, µ, ν is also a many-particle eigenstate with energy m N m e m . Since this energy depends only on N , it can occur for several values of λ, µ and ν. In particular, since the single-row diagram λ = (N ) is compatible with any distribution N over the single-particle eigenstates, the spectrum of a given block λ is always contained in the one of N indistinguishable bosons.
Actually, the above results only require that the eigenstates of the many-body Hamiltonian are Fock states in a given mode basis. This is for example also the case in the infinite interaction limit of the BHH [U/J → ∞ in Eq. (46)], where Fock states in the site basis are many-particle eigenstates with energy U 2 M m=1 N m (N m − 1). In the intermediate case where the ratio U/J is finite, we find numerically that eigenstates belonging to different irreps λ are in general not degenerate.
For example, in Fig. 15, we plot the energy spectrum of N = 3 distinguishable particles in a three-site BH lattice, using the parametrization J = (1 − η)J 0 , U = ηJ 0 to cover the full range from no interaction (η = 0, U = 0, J = J 0 ) to infinitely strong interaction (η = 1, U = J 0 , J = 0). The 10 many-particle energies belonging to the bosonic sub-space λ = (3) are shown in blue, the 8 energies from the mixed-symmetry subspace λ = (2, 1) are colored in red (they are doubly degenerate for a system of distinguishable particles) and the single energy from the fermionic subspace λ = (1, 1, 1) is plotted in green (compare with Fig. 14). In accordance with the above discussion, degeneracies occur in the non-interacting limit U/J = 0 and for infinitely strong interaction U/J → ∞, otherwise the spectra of the various subspaces are different, except for accidental degeneracies.

C. Spectral signatures of partial distinguishability
We now apply our above analysis of the many-particle spectrum to decipher the frequency components that appear in the time evolution of EVs, for various levels of (in)distinguishability of the initial state.
For non-interacting particles, we have seen that the spectra associated with different irreps λ are all contained in the bosonic one, with λ = (N ). Therefore the set of frequencies contributing to the dynamics is independent of the decomposition of the initial state over the invariant subspaces and, consequently, also independent of its degree of (in)distinguishability. From the Heisenberg evolution of the annihilation operators associated with the single-particle eigenbasis, one easily concludes that the time-dependent EVs of kPOs can contain frequency components obtained as sums of k singleparticle frequencies ω mn = e m − e n . Note, however, that the amplitudes of the individual frequency components do depend on the state's DOI. This can for instance be seen in Fig. 7, where the time-dependent EVs of states with different DOI values are shown. All three curves share the same oscillation frequencies, only the amplitude of the oscillations differs. For systems with non-vanishing interaction strength, we learned in the previous section that subspaces associated with different irreps λ have a unique spectrum in general. The maximum number of frequencies coming from a subspace of dimension d = dim(λ, M ) is equal to d(d − 1)/2. Hence, the number of frequencies appearing in the time-dependent EV of a species-blind observable O increases with the number of subspaces involved in the decomposition (60) of the initial Fock state.
This can for instance be seen in the left panel of Fig. 16, where we plot the EV of n 2 1 [t] in a double-well system for three different initial states of N = 8 particles that interact with a strength U = 0.3J. As we go from the completely indistinguishable state (black curve) to the completely distinguishable one (red curve) by changing the species of the particles in the second mode, we observe both a decrease of the average value of n 2 1 [t] (consistently with the results in Sec. IV B) and a suppression of the oscillations around this average, which indicates the presence of more contributing frequencies.
A discrete Fourier transform (DFT) of the signals, shown in Fig. 16, confirms our analysis. In the given frequency inter- val, we find five frequencies that are common to the EVs of all states, five frequencies that are only common to the intermediate state (yellow curve) and the maximally distinguishable one (red curve), and one frequency that appears exclusively for the maximally distinguishable state. This is easily understood in terms of the irreps λ which appear in the decomposition (60) of the initial state: the second row of the Young diagram cannot be larger than the number of particles in the minority species. The fully indistinguishable state therefore only probes the symmetric irrep λ = (8), at the origin of the five frequencies common to all states. The two other states have finite weight in two additional subspaces, λ = (7, 1) and λ = (6, 2), from which the five frequencies that they share originate. The frequency exclusive to the maximally distinguishable state must come from the three dimensional subspace associated with λ = (5, 3), since the last subspace λ = (4, 4) is one-dimensional for M = 2 and does not contribute any frequency. Note that there are also higher frequencies outside the considered frequency window, not shown in the plot, to which the same rules apply.
In the above example, we clearly see a correlation between the number of frequencies contributing to the dynamics and the species distribution S: the more homogeneous the latter, the more different subspaces are probed, resulting in more frequencies. In this example, the DOI value I of the initial states (see legend of Fig. 16) also correlates with the homogeneity of species distribution and thus with the number of frequencies.
In general, however, the DOI does not only depend on the total species distribution, but also on how the particles of each species are initially distributed over the modes [60]. For example, the three states illustrated in the right panels of Fig. 16 all have the same DOI value I = 1/2, but different species distributions S. If we again plot n 2 1 [t] for each initial state, the EVs have the same time average, consistently with the DOI values. However, the DFT, shown in the right panel, reveals different sets of frequencies for each state. In accordance with our expectation, the EV of the state with S = (5, 3) (yellow curve) has one more frequency than the one with S = (6, 2) (black curve), which originates from the threedimensional subspace λ = (5, 3). Surprisingly, the state with a balanced species distribution S = (4, 4) (red curve) shows a smaller number of frequencies, although it could in principle probe the same subspaces as the other states (as well as the one-dimensional subspace λ = (4, 4), which does not contribute any additional frequency). However, by calculating the expansion (60) for this Fock state, we find that it has zero weight in the subspaces λ = (7, 1) and (5, 3), which we could not predict from the species distribution alone. Nevertheless, the frequency content allows to distinguish the three initial states with the same DOI value. Such spectral signatures therefore provide valuable complementary information on the initial state's distinguishability.

VI. CONCLUSION AND OUTLOOK
Many-particle interference is traditionally considered in the context of multiport optical interferometry, of which the Hong-Ou-Mandel experiment is the archetype. With this paper, we shifted the focus to the time-continuous, and possibly interacting, evolution of many-body systems, as encountered in condensed matter or cold atom physics. Accordingly, we looked for signatures of many-particle interference in the time-dependent expectation values of few-particle observables accessible to experiments on these platforms. In order to identify these signatures, we systematically varied the distinguishability of the system's constituents by assigning them to various mutually distinguishable species. Although these species are introduced as a convenient theoretical tool, real systems typically posses degrees of freedom which can act as distinguishing labels. In photonic interference, this can be the polarisation or the spectro-temporal form of the photons; in cold atom systems, the internal (electronic) state. These can therefore be exploited to experimentally observe the effects predicted in this paper, notably the sensitivity of density fluctuations to distinguishability, and the fact that single-particle observables are also affected by many-particle interference in interacting systems.
In particular, cold atoms trapped in optical lattices allow the realization of the multi-species Bose-Hubbard model which we used to illustrate our results, with hyperfine states playing the role of species. In these systems, high-resolution microscopes offer single-site resolution [70][71][72] as well as the possibility of counting the number of atoms per site up to N = 3 [73]. The controlled population of three different internal hyperfine states has also been recently demonstrated experimentally for trapped fermionic [74] and bosonic [75] atoms in a one-dimensional optical lattice. The species-blindness for the tunneling part of the Bose-Hubbard Hamiltonian is given by default, since the tunneling probability is to a good approximation independent of the internal hyperfine state, as shown e.g. in Refs. [74,75]. Although the interaction strength, i.e. the s-wave scattering length, does in general depend on the hyperfine states of the interacting atoms, measurements for 87 Rb reveal [76] that the inter-and intraspecies scattering lengths are very similar to each other. For other elements, mismatches of the scattering lengths can in principle be minimized with the help of Feshbach resonances.
Of course, understanding the effects of partial distinguishability in generic many-body systems, beyond the idealized scenario of initial Fock states and species-blind dynamics, remains an outstanding challenge. In particular, the states encountered in realistic systems can be much more complex than the Fock states considered here, where each particle is assigned to a given mode and species. We explained in Sec. III B under which circumstances our results can be easily extended to more general states, but we think that a study of manyparticle interference with states for which such an extension is not straightforward would be of particular interest. Indeed, this would shed light on the links between coherence, entanglement and indistinguishability, since such states display additional correlations that interplay with those induced by indistinguishability.
Another promising future research direction is to further explore the regime of strongly interacting particles. In our perturbative diagrammatic approach, observables are "dressed" by interaction vertices, such that many-particle evolutions involving more and more particles contribute to their expectation values. In the last sections, we used symmetry arguments to identify signatures of distinguishability that are independent of the interaction strength. It is however still unclear how these two approaches relate, and more theoretical efforts will be needed to further clarify the ways in which interactions and many-particle interference contribute together to the complexity of many-body dynamics.