Non-Hermitian physics is a rapidly developing field combining interesting fundamental concepts with important potential applications [1,2]. The research is particularly well advanced in photonics [3], thanks to the well-developed optical techniques of controlling the effective (both real and imaginary) potentials and of the measurements of the wave functions. The key step was the discovery of the parity-time (PT)-symmetry-breaking transition for the eigenstates occurring at exceptional points of non-Hermitian Hamiltonians [4]. Such a transition makes the system chiral, allowing the nonreciprocal optical transport without external magnetic field [5,6] useful for realizing optical isolation [7,8]. The PT-symmetric transition occurs at an exceptional point, where the matrix representing the non-Hermitian Hamiltonian becomes defective, and its eigenvectors coalesce. The most well-known example of such a transition is the damped harmonic oscillator, whose frequency becomes strictly zero above critical damping (overdamped regime).

After several initial attempts [9–11], a lattice with a non-Hermitian PT-symmetric potential was considered in a seminal work [12] that has demonstrated the modification of the dispersion of the lattice modes, with the formation of nondispersive (flat) regions bounded by exceptional points, for a sufficiently large degree of non-Hermiticity. This work has spawned an important activity focused on the studies of various lattices in both one- (1D) and two-dimensional (2D) cases, mostly within the tight-binding approximation (with non-Hermitian on-site potentials or coupling terms), in photonics [13] and other systems [14]. These studies often analyzed the interplay of non-Hermiticity with topology [15–18] or interactions [19,20]. In those works, the consideration was usually restricted to a single band (or a pair of branches, as in the case of the Dirac Hamiltonian), and thus usually a single region with PT-symmetry breaking was observed. However, it turns out that such a restricted consideration is not always correct, so going beyond the tight-binding limit (as in the original case of Ref. [12]) or considering several branches becomes sometimes necessary.

One of the tools for the experimental studies of periodic potentials is the Talbot effect [21,22], known as a self-imaging process. To date, the non-Hermitian Talbot effect was considered only theoretically in 1D photonic lattices [23–25]. As to 2D photonic lattices, such as hexagonal lattices, the Talbot effect was mostly accomplished based on the real part of the index [26,27], that is, in a Hermitian configuration.

In this work, we experimentally demonstrate the dynamic behaviors of a wave packet in a non-Hermitian photonic graphene structure realized in a three-level $\mathrm{\Lambda }$-type atomic vapor configuration. The observed PT-symmetry-like phase transition in the current 2D honeycomb optical lattice is achieved with uniform imaginary part on lattice sites, which is extremely different from previously demonstrated non-Hermitian dynamic behaviors in 1D optical lattices relying on gain and loss in atomic systems [28,29]. Also, we observed the Talbot effect of the discretized probe field at certain propagation distances outside the cell. Different from previous works, here we demonstrate the Talbot effect of a 2D hexagonal lattice resulting from a non-Hermitian configuration with both $n_R$ and $n_I$ modulated into honeycomb profiles. We show that the non-Hermitian character of the lattice allows switching the lattice symmetry between the honeycomb and simple hexagonal (“inverted-honeycomb”) lattices. This occurs thanks to the PT-symmetry breaking in the $\mathrm{p}$-band, giving rise to the lowest decay of the corresponding modes. Thus, the degree of non-Hermiticity of the established complex honeycomb photonic lattice becomes a new “control knob.” This can be promisingly exploited to develop an all-optical-controllable Talbot–Lau interferometer [30] and improve the related applications, such as imaging cold atoms [31,32] and atom lithography [33].

