Momentum space AC Josephson effect and intervalley coherence in multilayer graphene

Domain: physics
Published: 2025-12-24
Editorial outcome: minor_revision

Abstract
Electrontransportdrivenbythephasecoherenceandinterferenceofquantummany-bodywavefunctionsisafascinatingphenomenonwithpotentialtechnologicalsignificance.Superconductivity,forexample,enablesdissipationlesstransportthroughmacroscopicphasetwisting,whilecharge-densitywavesshowbroadbandnoiseandAC-DCinterferenceoncetheirphaseisdepinned.Weidentifyananalogousphase-coherentdynamicsintheintervalley-coherent(IVC)statefrequentlyobservedinrhombohedralmultilayergraphene.Weshowthatastaticmagneticfielddrivesanoscillatingintervalleycurrent,whichproducesanoscillatingorbitalmagnetizationandtherebyadetectableACHallresponse.ThismechanismissimilartotheJosephsoneffectobservedinsuperconductorsbutnowhappeninginmomentumspace.Inthisanalogy,thestaticmagneticfieldactsastheDCvoltage,whiletheoscillatingintervalleycurrentistheACJosephsoncurrent.Weprovidemicroscopiccalculationsofallthephase-numberfree-energyparametersinrhombohedraltrilayergraphene,establishingthiseffectasanexperimentalsignatureofIVCorder.

Introduction
Magnetism arising from spin and valley ordering in multilayer graphene systems, both with and without moiré potential, has become a well-established phenomenon following the discovery of the integer quantum anomalous Hall effect in twisted bilayer graphene. There are three primary types of spin-valley ordered states: spin-polarized ferromagnets, valley-polarized orbital ferromagnets, and the intervalley-coherent (IVC) state. Spin-polarized ferromagnets typically expand their phase boundaries under in-plane magnetic fields, whereas valley-polarized orbital ferromagnets are known for exhibiting a strong anomalous Hall effect. These states, characterized by Ising-like order parameters at temperature below the magnetic anisotropic energy, also generate magnetic stray fields detectable via nano-SQUID. In contrast, the intervalley-coherent (IVC) state, characterized by its XY-like order parameter, remains elusive and is typically identified by eliminating other symmetry-broken states and through comparison with mean-field theory. Theories suggest that IVC fluctuations could be the microscopic origins of superconductivity in multilayer graphene and homobilayer moiré transition metal dichalcogenides. The IVC state was previously thought to exhibit no prominent transport signatures and can mainly be detected through STM experiments. We note that the difference in the dielectric screening environment can influence the observable magnetic ground states in STM versus transport experiments. In this study, we report a unique transport phenomenon associated with the IVC state: an out-of-plane static magnetic field induces an oscillating orbital magnetization, leading to an AC Hall effect, as shown in Fig. 1. This effect, analogous to AC Josephson effect in a superconductor, can be detected using probes sensitive to finite-frequency Hall conductivity \sigma_{xy}(\omega), such as microwave impedance microscopy, and may serve as a definitive diagnostic tool for identifying the IVC state.

