8 not causally connected points) that is independent of the distance, and a spatially dependent (causal) part, that expands in the form of a wave front. In Fig. 7 b) we use g ′ = 0 . 5, where at low energies the spectrum of excita- tions is gapped. Although the dispersion in this model is linear (relativistic) only at sufficiently high energies above the gap we still see a clear light-cone structure. The expansion velocity of the “light cone” in this case is consistent with twice the maximum of the group velocity v$_{gr}$ ( k ) = dϵ/dk given by v$_{gr,max}$ = { 2 Jg for | g | < 1 2 J for g | ≥ 1. (30) So for g ′ = 1 we find 2 v$_{gr,max}$ = 4 J , and for g ′ = 0 . 5 we find 2 v$_{gr,max}$ = 2 J . These are indeed the expansion velocities that we see in Fig. 7. FIG. 8: Long-time behavior of the correlation function for T = 1, κ = 8, and V = 80. The correlation function first develops algebraic scaling, so the system forms a metastable quasi-superfluid state. On longer time scales the correlation function shows exponential decay. The coherence is lost due to dynamical vortex unbinding. V. DYNAMICAL VORTEX UNBINDING In this section we address the important question of how the supercritical state relaxes to the ground state, i.e. the second stage of the dynamics. As we mentioned in the introduction the anticipated mechanism for this re- laxation is vortex unbinding. This process is intrinsically nonlinear and requires a more sophisticated treatment than that of the noninteracting “light cone” dynamics. In this work we use two complimentary approaches. In Sec. V A we use a numerical implementation of the TWA to simulate the dynamics in the system. In Sec. V B we generalize a renormalization group approach to analyti- cally describe the process of relaxation in real time. description, we conclude that the dynamical transition that we observe is due to vortex unbinding. The exam- ple of a single run shown in Fig. 2 is consistent with this picture: Defects are created soon after the quench, but they only gradually separate on a much longer times scale. It is this process that we refer to as the reverse Kibble-Zurek mechanism. To better characterize the process of vortex unbinding further we fit the correlation function G ( x, t ) to either algebraic or exponential fitting functions. Such choice is motivated by the two possible regimes of the equilibrium system and is supported by the analytic renormalization group results presented in the next section. The alge- braic fitting function we use is c ( L/π | sin( πx/L ) | ) − τ/ 4 and the exponential function is c exp( −| sin( πx/L ) /x$_{0}$ ). Note that in the fitting functions we use the conformal distance L/π | sin( πx/L ) | ), which is more appropriate in finite systems with periodic boundary conditions (see e.g. Ref. [33]). In equilibrium the algebraic exponent τ would be the relative temperature T/T$_{c}$ . Any value above 1 is therefore supercritical. The parameter x$_{0}$ defines the length scale of the exponential decay. The parameter c in both functions gives an overall scale. Using these fitting functions we analyze four different situations corresponding to the same initial temperature T = 1 and the same parameter κ = 8, but with dif- ferent initial couplings V between the planes. The first (I) case corresponding to V = 80 is identical to the one plotted in Fig. 8. The other three curves correspond to V = 70 , 50 , 20 (II–IV). In Fig. 9 we show the expo- nent τ extracted from the fit as a function of time for these situations. In all of them at short times G ( x, t ) develops algebraic scaling when the light-cone dynamics reaches the system boundaries. For the cases I–III the emerging scaling exponent τ is well above the critical ex- A. Numerical approach Within TWA we need to solve the full nonlinear equa- tions of motion (3) subject to the initial conditions dis- tributed according to the Wigner function (4). Then the equal-time correlation functions or other observables are found by averaging the Weyl symbol of the corresponding observable computed at time t over the fluctuating ini- tial conditions. Note that since we are interested only in phase-phase correlation function the corresponding Weyl symbol is obtained by simply substituting the Heisenberg quantum operator corresponding to the phase with the classical phase $^{31}$. In Fig. 8 we show the result of such simulations. We can observe how the metastable super- fluid state relaxes to the disordered state. For that, we show the correlation functions of the system on a much longer time scale than in Fig. 3. The exponent of the al- gebraic scaling gradually decreases. Eventually the cor- relation function is more accurately approximated by an exponential fitting function, signalling that the thermal Bose gas phase has been reached. Because this is the phase of deconfined vortices, and because the intermedi- ate superfluid phase is well described by a phonon-only