Figure 1 shows the scheme for experimentally demonstrating the evolution dynamics of the probe wave packet passing through the laser-induced complex photonic lattice. The probe field ${\mathit{E}}_1$ (wavelength ${\lambda }_1=780.2\mathrm{nm}$, frequency ${\omega }_1$, Rabi frequency ${\mathrm{\Omega }}_1$) and the coupling field ${\mathit{E}}_c$ with a spatially periodic intensity profile propagate along the same $z$ direction to drive a three-level $\mathrm{\Lambda }$-type ${}^{85}\text{Rb}$ atomic system [Fig. 1(b)]. Three coupling beams ${\mathit{E}}_2$, ${\mathit{E}}_2^{\prime }$, and ${\mathit{E}}_2^{″}$ (${\lambda }_2=780.2\mathrm{nm}$, ${\omega }_2$, ${\mathrm{\Omega }}_2$, ${\mathrm{\Omega }}_2^{\prime }$, and ${\mathrm{\Omega }}_2^{″}$, respectively) from the same single-mode tunable external cavity diode laser (ECDL2) are symmetrically arranged with respect to the $z$ axis [see Fig. 1(c)]. They intersect at the center of the atomic vapor cell (with a temperature of 120°C) to construct a hexagonal lattice with a period of $d\approx 111.8\mathrm{\mu m}$ in the $x\text{−}y$ plane. Here ${\mathrm{\Delta }}_i={\omega }_i-{\omega }_{ij}$ is the detuning between the resonant transition frequency ${\omega }_{ij}$ and the frequency ${\omega }_i$ of field ${\mathit{E}}_i(i=1,2)$, and Rabi frequency is defined as ${\mathrm{\Omega }}_i={\mu }_{\mathit{ij}}{E}_i/h$, where ${\mu }_{\mathit{ij}}$ ($i$, $j=1,2,3$) is the dipole momentum between transition $|i⟩\leftrightarrow |j⟩$ and $E_i$ is the amplitude of the electric field corresponding to ${\mathit{E}}_i$. A near-parallel probe beam ${\mathit{E}}_1$ with a Gaussian intensity profile from another ECDL1 is injected into the medium to cover the induced 2D lattice. With the laser frequencies tuned close to the two-photon resonance condition ${\mathrm{\Delta }}_1-{\mathrm{\Delta }}_2=0$, an electromagnetically induced transparency (EIT) window can be effectively generated to “write” a spatially periodic refractive index felt by the probe field. EIT refers to a phenomenon in which the resonant absorption of the atomic gas experienced by the probe field is significantly suppressed within a certain frequency width under the action of a strong coupling field [34]. When the probe detuning ${\mathrm{\Delta }}_1$ is set as $-100\mathrm{MHz}$, a discretized probe pattern can be observed at the output plane of cell around the EIT window. The output probe beam is imaged onto a charge-coupled device (CCD) camera [Fig. 1(a)] via an imaging lens.

The complex refractive index in an atomic EIT configuration is expressed as $n={(1+\chi )}^{1/2}\approx 1+\chi /2$, where $\chi ={\chi }^{\prime }+i{\chi }^{\prime \prime }$ is the susceptibility for describing the optical response of the probe field. The susceptibility of the system at probe frequency is [35] $$\chi =\frac{iN{|{\mu }_{31}|}^2}{\hslash {\epsilon }_0}\times {\left[({\mathrm{\Gamma }}_{31}+i{\mathrm{\Delta }}_1)+\frac{{|{\mathrm{\Omega }}_c|}^2}{{\mathrm{\Gamma }}_{32}+i({\mathrm{\Delta }}_1-{\mathrm{\Delta }}_2)}\right]}^{-1}.$$(1)

In the above expressions, ${\mathrm{\Gamma }}_{\mathit{ij}}$ is the decay rate between states $|i⟩$ and $|j⟩$ related to the longitudinal and reversible transverse relaxation times of each state, $|{\mathrm{\Omega }}_c{|}^2={|{\mathrm{\Omega }}_2+{\mathrm{\Omega }}_2^{\prime }+{\mathrm{\Omega }}_2^{\prime \prime }|}^2$ represents the Rabi frequency for the transition |$2⟩\leftrightarrow \text{|}3⟩$ driven by the hexagonal coupling field, and $N$ is the atomic density. The calculated $\chi$ is shown as in Figs. 1(d)–1(g), where both ${\chi }^{\prime }$ and ${\chi }^{\prime \prime }$ exhibit honeycomb distributions. As a result, the non-Hermitian photonic graphene is achieved, and the ratio ${\chi }^{\prime }/{\chi }^{\prime \prime }$ determining the degree of non-Hermiticity can be manipulated by controlling the laser detunings.