Results
Momentum-space AC Josephson effect
In multilayer graphene with a low electron density of about 10^{-4} electrons per carbon atom (typical in most experiments), electrons occupy a tiny fraction of the Brillouin zone, specifically around the opposite corners, which are denoted as valleys \tau: K, K' = \pm K. To make the connection to the BCS wave function more transparent, we now perform a particle-hole transformation in one of the valleys and denote its transformed vacuum by |\Psi_{vac}|. Then, the IVC ground state can be written analogously to a BCS superconductor, except that the pairing occurs between electrons and holes in opposite valleys, rather than between electrons with opposite spins and momenta:
|\Psi_{\phi}| = \prod_{k} (u_{k} + v_{k} e^{i\phi} c^{\dagger}_{K+k} c_{K'+k}) |\Psi_{vac}|.
Here, k runs over a small subset of wave-vectors (reflecting the low electron density) and u_{k} = \cos\theta/2 and v_{k} = \sin\theta/2 are real coefficients representing valley-mixing angle. The particle-hole conjugated vacuum is
|\Psi_{vac}| = \prod_{k} c^{\dagger}_{K'+k} |0|.
The operator c^{\dagger}_{K+k} c_{K'+k} acting on |\Psi_{vac}| creates an exciton with wave-vector 2K. Thus, |\Psi_{\phi}| describes a condensate of such 2K excitons with \phi representing the global phase of the condensate. The spin degree of freedom is not important to our discussion, and it suffices to treat the electrons and holes as having the same spin projection. Analogous to the BCS state in superconductors, the state |\Psi_{\phi}| possesses a well-defined global phase but an indefinite exciton number. The wave function with a definite number of excitons can be obtained by integrating |\Psi_{\phi}| over the global phase after multiplying by an appropriate phase factor, effectively “localizing" the state in particle-number space:
|\Psi_{N\tau}| = \int \frac{d\phi}{2\pi} e^{iN\tau\phi} |\Psi_{\phi}|.
Here |\Psi| is not normalized. This exciton picture of electron-hole pairing at distinct wave-vectors has been widely used to describe one-dimensional charge-density-wave (CDW) metals like NbSe_{3} (e.g., ref. 55). However, there is an important difference between 1D CDW and our 2K exciton condensate in higher dimensions. In 1D, the number eigenstate with N excitons carries a net velocity of 2N\hbar k/m* (m* is the effective mass) and thus describes a current-carrying state. In our case, however, if the electron state |\Psi_{N\tau}| does not carry any net velocity, it generally does not produce a measurable net current. However, two-dimensional materials with low-energy states at the honeycomb Brillouin zone corners, such as rhombohedral multilayer graphene, are special. In this system, Bloch states at the \pm K points with the same energy acquire opposite orbital angular momenta when the Dirac points are gapped by a displacement field, interactions, or both. In this case, an exciton formed from an electron and a hole in opposite valleys has the same chirality—unlike the textbook situation where electrons and holes bend in opposite directions, here they bend in the same direction. As a result, the exciton number eigenstate |\Psi_{N\tau}| carries a finite orbital moment with measurable consequences. In particular, an external magnetic field B shifts the energies of opposite valleys in opposite directions, effectively acting as a valley-contrasting chemical potential:
F_{B} = -M B = -\frac{1}{2} g_{v} \mu_{B} N\tau B,
where M = \frac{1}{2} g_{v} \mu_{B} N\tau is the net orbital moment. We substantiate this relationship through a detailed microscopic calculation in the following sections.
Another unusual aspect of a 2K exciton condensate is the presence of intrinsic intervalley tunneling from Umklapp scattering. While 4K is not a reciprocal lattice vectors, 6K is. Consequently, lattice translation symmetry prohibits two-body Umklapp scattering but allows three-body Umklapp scattering. This process is analogous to the Josephson coupling in Josephson junctions, which maintains phase coherence across a junction by tunneling Cooper pairs. Similarly, three-body Umklapp scattering enforces phase coherence in the IVC exciton condensate by tunneling three electrons from one valley to another. This interaction introduces an intrinsic potential energy that maintains phase coherence:
F_{ump} = -\gamma \cos(3\phi).
Equations (4) and (5), along with the superfluid kinetic terms,
F_{0} = \frac{\rho}{2} (\nabla\phi)^{2} + \frac{1}{2\chi} N^{2}_{\tau},
where \rho and \chi are the superfluid stiffness and valley susceptibility, define the unusual transport phenomena in IVC states. They define a phase-number action, F = F_{IVC} + F_{0} + F_{B}, that applies to any 2K exciton-condensate with electrons and holes have the same chirality:
F_{IVC} = \frac{\rho}{2} (\nabla\phi)^{2} + \frac{1}{2\chi} N^{2}_{\tau} - \gamma \cos(3\phi) - \frac{1}{2} g_{v} \mu_{B} N\tau B.
The phase coherent transport of superfluid can be obtained by treating N\tau and \phi as dynamical conjugate variables:
\hbar \frac{d\phi}{dt} = \frac{\partial F_{IVC}}{\partial N\tau} = \frac{N\tau}{\chi} - \frac{1}{2} g_{v} \mu_{B} B,
\hbar \frac{dN\tau}{dt} = -\frac{\partial F_{IVC}}{\partial \phi} = \rho \nabla^{2}\phi - 3\gamma \sin(3\phi).
Equations (8) and (9) represent the main results of our article, describing phase-coherent dynamics in IVC. We focus on a uniform solution, then this differential equation has stable fixed points at (N*\tau, \phi*) = (\hbar \omega_{B} \chi, 2n\pi/3) and saddle-points at (N*\tau, \phi*) = (\hbar \omega_{B} \chi, (2n+1)\pi/3) with n \in Z, as indicated in Fig. 2a. Here \omega_{B} = g_{v} \mu_{B} B/2\hbar is the orbital Larmor frequency. Evaluating the free energy at the saddle point gives the critical energy F* = \gamma + \frac{1}{2} \chi (\hbar \omega_{B})^{2}. When the free energy at the initial condition \phi(0) = N\tau(0) = 0 exceeds F*, i.e., B > \frac{g_{v} \mu_{B} \chi}{\gamma}, the phase \phi keeps precessing, producing a Josephson-like effect shown in Fig. 2b, c. The static magnetic field acts like a chemical potential for the exciton, inducing a nearly linear, coherent phase rotation of the condensate over time. This phase dynamics then subsequently generates an oscillating exciton number N\tau(t), through intervalley tunneling, leading to an oscillating orbital magnetization, M(t) = \frac{1}{2} g_{v} \mu_{B} N\tau(t). If the free energy at the initial condition \phi(0) = N\tau(0) = 0 is smaller than F* (i.e., B < \frac{g_{v} \mu_{B} \chi}{\gamma}), the system instead remains trapped near a stable fixed point \phi* = 2n\pi/3, exhibiting small oscillations around it. In both regimes, the dynamics inevitably lead to an oscillating orbital magnetization.
Table 1 summarizes the parallels between the AC Josephson effect in a superconductor and the analogous phenomenon in our 2K exciton condensate. For small N\tau(t)/\chi and B > \frac{g_{v} \mu_{B} \chi}{\gamma}, the solution is straightforward: \phi(t) = \omega_{B} t and M(t) = M_{0} [1 - \cos(3\omega_{B} t)] where M_{0} = \frac{1}{2} g_{v} (\gamma/\hbar \omega_{B}) \mu_{B}. Note the amplitude (frequency) of M(t) decreases (increases) linearly with external field B. The minor nonlinear feature in Fig. 2b results from the finite N\tau(t)/\chi.

Rhombohedral trilayer graphene (RTG)
We now demonstrate how to compute all the parameters in F_{IVC} for IVC states in RTG. We start with the many-body Hamiltonian H = H_{b} + V, where
H_{b} = \sum_{k, \sigma, \sigma', \tau} (\tau k_{x}, k_{y}) c^{\dagger}_{\tau, \sigma, s} (k) c_{\tau, \sigma', s} (k).
The interaction is modeled as a density-density type:
V = \frac{1}{2} \sum_{q} V_{q} n_{q} n_{-q},
where n_{q} = \sum_{k, \sigma, \tau, s} c^{\dagger}_{\tau, \sigma, s} (k+q) c_{\tau, \sigma, s} (k) and V_{q} = 2\pi k_{e} / (\epsilon_{r} |q|) is the Coulomb potential. The anisotropic interaction is suppressed in the strong electric-displacement limit, and it is not important for discussing IVC transport. Here, k_{e} = 1.44 eV nm is the Coulomb constant, \epsilon_{r} = 15 is the screening constant, and d = 50 nm is the distance from the gate to the material. Next, we solve this Hamiltonian using the Hartree-Fock approximation, H_{MF}, where the matrix elements of H_{MF} are determined self-consistently. For a detailed discussion of the self-consistent Hartree-Fock calculations, including the discussion of Fast Fourier Transformation trick to obtain the self-energy, see ref. 63.
The ground state thus obtained at (n_{e}, U) = -2.5 \times 10^{11} cm^{-2}, 15 meV is an intervalley coherent (IVC) state, consistent with previous mean-field calculations. In these calculations, we did not perform a particle-hole conjugation in one of the valleys, following the convention in prior works to allow for direct comparison. A line cut of band structure for the IVC running wave is shown in Fig. 3.
Next, we introduce a Lagrange multiplier into the Hamiltonian H:
H_{\lambda} = -\lambda (|K><K| - |K'><K'|).
The resulting band structure with \lambda = 0.02 meV is shown in Fig. 4a. The introduction of \lambda results in a converged solution with a net valley imbalance N\tau = \rho_{\lambda} = \int \frac{d^{2}k}{(2\pi)^{2}} \sum_{n} f(\epsilon_{\lambda, n, k}) |\psi_{\lambda, n, k}|^{2}.
(Here \epsilon_{\lambda, n, k} and \psi_{\lambda, n, k} are the mean-field quasiparticle energy and eigenstate). The resulting band structure with \lambda = 0.02 meV is shown in Fig. 4a. The introduction of \lambda results in a converged solution with a net valley imbalance N\tau = \rho_{\lambda}. The energy increase is E = \frac{1}{2} Tr[(H_{b} + H_{MF}) \rho_{\lambda}], as shown in Fig. 4b, c. Importantly, this state we simulated with finite \lambda represents a many-body state where the energy increase due to depletion of excitons from the condensate, N^{2}_{\tau} = (2\chi), is counterbalanced by the energy lowering of the orbital magnetization -MB. (In valley-pseudospin language, this corresponds to an excited state with a slight valley imbalance, where the valley pseudospin order parameter deviates a little bit away from the XY-plane). This scenario represents a stationary solution of Eq. (8) where \phi = 0 and N* = \frac{1}{2} \chi g_{v} \mu_{B} B. Here M is computed using the so-called modern theory for orbital magnetization:
M = \frac{e}{\hbar} \int d^{2}k \sum_{n} f(\epsilon_{\lambda, n, k}) \text{Im} <\psi_{\lambda, n, k} | \frac{\partial H_{MF}}{\partial k_{x}} | \psi_{\lambda, n', k} > <\psi_{\lambda, n', k} | \frac{\partial H_{MF}}{\partial k_{y}} | \psi_{\lambda, n, k} > / (\epsilon_{\lambda, n, k} - \epsilon_{\lambda, n', k})^{2}.
The linear and quadratic fitting shown in Fig. 4b, c yield the parameters for the phase-number action in RTG as g_{v} \sim 3.97, \chi^{-1} \sim 0.126 meV. They are used to generate Fig. 2.

3-body Umklapp scattering
At dilute doping where charge carriers are concentrated in the two opposite valleys, the Umklapp process that transfers electrons between these valleys necessitates at least a three-body scattering process. This is because a two-body scattering process does not provide sufficient momentum to match a reciprocal lattice vector. Such large-momentum transfer processes means a delta-function-like potential in the K-p continuum approximation. At second order in Coulomb potential and electron-phonon coupling, only one diagram contributes to this three-body scattering, as shown in Fig. 5. Under a lattice translation a, c \rightarrow c e^{iK \cdot a}, this diagram picks up a phase e^{i6K \cdot a} = 1 since 6K is a reciprocal lattice vector. This scattering process begins with three electrons all initially in the same valley (K). These electrons undergo two distinct scattering events to transition to the opposite valley K' (modulo G). First, two of the electrons scatter, moving one to the final valley K' (modulo G) and the other to an intermediate state at \Gamma point with a matrix element V(K). Subsequently, the electron at the \Gamma point scatters with the remaining electron at K, transferring both to valley K' (modulo G) with matrix elements V(K). The strength of this scattering amplitude can therefore be estimated as follows:
\gamma \approx \frac{|V(K)|^{2}}{|\epsilon_{K} - \epsilon_{\Gamma}|} \approx 0.2 \mu eV/charge.
In the last line, we used V(K) \approx 2\pi k_{e} / |K| = 1.33 meV at |n_{e}| = 2.5 \times 10^{11} cm^{-2} and |\epsilon_{K} - \epsilon_{\Gamma}| = 9 eV.
This type of lattice-scale interaction receives higher-order contributions from both electron-electron and electron-optical phonon interactions. We and others have estimated the strength of two-body lattice-scale interactions, finding them significantly enhanced by the large density of states at graphene’s \pi band. In fact, without at least an order of magnitude enhancement in two-body lattice-scale interactions, their magnitudes would be insufficient compared to the Zeeman energy to explain the experimental findings. Figure 8 in ref. 51 attempts to quantify this renormalization. Therefore, by treating \gamma as a phenomenological constant and leveraging our better theoretical control of other constants in F_{IVC}, we propose detecting the oscillating magnetization induced by a static magnetic field as a practical experimental method to determine this parameter.

Damping
The orbital magnetization oscillation cannot persist indefinitely; it eventually dissipates energy and settles into a partially valley-imbalanced state. Below, we discuss three distinct channels through which the order-parameter dynamics can lose energy: finite-temperature effects, atomic-scale disorder, and intrinsic damping.
(1) Finite temperature effect: the universality class of F_{IVC} is XY-like, and its topological defects are phase singularities known as merons. Below the Kosterlitz-Thouless (KT) transition temperature T_{KT}, or when \gamma is sufficiently strong, merons and antimerons are bound, preventing phase slippage. However, for T > T_{KT}, unbound merons cause valley imbalance through phase slippage. This binding is crucial for the observability of our proposed phase-coherent transport phenomena. We have estimated T_{KT} based on a microscopic Hartree-Fock calculation of IVC domain wall. This is orders of magnitude higher than the lowest temperatures accessible in many 2D transport experiments, making our proposed phenomena experimentally feasible.
(2) Atomic-scale disorder: atomistic defect in graphene such as missing carbon atoms or sp^{2} \rightarrow sp^{3} hybridization distort the nearby electronic charge and bond densities. Depending on the defect symmetry, the resulting pattern can range from a uniform Kekulé bond order to a Kekulé vortex with a winding phase profile. The latter is also directly visualized by STM in monolayer graphene’s zeroth Landau level. For concreteness, let us consider a uniform Kekulé pattern generated by a divacancy. Ab initio calculations estimate the structural energy cost at 7–8 eV per atom. Unlike conventional charge-density-wave systems, which occur near commensurate fillings (e.g., 1/2, 1/3, 1/4), gated graphene has carrier densities that are orders of magnitude smaller. For n_{e} = 10^{12} cm^{-2}, with unit-cell area A_{uc} = \sqrt{3}a^{2}/2 \approx 5.24 \mathring{A}^{2}, the filling fraction is only \nu = 1.3 \times 10^{-4} \pi-band electrons per carbon atom. Thus, the probability of a carrier occupying a defect site is vanishingly small. The corresponding defect-formation energy per electron is V_{im,j} \approx 8 \times 1.3 \times 10^{-4} eV \approx 1 meV. Within a continuum description, disorder pinning contributes to the free energy as localized potentials acting on the phase field at defect positions r_{j}:
F_{imp} = \sum_{j=1}^{N} V_{im,j} \cos(\phi(r_{j}) - \theta_{j}).
where \theta_{j} is the locally preferred IVC phase angle and V_{im,j} \approx 1 meV the effective pinning strength. While the atomic-scale defect energy (several eV) will still be strongly reduced by the large device area, atomic-scale pinning can be neglected. Recent STM studies have revealed bilayer graphene devices with defect-free regions exceeding 30 \times 30 nm^{2}, confirming the experimental feasibility of accessing the clean limit. On the transport side, transverse magnetic electron focusing experiments and the observation of family of unconventional superconductivity, provide further evidence that multilayer graphene can indeed reach the very clean limit with mean-free path similar to the device dimension.
(3) Intrinsic damping: when extrinsic defects and finite-temperature effects are negligible, oscillations are ultimately limited by intrinsic damping. A useful perspective is to view the dynamics as an analog of ferromagnetic resonance in magnetism. On the valley Bloch sphere, the poles correspond to time-reversal-broken states with majority occupation in valley K or K', while the equator (the XY-plane) represents IVC states. If the ground state is polarized at the north or south pole, there is no natural external field that couples directly to the IVC order parameter since valley is fundamentally just crystal momentum. Coherent rotations about the poles—valley-wave excitations—are very difficult to excite, since an external probe would need to supply the large momentum transfer between K and K' (c^{\dagger}_{K} c_{K'}). By contrast, if the ground state is an XY-valley ferromagnet (an IVC state), a uniform out-of-plane magnetic field B acts as a valley-contrasting chemical potential that couples directly to valley-imbalance. In this analogy, the term N^{2}_{\tau}/2\chi provides the uniaxial valley anisotropy energy favoring N\tau = 0, while the three-body Umklapp term F_{ump} provides the in-plane valley anisotropy energy to pin the IVC phase when B = 0. At finite B, the dynamics of the orbital magnetization closely resemble ferromagnetic resonance: an out-of-plane B field drives oscillations of the orbital magnetization and produces strong absorption of an ac field B_{ac} at three times the orbital Larmor frequency. The intrinsic damping can then be interpreted as a valley analog of the Kamberský mechanism, where time-dependent valley-order parameter dynamics generates particle-hole excitations near the Fermi surface, and quasiparticle decay channels transfer energy from the collective mode to the lattice. The valley analog of Gilbert parameter modifies the order-parameter dynamics as follows:
\dot{\phi} = \frac{N\tau}{\hbar\chi} - \frac{1}{2\hbar} g_{v} \mu_{B} B + \alpha_{\phi} \dot{N\tau},
\dot{N\tau} = -\frac{3\gamma}{\hbar} \sin(3\phi) + \alpha_{N\tau} \dot{\phi}.
Estimating the damping parameters is more subtle than for the other free-energy terms, since they depend on the quasiparticle lifetime \hbar\tau^{-1}. For a very conservative estimate of \hbar\tau^{-1} \approx 0.1E, we estimate the damping parameter to be \alpha_{\phi} \approx 4 \times 10^{-3} and \alpha_{N\tau} \approx 10^{-8}. The difference in scale reflects the microscopic origin of the two terms: \alpha_{N\tau} is suppressed because it requires processes that violate valley number conservation and is therefore at least second order in the three-body Umklapp term, while \alpha_{\phi} is contributed bandstructure details such as trigonal warping. Since \alpha_{N\tau} is practically zero the damping dynamics is mainly controlled by \alpha_{\phi}. Figure 6 shows numerical solutions of Eqs. (17) and (18) for an external field B = 0.15 T with damping. The system undergoes thousands of ac cycles before dissipating enough energy to reach the critical energy surface (red dashed line in Fig. 2a). Beyond this threshold, the trajectory spirals slowly into a stable fixed point and the oscillations cease. We note that the orbital magnetization at the stable fixed point agrees with the steady-state solution of the constrained Hartree-Fock mean-field calculation with a Lagrange multiplier \lambda = \frac{1}{2} g_{v} \mu_{B} B, where B is the initial condition used in solving the dynamical equations.

Methods
We use constrained Hartree-Fock methods with Lagrange multipliers to compute the valley stiffness \rho and the anisotropic susceptibility \chi^{-1} associated with collective valley rotations. The orbital magnetization is obtained from the mean-field wave functions and eigenenergies of the constrained Hartree-Fock calculation, following the justification in refs. 101–103. The strength of the three-body Umklapp scattering term is estimated using second-order perturbation theory. The damping parameter is evaluated using Kambersky’s theory in the intraband limit. The nonlinear equations of motion are solved using standard adaptive Runge-Kutta methods.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have answered my questions, and that of the other Referees, very satisfactorily. The corrections, clarifications\nand new material added to the main text make the work stronger in my opinion. In particular, I have found the discussion of\nthe dynamics in the presence of damping, and how to recover mean-field results in the long-time limit, to be very illuminating\nand of clear importance to future experimental efforts aimed at measuring the proposed orbital magnetization oscillations.\nI enthusiastically support publication in Nature Communications in its current form, with very minor comments appended\nbelow.\n(1) When referring to the damping numerics on line 383, I assume the authors mean Fig 5 rather than Fig 1?\n(2) The authors have added a comment on the importance of C_2 T symmetry in response to one of the other referees, but I\ndon't think that's correct. I got initially quite confused about this statement because obviously a perpendicular displacement\nfield D cannot break an in-plane rotation symmetry.\nThe resolution is that C_2 rotations are not a symmetry of rhombohedral graphene due to its uniform ABC stacking order (as\nopposed to TBG for example). True (three-dimensional) inversion symmetry however is a symmetry of the rhombohedral\nstructure. Applying a displacement field breaks the 3D inversion (rather than breaking C_2), which permits a Berry curvature\ndistribution to arise in each valley (as long as it is compensated between the two valleys to preserve time reversal).\nI leave it to the author's discretion how much (if any) of this discussion is actually useful, but they should at least correct the\nstatement as written.\n(3) I am not sure how the authors are doing on length limits with the new material added but, if necessary, I think shortening\nthe Damping section (with details relegated to the Supplement) could be advisable: at the moment it is quite long compared\nto other sections of the main text.", "reviewer_3": "The authors have satisfactorily addressed all my comments. Given that this paper contains an interesting and original idea\nabout a timely topic, I recommend publication."}, "author_response": ""}]

Decision letter