Probing Atomic Majorana Fermions in Optical Lattices

We introduce a one-dimensional system of fermionic atoms in an optical lattice whose phase diagram includes topological states of different symmetry classes. These states can be identified by their zero-energy edge modes which are Majorana fermions. We propose several universal methods of detecting the Majorana edge states, based on their genuine features: zero-energy, localized character of the wave functions, and induced non-local fermionic correlations.


Introduction
There is at present growing interest in topological phases of matter -fermionic insulators and superconductors with a gapped excitation spectrum and nonlocal topological order (TO). TO is characterized by a topological invariant taking discrete values [1,2,3,4,5,6]. An immediate consequence of TO is the existence of robust zeroenergy modes localized at defects, edges, and interfaces between different topological phases [7,8,9]. The number of these states and their properties are directly linked to TO and, therefore, can be viewed as a probe of the topological state providing access to the corresponding phase diagram. A striking example are zero-energy Majorana fermions, as discussed in the context of topological quantum computing [10,11,12].
Most of the developments and proposals to realize and study topological phases in general, and Majorana fermions in particular, have been related to condensed matter systems [13,14,15,16,17], with a first observation of Majorana edge modes in a hybrid superconductor-semiconductor nanowire device recently reported in Ref. [18]. On the other hand, there are several promising proposals to realize topological phases with quantum degenerate gases of atoms and molecules [19,20,21,22,23,24], with questions of detection of such phases presently in the focus of interest [25,29,30,31]. Below we will outline a measurement scenario for topological phases and phase transitions by monitoring signatures of Majorana fermion edge-states in atomic systems. Our analysis builds directly on recent experimental advances such as single site addressing and measurements in optical lattices [26,27] in combination with traditional time of flight and spectroscopic techniques, and complements the recent proposals to detect topology in the bulk as proposed in [29,30,31].
To illustrate our ideas for detection we introduce a simple, but experimentally realistic example of a zig-zag chain (c.f. Fig. 1a). This model is an extension of the familiar Kitaev model of spin-less fermions coupled to a BCS-reservoir with the additional feature of next-to-nearest neighbor couplings (see below), and can be realized with cold atoms generalizing ideas outlined in Ref. [23]. The model has a remarkably rich phase diagram, allowing different topological states in two symmetry classes supporting Majorana edge modes, and provides an ideal playground to demonstrate the presence of TO and the topological phase transitions via related Majorana edge states. We propose a fully reproducible preparation of an initial state for Majorana fermions and detection with atomic measurement techniques, including (i) zero-energy, (ii) localization near the edge and (iii) induced non-local fermionic correlations, which, together, allow for an unambiguous probe of Majorana fermions. While the present example is 1D, the techniques described below have a straightforward extension to 2D, as illustrated by discussion of a p x + ip y superfluid.

The Zig-zag Chain
In the following Section we discuss the Zig-zag Chain, a one-dimensional (1D) model system of spinless fermions which allows for a variety of topological phases, thus providing an ideal playground for testing our detection techniques. We start in Subsec. 2.1 introducing the model and presenting its topological phase diagram. The reader interested in the derivation of these results is referred to Subsec. 2.2 . Further, we review the ideas of [23] in Subsec. 2.3, explaining how the zig-zag chain can be realized in a cold atom implementation.