The interaction of the probe beam with the induced refractive index pattern is described within the paraxial approximation, mapping the Helmholtz equation to the Schrödinger equation for the field’s amplitude [36,37]. The susceptibility determines an effective potential in this equation, with real and imaginary parts proportional to each other, but with opposite signs. Strictly speaking, the system is not PT-symmetric, contrary to the previous works [12], but it still exhibits a similar non-Hermitian transition while being easier to implement. The transition is determined by $$\frac{\nabla \mathrm{\Gamma }}{\nabla U}=\frac{{\hslash }^2/2m{l}^2}{\mathrm{\Delta }\mathrm{\Gamma }},$$(2)where $U$ is the potential of the formed system and $l$ is a characteristic scale (lattice parameter). The ratio of the gradients of the imaginary and real parts of the potential is equal to the ratio of the confinement energy to the maximal difference in decay $\mathrm{\Delta }\mathrm{\Gamma }=l\nabla \mathrm{\Gamma }$ (between the minimum and the maximum). If $\nabla \mathrm{\Gamma }$ is increased beyond the value given by this expression, the system exhibits a non-Hermitian behavior (localization at the minima of the decay forming a hexagonal lattice). On the contrary, if $\nabla U$ is increased, the system recovers the Hermitian behavior (localization at the minima of the potential, forming a honeycomb lattice). This allows one to control the effective symmetry of the lattice in experiments via the degree of non-Hermiticity ${\chi }^{\prime }/{\chi }^{\prime \prime }=U/\mathrm{\Gamma }$, thanks to the direct correspondence between the susceptibility and the effective potential in the paraxial approximation.

To demonstrate this transition, we need to start by studying a quantum particle in an arbitrary complex potential (not periodic). The spatial dynamics of a wave packet in a complex potential can be described using the Schrödinger equation, $$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2m}\mathrm{\Delta }\psi +(U-i\mathrm{\Gamma })\psi ,$$(3)where $\psi (x,y)$ is the wave function corresponding to the amplitude of the electric field, $m$ is the particle mass, $U(x,y)$ is the real part of the potential, and $\mathrm{\Gamma }(x,y)$ is the imaginary part of the potential (we also define the decay rate in units of inverse time $\gamma =\mathrm{\Gamma }/\hslash$ for convenience).

