Exploring new states of matter with a photonic emulator

We implement the equation of motion of the large-N Gross-Neveu model from strong interaction physics in photonic waveguide arrays and study one of its paradigmatic multi-fermion bound state solutions in an optical experiment. The present study constitutes an important first step towards waveguide-based simulations of phenomena relevant for high-energy physics.

In recent years, earlier speculations about the existence of new states of fermionic matter in the form of inhomogeneous phases where translational invariance is spontaneously broken have turned into a firm theoretical prediction for certain low-dimensional, exactly solvable field theories such as the Gross-Neveu model [1][2][3]. Similar phenomena have been predicted and extensively studied in a wide variety of research fields, ranging from condensed-matter systems (so-called FFLO phases in superconductors in large magnetic fields) [4,5], via ultracold atomic gases, nuclear physics, and the interior of neutron stars, to quark matter at highest densities; see [6][7][8] for reviews. While exact theoretical solutions are now available for 1D models much less is known for the relevant cases of 2D layered structures or in full 3D. This is not merely due to technicalities, but due to the fact that a deeper understanding of the Dirac equation in inhomogeneous phases is conceptually lacking; e.g. higherdimensional analogues of the 1D Peierls instability are still searched for [9]. Moreover, to date fermionic matter characterized by a spontaneously broken translational symmetry is extremely difficult to realize in an experiment. However, testing such phenomena experimentally would provide unprecedented insights into the very foundations of theory and trigger conceptually new theoretical approaches for the description of these systems.
Modern optics and photonics is driven by the fact that photons can be coherently controlled in space and time at the highest precision level. This goes hand in hand with recent developments in precision fabrication and design of optical systems. A prominent example is given by waveguide arrays [10] that can be designed as photonic analogues of systems governed by complex wave equations, including even quantum mechanical equations such as the relativistic Dirac equation. Appropriately designed photonic waveguide arrays with alternating refractive indices of adjacent lattice sites have already proved their capability to emulate a variety of relativistic phenomena in a * felix.karbstein@uni-jena.de wide range of parameter regimes, including Zitterbewegung [11], pair creation [12], particles with random mass [13], ultra-strong magnetic fields [14] and even tachyons [15]. Since the Dirac equation governs the dynamics of almost all known matter particles in the universe at the microscopic level, photonics has the potential to explore states of fermionic matter in an unprecedented way.
In this work, we present a photonic emulator for the physics of relativistic fermion systems and apply it to the massless Gross-Neveu model in the large-N limit [1]. Remarkably, at low temperatures T the latter favors a ground state where translational symmetry is spontaneously broken. This manifests itself in a spatially inhomogeneous scalar condensate, or equivalently a coordinate dependent fermion mass. As exact solutions are available for both the condensate shape and the full Dirac spectrum [2,3], one of those can be used as a paradigmatic example for mapping the Dirac equation in an inhomogeneous condensate onto photonic waveguide arrays.
The equation of motion of the fermion field ψ (i) follows from the Euler-Lagrange equation, yielding As the GN model is a relativistic QFT, and thus genuinely a multi-particle theory, the fermion fields are not single-particle wavefunctions but second quantized field operators featuring infinitely many positive and negative energy states. This implies that also the bilinear combination of the fermion fields inside the parentheses in Eq.
(2) is operator valued. In the 't Hooft limit, defined by sending N → ∞ while keeping N g 2 = const., this bilinear combination can be replaced by its expectation value in the considered multifermion state [16]. Corrections are parametrically suppressed by inverse powers of N , rendering this replacement exact. In this large-N limit, the equation of motion (2) reduces to with scalar potential S(x) given by Here, we focus on static configurations for which the expectation value ψ(n) ψ (n) is time-independent. Equation (3) corresponds to the Dirac equation describing fermions with a prescribed coordinate dependent mass S(x). The difficulty in solving it arises from the selfconsistency condition (4), that is, S(x) itself is defined in terms of solutions of the Dirac equation (3). Aiming at the ground state of the system for given fermion density and temperature, one has to solve Eqs. (3) and (4) with infinitely many single particle states occupied: apart from all negative energy states, additional energy levels are to be populated until the prescribed value for the fermion density is reached. The infinite sum over the negative energy states is divergent and requires regularization. However, it turns out that all dependences on the regularization, as well as the bare coupling constant g can be traded for the physical fermion mass m. The latter is the only parameter of the translationally invariant vacuum characterized by a filled Dirac sea and all positive energy levels empty [16]. The physical fermion density in the system is also measured relative to the vacuum state. Multi-fermion bound states of N f ≤ N fermions in excess to the vacuum are often referred to as baryons [17] due to their analogy to baryons in hadron physics. Baryon number one is assigned to a state with N f = N , and the fermion density in the system is conventionally parameterized by the baryon density ρ. The similarity of the model to the strong interactions has triggered a substantial amount of analytical and numerical studies in recent years [18][19][20][21][22][23][24][25][26][27].
In Fig. 1, the phase diagram of the GN model in the (ρ, T ) plane is shown [28]. In the regions where S(x) = const. and S(x) = 0 the ground state of the system is translationally invariant. In the gapped phase, an effective mass for the fermions is spontaneously generated by chiral symmetry breaking, while they remain massless in the gapless phase. Finally, in the phase where the scalar condensate S(x) exhibits an explicit dependence on x the ground state breaks translational symmetry and is characterized by a coordinate-dependent fermion mass. This is the regime of new states of matter [2,3,28]; cf. also [29,30] for applications in condensed matter physics.
It can be shown that the single-kink scalar potential corresponds to one of the cases for which Eqs. (3) and (4) can be solved analytically [31]. Equation (5)  The total fermion number associated with this object is N f = n−N/2 [32]. As it interpolates between asymptotic states with positive and negative physical fermion mass m for x → ±∞, respectively, it amounts to a manifestly relativistic object which has no non-relativistic limit [16]. Now, the idea of our work is to emulate the Dirac equation by a waveguide array implementing a spatially inhomogeneous potential S(x), which is a self-consistent solution of Eqs. (3) and (4). In this manner, multifermion bound-state formation, which is closely related to the phenomenon of translational symmetry breaking [7], can be directly probed. Such an experiment illustrates that photonic platforms have the potential to serve as a viable laboratory for the physics of relativistic selfinteracting fermion systems; in this way, a new and flexible tool to search for new states of matter in photonic experiments becomes available.
Equation (3) can be rewritten as FIG. 1. Phase diagram of the theory. New states of matter characterized by S(x) = const. arise for small temperatures [28].
We choose α = σ 1 and β = σ 3 , with Pauli matrices σ i , and discretize the spatial coordinate x on a lattice. Obviously, the off-diagonal matrix α couples the equations for the upper and lower spinor components. Following [33], we decompose a one-dimensional lattice with lattice spacing d into two independent sublattices of lattice spacing 2d to accommodate the two independent components of the spinor field: the upper component resides on the even sites, the lower component on the odd sites. In turn, the pair of lattice sites 2n and 2n − 1, with n ∈ Z, is associated with the same spatial coordinate x. When expressing the components of the spinor field as ψ 2n−1 (t) in terms of complex amplitudes a where we introduced the dimensionless time z = t/(2d) and the dimensionless potential S N = 2d S(N d). The evolution of the amplitude a (i) N with z is coupled to the neighboring amplitudes a (i) N ±1 . Equation (7) is amenable to a waveguide implementation [11]. As each of the distinct fermion species i ∈ {1, . . . , N } fulfills the same equation, it suffices to simulate it for a single fermion species. The coordinate z can be mapped onto the longitudinal coordinate of the waveguide array, and a (i) N to the amplitude of the light coupled into the N th waveguide. Here, it experiences the waveguide-specific index of refraction S N , as shown in Fig. 2(b). Using our conventions for the Dirac matrices, the analytical solution for the valence spinor associated with the single kink potential (5) is given by [16,31] where ϕ is an arbitrary global phase, which drops out in observables such as ψ † 0 ψ 0 . For our experiments, we fabricate various waveguide lattices using the direct-laser writing technology [35]. The length of the waveguides is l = 100 mm, the lattice constant is d = 16 µm, and the average refractive index of each full lattice is δn = 6 × 10 −4 . In order to implement the last term in Eq. (7), the two sublattices are realized by fabricating an alternating sequence of waveguides with high and low refractive index change [11]. This tuning is accomplished by varying the ratio between the writing velocities of adjacent waveguides in the sublattices; the magnitude of the writing velocity difference is proportional to the scalar potential S N . Importantly, changing the writing speed leaves the intersite hopping of κ = 0.14 mm −1 essentially unchanged [36]. A fluorescence microscopy technique [37] enables us to map the flow of light from the top of the sample and, thus, to visualize the spinor wave packet evolution. The array is excited by a broad Gaussian beam at a wavelength of λ = 633 nm with a spot size of ∼ 80 µm in the transverse direction, covering approximately 5 waveguides. A representative example of dynamics in such a fabricated waveguide structure is shown in Fig. 3. The refractive index detuning between the sublattices far away from the kink at the center of the structure -that is, at S N (±∞) -is ∆n ≈ 2 × 10 −4 . Upon excitation in the center of the array, i.e., directly at the kink, clearly a bound dynamics is observed in the form of an oscillatory motion of the light beam (see Fig. 3(a)) matching the analytical predictions for the valence level of the kink potential (5). Our observation is supported by a numerical integration of Eq. (7), that yields dynamics that resembles our experimental data (see Fig. 3(b)). The implemented (normalized) refractive index distribution of the waveguide lattice is shown in Fig. 3(c), emulating the scalar potential (5) in the center of the structure. A plot of the respective eigenvalues is shown in Fig. 3(d), showing the mid-gap state (8) that is localized at the kink.
Moreover, the eigenvalue diagram 3(d) shows that not only the mid-gap bound state as predicted for the GN valence state is emerging, but that also the two states at the inner edges of the bands slightly penetrate into  gap and, hence, become localized bound states as well. This behavior is a result of the intrinsic lattice discretization and, hence, a characteristic feature of our photonic emulator. When launching light into the center of the lattice, all three states are excited; as a result, one observes a beating between the states. However, this beating is spatially confined in the very vicinity of the kink. Interestingly, the beating length l B can be used to estimate the width of the band gap in the experimental system, as it is connected to the energy difference ∆E of the states by l B = 2π/∆E [38]. As the other localized states reside close to the edge of the band gap, ∆E is a measure for the band gap width. To explore this feature, we conducted several experiments similar to that in Fig. 3 with different normalized refractive index detuning ∆n at the edges of the waveguide lattices (i.e., at S N (±∞)). In Fig. 4, we plot the experimentally determined gap widths as a function of ∆n. Evidently, the band gap width increases almost linearly for increasing detuning. Equation (5) implies that the width of the band gap is directly related to the fermion mass m, which also determines the spatial localization of the valence spinor (8). Hence, our emulator allows for an indirect determination of both the fermion mass and the spatial localization of the valence spinor (see the insets in Fig. 4), via observing the oscillation period of the light wave.
In conclusion, we reported an implementation of a photonic emulator for phenomena related to bound-state formation and translational symmetry breaking in the ground state of the large-N Gross-Neveu model. Our results suggest that waveguide optics could provide an experimentally accessible classical emulator to test complex predictions concerning the spontaneous breaking of translational invariance and the formation of spatially inhomogeneous phases in the phase diagram of the theory. A particularly interesting future application of our photonic emulator will be the study of multi-fermion bound state formation in quantum field theories in the non-relativistic limit. Here, the Dirac sea involving infinitely many filled negative-energy states can be integrated out. This results in a no-sea effective field theory featuring fermion fields of manifestly positive energy only [39], for the study of which our photonic emulator should be ideally suited. The no-sea analogue of the Gross-Neveu model as well as many other 1+1 dimensional field-theories with four-fermion interactions can even be solved analytically which provides various benchmark solutions for such studies [40]. In particular for phenomena in 2D and 3D settings, where the search for new states of matter is extremely challenging, photonics may provide a new promising route for exploring these phenomena experimentally.
The authors would like to thank C. Otto for preparing the high-quality fused silica samples employed in this work. This work has been carried out within the framework of the ACP Explore project "Enlightning New States of Matter" of the Abbe Center of Photonics (ACP). AS acknowledges financial support from the European Research Council (grant EPIQUS), the Al-  Fig. 2 (b).