The Model System and its Properties
We consider single-component fermions on a finite 1D chain of size L with Hamiltonian H = H 1 + H 2 + H µ (c.f. Fig. 1a). Here, H µ = −µ L j=1 a † j a j , and where a i and a † i are fermionic operators (i = 1, . . . , L), J α ≥ 0 and ∆ α = |∆ α | e iφα are nearest-neigbour (NN, α = 1) and next-to-nearest-neighbour (NNN, α = 2) hopping and pairing amplitudes, respectively, and µ is the chemical potential. As shown in [23] and reviewed in Subsec. 2.3, the pairing terms in H α are obtained by a Raman induced dissociation of Cooper pairs (or Feshbach molecules) forming an atomic BCS reservoir, with µ a Raman detuning. In Fig. 1a we represent the chain as a zig-zag, which leads naturally to a cold atom implementation with optical lattices allowing control of the relative strength of NN and NNN amplitudes by changing the zig-zag geometry.
Although H can be viewed as a Hamiltonian of two coupled Kitaev chains [32] (with odd or even sites), the resulting topological phase diagram is substantially richer. Since the coupling parameters can be changed in a real-time experiment, this setup gives rise to a platform for exploring topological properties of matter, including topological phase transitions and Majorana fermions. We summarize here only the main features of the topological phase diagram, and refer the interested reader to Subsec. 2.2 a) The topological symmetry class [1] of H depends crucially on the relative phase φ = φ 1 − φ 2 of the pairing amplitudes (one of the phases, say, φ 2 , can be gauged away by redefining the operators a j ). For φ = 0, π (mod 2π) it is the class BDI (Cartan) with time-reversal, particle-hole, and chiral symmetries, for all other values of φ it is class D with only particle-hole symmetry. Both classes allow topologically nontrivial states characterized by the integer-valued winding number ν ∈ Z in the class BDI and by the Chern parity number P C = 0, 1 ∈ Z 2 in the class D. The topological phase diagram contains states with ν = 0, ±1, ±2 for the BDI class, and with P C = 0, 1 for the D class. Fig. 2a summarizes results for the simplest case µ = 0 and |∆ α | = J α , where the phase diagram is characterized by the relative phase φ and the ratio J 2 /J 1 . We find that ν = 2, P C = 0 for J 2 > J 1 , else one has P C = 1 and ν = +1 or −1 for φ = 0 and π, respectively. Therefore, for both symmetry classes, the point J 1 = J 2 corresponds to a topological phase transition, where the bulk excitations become gapless (see Subsec. 2.2).
The presence of a nontrivial TO in the bulk leads to zero-energy modes, which in our case are Majorana fermions, at the interfaces between states with different TO. The number of these states is related to the difference of the corresponding topological invariants [8,9]. Using the hermitian Majorana representation, where c 2j−1 = a † j + a j and c 2j = (−i)(a † j − a j ) obey {c j , c l } = 2δ jl , the quadratic Hamiltonian reads H = i j,l A jl c j c l with real A = −A T . Diagonalizing A for a finite chain we identify the zero-energy modes and associated Majorana edge-mode operator γ E = k v Ek c k with (real) coefficients v Ek localized at the interface with some localization length l loc . [To be precise, the energy of these modes ∼ exp(−L/l loc ) approaches zero for L → ∞.] Thus, by detecting Majorana zero-energy edge modes and their number we can access the full topological phase diagram of our model: For the BDI class we have in total two such modes for J 1 > J 2 : γ L = c 1 and γ R = c 2L for the left and the right edge, respectively, and four modes for J 1 < J 2 . The two additional zero modes are exponentially decaying inside the bulk (see Fig. 2b and Subsec. 2.2). For the D class there are two zero-energy modes for J 1 > J 2 , which decay exponentially inside the chain, and no such modes for J 1 < J 2 (Fig. 2c).
The presence of Majorana zero-energy modes comes along with a degeneracy of the ground state: The dimension of the ground-state subspace is 2 N M /2 , where N M is the number of the Majorana modes (always even in our case). A state in this subspace generically has correlations between Majorana operators from different edges, resulting in non-local fermionic correlations (in contrast to local correlations in the bulk). For the case with only two Majorana modes (ν = 1 or P C = 1), the two ground states |G ± have even (+) or odd (−) number of fermions, respectively, and the nonlocal correlations G ± | γ L γ R |G ± = ±i are related to the fermionic parity, in full analogy to Ref. [32].