At a small scale (for example, smaller than the period of the lattice), we can apply the Taylor expansion to the potential. The gradient of the real part of the potential acts as a force $F=-\nabla U$. This force can accelerate the wave packet in the direction opposite to the gradient, with the acceleration given by ${\mathit{a}}_H=\mathit{F}/m$ (the subscript $H$ stands for “Hermitian”). However, the dynamics of the wave packet is also affected by the gradient of the imaginary part $-\nabla \mathrm{\Gamma }$. The center of mass of the wave packet can be found as $${\mathit{r}}_c=\int \mathit{r}f(\mathit{r})\exp [\gamma (\mathit{r})t]\mathrm{d}\mathit{r},$$(4)where $f(\mathit{r})={|\psi (\mathit{r})|}^2$ is the probability distribution of the wave packet. For sufficiently small times, we can apply series expansion both to the exponent and to the decay rate, which gives $${\mathit{r}}_c\approx \int rf(\mathit{r})[1+\gamma ({\mathit{r}}_0)t+\nabla \gamma (\mathit{r}-{\mathit{r}}_0)·\mathit{r}t]\mathrm{d}\mathit{r},$$(5)where ${\mathit{r}}_0$ is the initial position of the wave packet center of mass. The first term in the brackets gives just ${\mathit{r}}_c={\mathit{r}}_0$. The second term also gives ${\mathit{r}}_c={\mathit{r}}_0$, but it changes the normalization of the wave function appearing because of the overall gain or decay of the wave packet and is of no interest for us. The last term is responsible for the displacement of the wave packet center of mass due to the change of the shape of the wave packet: the regions with higher decay are losing particles faster than the regions with smaller decay. For a wave packet with initial Gaussian distribution with a width $\sigma$, one finds the following displacement ${\mathit{r}}_c\approx {\mathit{r}}_0+\nabla \gamma {\sigma }^{2}t$. It means that the center-of-mass velocity (the “non-Hermitian velocity”) arising due to the gradient of the imaginary part of the potential is ${\mathit{v}}_{\mathrm{NH}}=\nabla \gamma {\sigma }^2$. In a system where both gradients (of the real and imaginary parts of the potential) are present, the velocity changes over time as ${\mathit{v}}_c={\mathit{v}}_{\mathrm{NH}}+{\mathit{a}}_{H}t$. If the two gradients have opposite directions, the wave packet first moves in the direction of the non-Hermitian velocity, then stops and starts to move in the direction of Hermitian acceleration. This occurs at any ratio of $\nabla \mathrm{\Gamma }$ and $\nabla U$, but the critical time $t_c$, given by $v_{\mathrm{NH}}=a_H{t}_c$ (when the Hermitian behavior starts to dominate), can be different.

Now we come back to a periodic lattice, with the minima of the real potential corresponding to the maxima of the decay. This system does exhibit a transition between the Hermitian and non-Hermitian regimes, appearing thanks to the presence of a characteristic length $l$ (the lattice parameter). The required time for the wave packet to collapse at the nearest decay minimum is $T_0=l/2v_{\mathrm{NH}}$. We therefore have the following equality for the critical behavior of the system, corresponding to the equilibrium between the Hermitian and non-Hermitian contributions, $$v_{\mathrm{NH}}=a_{H}l/2v_{\mathrm{NH}}.$$(6)

Replacing the variables in this expression by their definitions, we finally obtain the qualitative transition threshold given by Eq. (2).

In the case of a honeycomb lattice with a non-Hermitian part of the potential, this effect leads to the “inversion” of the lattice: from a honeycomb lattice defined by the real potential, the system switches to a simple hexagonal (“triangular”) lattice, with particles localized at the centers of the hexagons of the honeycomb lattice, where the decay is the smallest. To understand this behavior better, we have performed numerical simulations of a 2D honeycomb lattice first with a purely real, and then with a complex potential as described above. We solve the Schrödinger equation, Eq. (3), numerically, using a narrow Gaussian wave packet as an initial condition [38].

The dispersion for the Hermitian case is shown in Fig. 2(a), with the energy plotted in the unit of $t$, the effective tight-binding hopping parameter of the $\mathrm{s}$-band. The band structure is obtained by Fourier transforming the solution of Eq. (3), and the false color map corresponds to its intensity (squared absolute value of the wave function). The dispersion exhibits an $\mathrm{s}$-band and a $\mathrm{p}$-band (with an admixture of the higher bands, because the effective potential is not deep enough to provide a strong localization of the upper on-site levels). Figure 2(b) shows the simulated dispersion in the presence of a sufficiently strong non-Hermiticity, above the critical value given by Eq. (6): a strong signal from the $\mathrm{\Gamma }$ point of the $\mathrm{p}$-band means that this point has the longest lifetime (smallest decay rate). The real-space images of the eigenstates corresponding to the $\mathrm{\Gamma }$ points of the $\mathrm{s}$- and $\mathrm{p}$-bands are shown in Figs. 2(c) and 2(d), respectively. These images are obtained by plotting the spatial distribution of the solution at selected energies. They confirm that switching between these two points indeed corresponds to the inversion of the lattice (from honeycomb to hexagonal). Taking multiple bands into consideration is essential in our configuration because the lowest energy and the longest-living states belong to two different bands.

Figure 3 shows the transmitted probe field through the coupling lattice [Fig. 3(a)] at different ${\mathrm{\Delta }}_1$ around the EIT window, where both ${\chi }^{\prime }$ and ${\chi }^{\prime \prime }$ are honeycomb lattices. For ${\mathrm{\Delta }}_1=-60$ and $-90\mathrm{MHz}$, the observed output probe exhibits honeycomb intensity patterns [Figs. 3(b1) and 3(b2)]. According to the calculated index distribution, the ratio ${\chi }^{\prime }/{\chi }^{\prime \prime }$ gets over 10, indicating that the real part dominates; the wave packet is localized at the minima (the sites of honeycomb lattice) of the real potential $\mathit{U}\sim {\chi }^{\prime }$, which is verified by the numerical simulation given in Figs. 3(c1) and 3(c2). Also, the output intensity at ${\mathrm{\Delta }}_1=-90\mathrm{MHz}$ is stronger than that at ${\mathrm{\Delta }}_1=-60\mathrm{MHz}$, since the single-photon absorption at ${\mathrm{\Delta }}_1=-90\mathrm{MHz}$ (closer to the two-photon resonance ${\mathrm{\Delta }}_1-{\mathrm{\Delta }}_2=0$) can be better suppressed, which is one key function of EIT. With ${\mathrm{\Delta }}_1$ tuned very close to the EIT resonance, both ${\chi }^{\prime }$ and ${\chi }^{\prime \prime }$ are nearly zero [28], which indicates that the lattice potential can be very weak, and the output probe cannot demonstrate a clear discretization, which is experimentally verified by Fig. 3(b3).

When ${\mathrm{\Delta }}_1$ is further decreased to ${\mathrm{\Delta }}_1=-110\mathrm{MHz}$ and $-115\mathrm{MHz}$, ${\chi }^{\prime }$ and ${\chi }^{\prime \prime }$ also exhibit a honeycomb profile [as shown in Figs. 1(g) and 1(f)], but with a very different ratio, ${\chi }^{\prime }/{\chi }^{\prime \prime }$, compared to ${\mathrm{\Delta }}_1-{\mathrm{\Delta }}_2>0$. The ratio ${\chi }^{\prime }/{\chi }^{\prime \prime }$ at ${\mathrm{\Delta }}_1=-110\mathrm{MHz}$ is about 0.5 and the imaginary part dominates; the wave packet is localized at the minima of the decay $\mathrm{\Gamma }\sim {\chi }^{\prime \prime }$, inheriting from them an inverted-honeycomb (that is, hexagonal) output distribution. Both the observed [Fig. 3(b4)] and simulated [Fig. 3(c3)] hexagonal output patterns consistently support our analysis. Here Figs. 3(b4) and 3(b5) demonstrate the output hexagonal intensity patterns of the probe beam traveling through the non-Hermitian photonic lattice with a honeycomb susceptibility distribution. The output probe at ${\mathrm{\Delta }}_1=-115\mathrm{MHz}$ possesses a similar hexagonal pattern to ${\mathrm{\Delta }}_1=-110\mathrm{MHz}$, but with a weaker intensity due to being far from the two-photon resonance.