Derivation of the Topological Phase Diagram
In the following we present a detailed derivation of the topological phase diagram presented in Subsec. 2.1. According to the general tenfold classification scheme [1,3], the topological class to which a given Hamiltonian belongs to, is determined by the invariance properties of the Hamiltonian under time-reversal, particle-hole (or charge conjugation), and chiral (or sublattice) symmetry. This class specifies then, whether the states with a nontrivial topological order could exist and provides the corresponding topological invariant to distinguish them.
For a translationally invariant 1D spinless Bogoliubov-de Gennes lattice Hamiltonian in the quasimomentum basis where with a k = L −1/2 k e ikj a j , the invariances are equivalent to the following conditions: for the chiral symmetries, respectively, where U T , U C , and Σ are unitary matrices satisfying U T U * T = ±1, U C U * C = ±1, and Σ 2 = 1. For the Hamiltonian H of Eq. (1) one has ξ k = −(J 1 cos k + J 2 cos 2k) − µ/2 and ∆ k = −i(∆ 1 sin k + ∆ 2 sin 2k), such that the matrix H k is  . Spatial distribution of the zero energy modes of H for J α = |∆ α | and φ = 0. a) If J 2 < J 1 we have two zero energy modes that are located at the boundary. b) If J 2 > J 1 , two additional zero energy modes with an exponential decay appear. c) These modes extend more and more over the lattice when we approach the transition point For generic ∆ 1 = |∆ 1 | e iφ 1 and ∆ 2 = |∆ 2 | e iφ 2 , only the particle-hole symmetry condition is fulfilled with U C = σ x (as it should be for a general Bogoliubov-de Gennes Hamiltonian), and, therefore, the Hamiltonian H belongs to the class D. For , the conditions for the time-reversal and chiral symmetries are also satisfied with U T = diag(e −iφ 2 , e iφ 2 ) and Σ = σ x U T , and the Hamiltonian H belongs to the chiral BDI class. Note that the phase φ 2 can be gauged away by a k → e −iφ 2 /2 a k such that Re(∆ k ) = 0 and the vector h k belongs to the yz-plane for all k. Geometrically speaking, for the BDI class the vectors h k for all k are in the same plane, which is the yz-plane for φ 2 = 0.
For the chiral BDI class in 1D, different topological states are classified by the  Figure 4. Topological phase diagram in the chiral class for the case |∆ 1 | = J 1 but still ∆ 2 = J 2 and µ = 0. We see that we can obtain regions with ν = 0 (white), 1 (striped) and 2 (grey).
integer valued winding number ν ∈ Z defined in the following way (see, for example, [7]): By gauging away the phase φ 2 and using a unitary transformation U = exp(iπσ y /4), the Hamiltonian can be transformed into a canonical form with q k = ξ k + ∆ k and ϕ q k being the phase of q k . Then Alternatively, the winding number can be defined using the unit vector n k = h k / h k in the yz-plane (φ 2 is gauged away) as number of times the vector n k wraps around the unit circle when k goes around the Brioullin zone (hereê x is the unit vector along the x-axis). Note that as long as h k = 0 (or, equivalently, |q k | = 0), which corresponds to the gapped spectrum of the Hamiltonian, the winding number is well-defined and, therefore, can change its value only when the gap closes. This gap-closing condition is the necessary one for the topological phase transition. For |∆ α | = J α and µ = 0, the gap closing condition corresponds to J 1 = J 2 , and straightforward calculations give ν φ=0 = 1 two out of four zero-energy modes start to proliferate inside the bulk of the chain and become gapped on the other side of the transition; see Fig. 3c. The topological phase diagram in this chiral class for the case |∆ 1 | = J 1 , but still ∆ 2 = J 2 and µ = 0, is shown in Fig. 4. In this case h k = (0, ∆ 1 sin k+J 2 sin 2k, J 1 cos k+ J 2 cos 2k), where we assume φ 2 = 0, and, hence, real ∆ 1 . The boundaries between different topological phases results from the condition h k = 0.
For the D class [when φ = 0, π(mod2π)], there are only two different classes of the topological states characterized by the Chern parity number P C ∈ Z 2 . To define this topological invariant [4], one has to extend the Hamiltonian H k and thus h k from a onedimensional Brillouin zone k ∈ [−π, π] ∼ S 1 (topologically equivalent to the circle S 1 ) to a two-dimensional (2D) Brillouin zone k, t ∈ [−π, π] ∼ S 1 × S 1 = T 2 (topologically equivalent to a two-dimensional torus T 2 ), H k → H k,t = h k,t · σ, such that the resulting Hamiltonian is gapped and belongs to the D class in 2D: The unit vector n k,t = h k,t /| h k,t | maps then the 2D Brillouin zone T 2 into a 2D sphere S 2 . Different topological classes of such mappings are distinguished by the integer-valued Chern number It turns out [4], that the parity P C of C does not depend on the chosen extension of the Hamiltonian, and, therefore, provides the topological invariant for 1D Hamiltonians in the D class. (We set P C = 1 for odd C and P C = 0 for even C. Topologically nontrivial states correspond to P C = 1.) The extension can be constructed using a hermitian matrixh(k, t) = i h i (k, t)σ i , which has a gapped spectrum and depends on two parameters k ∈ [−π, π] and t ∈ [0, π] laser BEC reservoir optical lattice laser RF pulse Figure 6. Creation of the pairing term ∆a † j a † j+1 + h.c. according to Ref [23]. A onedimensional system of trapped fermions in two spin states is coupled to a molecular BEC via a RF pulse. An additional Raman laser with strong driving projects out one of the spin components, leading to the desired pairing term.
such thath(k, 0) = H k andh(k, π) = | h k |σ z . Note thath(k, t) interpolates between the initial Hamiltonian H k and the "trivial" one | h k |σ z , both from the D class (for t = 0 or π,h(k, t) does not necessarily belong to the D class in 1D). Such a matrix always exists because it belongs to the A class of general hermitian matrices (with no symmetries) with a gapped spectrum that is topologically trivial. Then the matrix is in the D class and provides the desired extension. The interpolationh(k, t) can be obtained, for example, in the following way: h i (k, t) for t ∈ [0, π/2] corresponds to a rotation of h k = h i (k, 0) to the y-axes, such that h i (k, π/2) is parallel to the y-axes for all k, and thenh(k, t) for t ∈ [π/2, π] describes the rotation of h i (k, π/2) to the z-axis. Following this strategy for the case |∆ α | = J α with µ = 0 we find C = 0 for J 2 > J 1 and C = sign(sin φ) for J 1 > J 2 , where φ = φ 1 − φ 2 . The point J 1 = J 2 corresponds to the topological phase transition, at which the gap vanishes, see Fig. 2c.

AMO realization of the Zig-zag Chain
The zig-zag chain defined in Eq. (1) allows for an optical lattice realization, as we briefly explain in this Subsection (see also Fig. 6). While the hopping term a † j a j+α + h.c. arises naturally in an optical lattice setup, the pairing term a † j a † j+α + h.c. can be engineered via the coupling of the system to a BEC reservoir of Feshbach molecules, as explained in [23]. The main idea is to couple the two internal spin states of the trapped 1D system of fermions to a Feshbach molecule via an RF pulse. If we denote by (a † p,↑ , a † p,↓ ) the two internal states of the trapped atoms with momentum p, then the result of the RF pulse is an effective pairing term of the form ∆a † p,↑ a † −p,↓ + h.c. Driving the lattice fermions additionally with a Raman laser creates an effective magnetic field, which for strong b) 0 10 20 0 driving (i.e. large magnetic fields) projects out one of the spin components, such that we obtain the spinless pairing term of Eq. (1) [23].

Preparation of Majorana modes
The main difficulty in detecting the Majorana modes is that their signal comes only from the edges and, therefore, has poor signal-to-noise ratio. To overcome this, an experiment has to be performed several times thus requiring a reproducible preparation of a desired initial quantum state of an open wire. In our setup, this can be achieved by starting with a closed quantum chain in a fully paired [34,35] (and, hence, unique) state corresponding to the ground state of a Hamiltonian (see Fig. 1b)). This state has an even parity (all particles are paired) and only local correlations between the sites. Making use of local addressing [26,27] we can "cut" the chain by ramping up the chemical potential on the site L, This process affects only the four Majorana modes c 2L−2 , c 2L−1 , c 2L , c 1 , and the corresponding instantaneous eigenvalues are 0 and 2µ for the odd parity sector, and µ ± 4 + µ 2 for the even one, such that there are two degenerate ground states of the Kitaev chain for µ → ∞. Because a superposition of fermionic states with different parity is forbidden by superselection rules, the corresponding eigenstates never mix during the adiabatic ramp (cf. Fig. 7a), and, hence, the final state of the open chain has even parity. In Fig. 7b we depict the evolution of the Majorana correlation functions between nearest-( ic 1 c 2L , solid) and next-to-nearest ( ic 1 c 2L−1 , dashed) neighbors, as well as the occupation of the site L ( n L , dotted) during the adiabatic ramp.
For the "zig-zag" chain we can start with two decoupled closed chains (with even and odd lattice sites) when only ∆ 2 and J 2 are nonzero but ∆ 1 = J 1 = 0. We then cut The presence of the Majorana modes (solid) leads to oscillations of the signal S(x) (see main text) which are absent in the bulk (dashed). Having their origin in nonlocal fermionic correlations of modes located at the edges, the oscillation frequency is independent of the size of the shade (c) which allows to rule out a finite size effect encountered, e.g. in a system of freely hopping particles (d) (see main text). them as described above into two open chains each with two edge Majorana modes and in the even fermionic parity state with all fermions being paired. Now we can switch J 1 and ∆ 1 on adiabatically to create a single quantum chain with even parity. Note that during this process there is no single-particle transition between the two subchains (with even and odd sites) because such a transfer would create single-particle excitations in both subchains which cannot occur in an adiabatic process due to the conservation of energy. Thus, each of the subchains remains in the even parity state.
Further increase of NN amplitudes results in a topological phase transition with closing the single-particle excitation gap allowing for a single-particle exchange between the subchains. Then, only the (even) parity of the entire chain is preserved. We would like to add here that parity violating processes, such as three-body losses, are sufficiently weak in fermionic systems (0.1-1s) in order not to be an experimental obstacle for the realization of Majorana physics Note that a similar idea can be used to create the initial Majorana state in the dissipative setting by adiabatically emptying one site [24]. Note that the cutting procedure avoids the need to establish the edge-edge correlations via a sequential process, which takes polynomial time due to the Lieb-Robinson bound [36,37].