When the discretized probe images are efficiently obtained at the output plane of the cell ($z=0$), we monitor the propagation characteristics of the output field in free space at certain ratios of ${\chi }^{\prime }/{\chi }^{\prime \prime }$ by simultaneously moving the CCD camera and the imaging lens far away from the cell along $z$. Here, detunings ${\mathrm{\Delta }}_1=-90\mathrm{MHz}$ and $-110\mathrm{MHz}$ are two typical points where the Hermitian and non-Hermitian properties emerge. The captured probe images at different $z$ planes are shown in Fig. 4. The revival effect of honeycomb pattern obtained from a weak non-Hermitian case is shown in Fig. 4(a). For planes ($z={\mathit{Z}}_T/2\mathrm{mm}$) corresponding to the half-integer (fractional) Talbot effect, the inversion between hexagonal and honeycomb lattices occurs. The inversion between hexagonal and honeycomb patterns due to the fractional Talbot effect is similar to the half-Talbot effect in a 1D photonic lattice [27], where the output discrete pattern at $z={\mathit{Z}}_T/2$ can exhibit a half-period shift (along the transverse direction) compared to that at $z=0$, indicating the interchange of positions of bright and dark fringes. Here the observed inversion between hexagonal and honeycomb profiles also shows the interchange of the bright and dark sites on the pattern. Figure 4(b) demonstrates the Talbot effect of the output patterns under a strong non-Hermitian condition. Theoretically, the Talbot length for a hexagonal lattice [39] in a Hermitian case is ${\mathit{Z}}_T=m\times {3d}^2/(2{\lambda }_1)\approx m\times 32\mathrm{mm}$, where $\mathit{d}$ and ${\lambda }_1$ are defined as the period of the photonic lattice and wavelength of the incident field, respectively. Also, the honeycomb and inverted honeycomb lattices with the same spatial period can possess equal Talbot length [40]. The measured Talbot length for both cases is $\sim 30\mathrm{mm}$, which basically agrees with the calculation.

In order to quantitatively demonstrate the occurrence of the PT-symmetric-like transition, we plot the symmetry parameter obtained from the intensities at two different points associated with each of the two lattices: $I_1$ (honeycomb) and $I_2$ (hexagonal). The symmetry parameter written as $s=I_1/(I_1+I_2)$ gives $s=1$ for a purely honeycomb lattice and $s=0$ for a purely hexagonal one. Since the intensity scales as a square of the wave function, and the effective damping of a mode depends on its overlap with the imaginary part of the potential, which also scales as a square of the honeycomb wave function, one can expect the fourth root of the symmetry parameter to exhibit the typical square root scaling with the damping.

Figure 5 shows fourth root of the symmetry parameter $s^{1/4}$ as a function of relative non-Hermiticity, which is experimentally controlled by the detuning ${\mathrm{\Delta }}_1$. We present both the results of numerical simulations (black dots) and experimental measurements (red dots with error bars, extracted from the data of Fig. 3 and similar), together with the analytical curve corresponding to the expected square root scaling. The position of the transition point ($\approx 0.4$) depends on the properties of the particular potential, and on the effective mass of the particles involved; it is the only fitting parameter of the theoretical model. The good agreement between the two sets of data and the analytical scaling confirms our interpretation of the switching of the lattice symmetry as an example of the PT-symmetric-like transition.

In summary, we have demonstrated how instantaneously tunable and reconfigurable non-Hermitian photonic graphene is experimentally constructed relying on EIT in a coherent multilevel atomic medium. Such complex photonic graphene can be used as an ideal platform to demonstrate the underlying intriguing non-Hermitian, nonlinear, and quantum beam dynamical features (for example, topologic properties in a non-Hermitian system) in Dirac band structure geometry by exploiting the easy controllability of the dispersion properties and the Kerr nonlinearity in an EIT window [36,37,41]. We show a lattice symmetry switching from honeycomb to hexagonal due to the non-Hermitian transition. Contrary to previous works, this transition involves two different energy bands.

The non-Hermitian transition in this 2D artificial complex periodic structure is further evidenced by the Talbot effect. The observed complex Talbot effect can be promisingly used to design an all-optical controllable Talbot–Lau interferometer. What is more, this configuration can be easily extended to ultracold atomic samples and potentially serve as an imaging technique to exhibit the distribution of ultracold atomic gas with a better resolution due to the more intensive spots (compared with the 1D Talbot effect).

Acknowledgment. G. Malpuech and D. D. Solnyshkov acknowledge the support of the projects EU “TOPOLIGHT” and “QUANTOPOL,” ANR “Quantum Fluids of Light,” ANR Labex GaNEXT, and of the ANR program “Investissements d’Avenir” through the IDEX-ISITE initiative 16-IDEX-0001.