Time-of-Flight (TOF) imaging
Having prepared the atomic Majorana fermions following the protocol presented in Sec. 3, we present now two possibilities to detect them in AMO setups. The first approach that is discussed in this Section is based on Time-of-Flight imaging allows to detect the existence and number of Majorana fermions in an optical lattice setup. Further, we present a complementary approach based on a spectroscopic setup in Sec. 5. Non-local correlations can be detected in TOF imaging, as illustrated in Fig. 8a. We consider atoms of mass m released at time t = 0 from lattice sites R j = ja with lattice spacing a. At time t ma 2 / (far field approximation) the atomic density distribution at the detector is given by n(x) ∼ j,j e 2ix(R j −R j )m/( t) a † j a j revealing the initial correlations of the atoms in the lattice. Therefore, the long-range Majorana correlations R j − R j ≈ L M (distance between the edge modes) will result in rapid oscillations of n(x) as compared to slow oscillations originating from short-range bulk correlations (R j − R j ≈ a). The contrast of the Majorana signal can be enhanced by local addressing: we shade all but N b sites adjacent to each of the two edges (N b = 2 in Fig. 8a), thus measuring n(x) = 2N b S(x), where we have normalized by the number of sites contributing to the signal, 2N b .
The calculation of the density distribution at the detector in a TOF experiment requires the knowledge of the first moments a † j a j of the state in the lattice before the trapping potential is switched off. For J α = ∆ α , J 2 = µ = 0 (ideal Kitaev chain), the non-vanishing first moments in the Majorana language are c 2j c 2k ± = c 2j+1 c 2k+1 ± = δ kj , c 2j c 2k+1 ± = − c 2j+1 c 2k ± = −iδ kj , and the edge-edge correlations are given by c 2L c 1 ± = im ± (d) [m(d) = 0 if we consider the bulk only] with d being the distance between the two Majorana edge modes. When only N b sites on the right and on the left sides of the chain contribute to the signal (other sites are covered with a shade), the atom density distribution at the detector reads where N bulk (x) = 2N b sin 2 ( mxa t ) + cos( mxa t ) stems from the atoms in the bulk only and the two edge sites separated by L M = ad give rise to N edge (x) = cos( mxa t d)/2 multiplied by the Majorana correlation m(d) = −i γ L γ R . For the general case, we obtain the correlation matrix a † j a j via a numerical diagonalization of the Hamiltonian H. We see that the second rapidly oscillating term reflects the presence of the longrange Majorana correlations and allows to determine the occupation of the Majorana subspace. If the initial state of the chain is prepared as described above, then m(d) = −i G + | γ L γ R |G + = 1, and the amplitude of this terms reaches its maximum. The results of the calculations for the ideal chain are shown in Fig. 8b: The presence of the Majorana modes leads to oscillations that are absent in the bulk. Having their roots in long-range fermionic correlations, we can distinguish these oscillations from those resulting from a finite size sample by changing the size of the shade: A signal indicating the presence of Majorana fermions exhibits oscillations with the same frequency but reduced amplitude when we include more sites of the bulk (see Fig. 8c). On the contrary, oscillations induced by a finite size effect have in general a frequency that is sensitive to the shape of the shade (see Fig. 8d where we consider a system of freely hopping fermions).
By ramping the chemical potential, we can detect the location of the phase transition, since the oscillations disappear when crossing the transition point at µ = 2 (see Fig. 9a). Further, from the amplitude A of the signal S(x) we can deduce the number of Majorana modes, as depicted in Fig. 9b: The amplitude for four zero modes (J 2 > J 1 ) drops by a factor of two upon reaching the transition point J 2 = J 1 .
Finally, we also provide some results for experimentally realistic situations like a non-ideal chain and the influence of an external confinement of the atoms. In Fig.  10 a) we show TOF signals for the non-ideal case where J = ∆, proving that the oscillations will still be present in such a scenario. Further, we show the effect of an external harmonic trapping potential modeled by x a x for a system of L = 30 sites and ∆ = 1.2 in Fig. 10 b). In the inset we present the local density distribution of the atoms in the trap. We find that in the case of not too large inhomogeneities the TOF experiment can be used for the detection of Majorana fermions.

Spectroscopy
While TOF imaging allows to detect the existence and number of Majorana fermions in an AMO setup, the combination of a spectroscopic setup with TOF allows to measure not only the energy of the Majorana states but also their wave functions. (The use of spectroscopic RF absorption spectra to detect Majorana modes in 2D p-wave superconductors was considered in Ref. [25].) Consider a 1D lattice with two bands where We then couple upper and lower bands on the first and the last sites of an open chain via a time-dependent perturbation Fig. 11a) and measure the momentum distribution in the upper band as a function of the frequency Ω, e.g., via TOF imaging.
We can now rewrite the external perturbation in terms of quasiparticle operators a l and Majorana operators γ L , γ R using a 1 = (γ L +ã † 1 −ã 1 )/2, and a L = (ã † L−1 + a L−1 − iγ R )/2. Choosing the even parity ground state |G + as an initial state and using γ L |G ± = |G ∓ and γ R |G ± = ±i|G ∓ , a straightforward application of the standard time-dependent perturbation theory results in the following the upper band momentum distribution: where we have introduced the notation A 0 = [m + (d) − m − (d)] sin 2 k 2 (L − 1) + m − (d) and f t (x) = 4 sin 2 (xt/2 )/x 2 . Further, ε k = ∆E b + 2j cos k is the dispersion in the upper band and m ± (d) denote the initial occupation of the even/odd parity subspace in the lower band. The contribution n 0 corresponds to a parity flip in the lower band (no energy cost) and the creation of an excitation (particle) in the upper band (energy ε k ). Thus, the absorption peak at Ω = ε k in n 0 indicates the ground state degeneracy. The edge-edge Majorana correlation length is encoded in the oscillation period. The term n 1 results from the creation of excitations in both upper (energy ε k ) and lower (energy 2J) band. The corresponding absorption peak is located at Ω = 2J + ε k , providing the direct measurement of the pairing energy gap (2J here) from the distance between the two peaks. Thus, the topologically non-trivial phase leads to a clear signal at Ω = 2J (for µ = 0 one has to replace Ω → Ω − µ, cf. Fig. 11b). Note that this setup allows to detect the presence of the zero energy Majorana subspace independently of its purity. Further, the number of Majorana modes is encoded in the strength of the resonance signal for a fixed momentum k. In Fig. 11c we present the results for the Hamiltonian H with J α = |∆ α |, φ = 0, where the perturbation is on sites 1, 2, L − 1 and L. Starting from a state with four Majorana modes, we reduce the ratio J 2 /J 1 . At the transition point, the amplitude of the resonance signal n k=π/2 drops to half its magnitude indicating the disappearance of two of the Majorana modes. If we go away from the ideal case, the quasiparticle operators and the corresponding matrix elements can be determined numerically.
With a slight modification of the above setup one can also probe the localized character of the Majorana zero modes. Namely, we apply the perturbation to the first n L sites: V in the upper band is now where v jl are the coefficients of the Bogoliubov transformation to the quasiparticle operatorsã ν (ν labels the quasiparticle modes with energies E ν ) in the lower band, , which diagonalizes H. Let us now consider k = 0 and the external frequency being in resonance with the lowest energy (ν = 0 Majorana) level, Ω = Ω r = ∆E b + k=0 . In this case, only the Majorana mode is excited in the lower band, and one has n k=0 (Ω r , t) ∼ n L j=1 v j0 2 ≡ M(n L ). The localized character of the Majorana mode can now be established by looking at the dependence of n k=0 (Ω r , t) on the number of sites n L affected by the perturbation. To be more specific, the absolute change of M(n L ) when n L is increased by one, ∆M(n L ) = |M(n L + 1) − M(n L )|, is expected to be ∆M(n L > 1) = 0 for the ideal Kitaev chain and ∆M(n L ) ∼ exp(−n L /l loc ) for a generic one. This is reflected in a comparison of the normalized quantity ∆M(n L )/∆M(1) (red circles) with the normalized wave function of the Majorana mode, |v j0 | 2 / |v 20 | 2 (dashed) as shown in Fig. 12.

Probing Majorana fermions in 2D.
In the last Sections we have presented methods to prepare and detect Majorana fermions in a one-dimensional model system. Note, however, that the proposed detection methods are universal and can be applied to other systems. As an illustration, we consider a 2D p x + ip y superfluid, as it was proposed in [33,38] with Hamiltonian H S = H hopp +H pair +H µ , with a NN hopping H hopp , a pairing term H pair = x,y ∆ x,y a † x a † y +h.c. where ∆ x,y = ∆ x (δ y,x+e 1 +iδ y,x+e 2 ), and H µ = −µ x a † x a x . The fermionic operators a x are defined on a 2D (square) lattice and the vectors e 1,2 are two independent elementary lattice translations. We introduce two vortices (each with vorticity 1) via a position  dependent order parameter ∆ x = |∆|e iφ(x) , which can be written onto the superfluid by phase imprinting in the Raman process [26].
For |µ| 4J we have a topological ground state supporting two Majorana modes located in the vicinity of the two vortices (cf. Fig. 13 a). The comparison of the TOF signal with and without vortices (Fig. 13c) shows that the presence of the Majorana modes leads to an oscillatory behavior of the density n(x) at the detector. Again, the contrast of the signal is enhanced by shading all but a small region near the vortices (Fig. 13b). To detect the Majorana modes in the spectroscopic setup, we apply an external perturbation V (t) = V 0 x∈Λ a x b † x e −iΩt + h.c., where Λ includes the lattice sites inside a w × w small regions around the vortices (here: w = 3). The resulting momentum distribution in the upper band n kx,ky for k x,y = π 32 is depicted in Fig. 13d. The presence of the Majorana modes leads to a clear signal at Ω /J = ε π 32 , π 32 /J ≈ 1.19 that is absent in a vortex-free setting.

Summary
In this work we have presented various techniques that allow for an unambiguous detection of atomic Majorana fermions and topological order with standard quantum optical tools, as time-of-flight imaging and spectroscopic techniques within our framework. To this end, we have introduced a one-dimensional model system that allows for an AMO realization and provides, due to its rich topological phase diagram, and ideal playground to test our ideas. Since some of our detection schemes require the preparation of a ground state with a definite parity, we have further provided a protocol that allows to achieve this goal. In addition, our detection techniques are universal in the sense that they can be used to detect atomic Majorana fermions in any dimension and geometry, as we have illustrated by the example of a 2D p x + ip y superfluid.