\section{Introduction}

The research on Dirac materials~\cite{vafek2013, wehling2014} encompasses many areas of physics, ranging from a fundamental theory with ties to particle physics, to concrete sought-after technological applications in materials science. From a theoretical point of view, interacting
Dirac systems are particularly appealing as a playground for development of novel and applications of known  quantum-field-theoretical tools, providing an opportunity for their testing against numerical and experimental results.

Graphene\cite{Novoselov2005,Geim2007,castroneto2009}, in particular, at charge-neutrality serves as a prime example of a class of materials where Dirac nature of the low-energy excitations emerges naturally from the underlying honeycomb lattice structure, and is further protected by its symmetries.
The presence of sufficiently strong interactions can lead to spontaneous symmetry breaking signaled by the appearance of an ordered ground state.
There is a wide variety of possible ordering patterns determined by small variations of the short-ranged interaction parameters, which have been analyzed by different theoretical approaches, ranging from mean-field calculations to the renormalization group and numerical simulations\cite{sorella1992,herbut2006,honerkamp2008,herbut2009, raghu2008, juricic2009, assaad2013, chandra2013, parisen2014, janssen2015, sorella2015, wang2016, hasselmann2016, volpez2016,sanchez2016}.
A large repulsive onsite interaction $U$, for example, exhibits a continuous quantum phase transition from the semimetallic (SM) phase into an antiferromagnetic spin-density-wave (SDW) state, whereas a nearest-neighbor interaction $V_1$ would induce a conventional charge-density wave (CDW) state.
While there is an ongoing debate on whether a pure second-nearest-neighbor interaction $V_2$ can induce a topological Mott insulating state\cite{raghu2008,duric2014} or whether it favors a sublattice charge-modulated state\cite{grushin2013,daghofer2014, capponi2015, motruk2015},
it is believed that the interplay between balanced $V_1$ and $V_2$ gives rise to yet another ordered ground state -- the Kekul\'e valence bond solid (VBS)~\cite{hou2007, roy2010c, weeks2010, garciamartinez2013}. Alternatively, this ordering can also emerge from electron-phonon interactions\cite{nomura2009,kharitonov2012,Classen2014b}.
The translation symmetry breaking Kekul\'e VBS corresponds to a particular dimerization pattern of the fermions on the honeycomb lattice with a bond-dependent nearest-neighbor hopping amplitude (see Fig.~\ref{fig:kekule}), and has recently been accessed even experimentally~\cite{Gutierrez2016}.

\begin{figure}[t!]
\includegraphics[height=0.45\columnwidth]{fig01a_kekule.pdf}
\hspace{0.2cm}
\includegraphics[height=0.45\columnwidth]{fig01b_freeenergy.pdf}
\caption{Left panel: An example of the Kekul\'e dimerization pattern. Shorter bonds are marked by thicker red lines and longer bonds by the thinner gray lines. Right panel: LG free energy across the phase transition, where the dotted (dashed) line corresponds to the disordered (ordered) phase, with $r>r_c\ (r<r_c)$.}
\label{fig:kekule}
\end{figure}

The Kekul\'e ordering pattern can be described by a $Z_3$ symmetric order parameter that describes breaking of the lattice translational symmetry and the $C_3$ rotational symmetry\cite{li2015}. Due to $Z_3$ symmetry the cubic terms in the order parameter in the Landau-Ginzburg (LG) free energy are allowed, which generically prevent a phase transition toward the ordered state from being continuous\cite{Golner1973}.
The LG free energy $F_{\text{LG}}$ takes the form
\begin{align}\label{eq:LGfe}
	F_{\text{LG}}[\phi]=F_0+r\, |\phi|^2 + g (\phi^3+\phi^{\ast 3})+\lambda |\phi|^4\,,
\end{align}
where $\phi$ represents the complex-valued order parameter of the VBS, $F_{\text{LG}}$ is symmetric under the field transformation $\phi \to \exp(2\pi i\, n/3)\phi, n\in \mathbbm Z$ and the parameter $r$ is used to tune through the phase transition by crossing a critical value~$r=r_c$.
In the case when the cubic coupling $g$ vanishes, the LG theory suggests that the phase transition is of second order, since at $\phi=\phi_0\neq0$ a non-vanishing minimum of the LG free energy emerges continuously from the origin at $\phi=0$.
For a finite~$g$, on the other hand, the minimum of $F_{\text{LG}}$ jumps discontinuously from $\phi=0$ to a $\phi_0\neq0$ as $r$ goes through its critical value~$r_c >0$, i.e. a first-order transition is expected (see Fig.~\ref{fig:kekule}). This simple consideration is not sufficient, however, to describe the Kekul\'e phase transition when massless fermions are present and coupled to the $Z_3$ order parameter.

We study the quantum phase transition of massless Dirac fermions toward the Kekul\'e VBS within the appropriate Gross-Neveu-Yukawa theory, and in $D=3+1-\epsilon$ dimensions. We find that the presence of (neutral) Dirac fermions tends to reduce the first-order character of the phase transition\cite{li2015}, and that closer to $D=2+1$ dimensions the phase transition may even be rendered continuous\cite{Roy2013c}. However, the values of parameter $\epsilon$ where this happens are not small, and this possibility seems beyond the realm of control of the $\epsilon$ expansion. Sizable corrections which may even overturn the conclusion based on the lowest order calculation
can generally be expected beyond the one-loop order. Adding a $U(1)$ gauge field that would couple exclusively to fermions, as the Kekul\'e order parameter is electrically neutral, on the other hand, significantly further reduces the discontinuous nature of the transition, and we find that in its presence the transition may become continuous even at small values of $\epsilon$. The resulting quantum critical behavior in this case is discussed in some detail.
We also present an argument that in the limit of a large number of fermion species the phase transition becomes continuous only for $\epsilon >1$,
even if higher-loop terms are included.

The paper is organized as follows:
In Sec.~\ref{sec:model}, we introduce a continuum model for Dirac fermions coupled to a $U(1)$ gauge field.
We also discuss the parametrization of the Kekul\'e order parameter and its coupling to the fermionic sector.
In Sec.~\ref{sec:rge}, the renormalization group equations for the full model are presented and their fixed-point solutions are studied in Sec.~\ref{sec:FPana}.
Conclusions are drawn in Sec.~\ref{sec:conc}.

\section{Model}\label{sec:model}

The motion of non-interacting fermions on the honeycomb lattice can be most simply described in terms of a tight-binding Hamiltonian
\begin{align}\label{eq:h0}
H_0=-t\sum_{\vec{R},i}\left[u^\dagger(\vec{R}) v(\vec{R}+\vec{\delta}_i)+\text{h.c.}\right]\,,
\end{align}
with nearest-neighbor hopping amplitude $t$.
Here, $u$ and $v$ are the electron annihilation operators
at the two triangular sublattices of the lattice and the sum runs over the sites $\vec{R}$ of the first triangular sublattice with position vectors $\vec{R}_1^T=a\left(\sqrt{3}/2,-1/2\right)$ and $\vec{R}_2^T=a\left(0,1\right)$.
The lattice spacing $a$ is set to $a=1$ and the three nearest-neighbor vectors $\vec{\delta}_i$, located on the second sublattice, explicitly read $\vec{\delta}_1^T=\left(1/(2\sqrt{3}),1/2\right)$, $\vec{\delta}_2^T=\left(1/(2\sqrt{3}),-1/2\right)$ and $\vec{\delta}_3^T=\left(-1/\sqrt{3},0\right)$.
Diagonalization of the free Hamiltonian $H_0$ provides a spectrum with two energy bands and band dispersion $\epsilon_{\vec{k}}=\pm |\sum_{i=1}^3\exp(\vec{k}\cdot\vec{\delta}_i)|$.
At the corners of the Brillouin zone, given by the two points $\vec{K}^T=(2\pi/\sqrt{3},2\pi/3)$ and $\vec{K}^\prime=-\vec{K}$, the two energy bands touch linearly and isotropically. Here, the energy dispersion gives rise to two inequivalent Dirac points.

Retaining only the Fourier modes near the band touching points a continuum low-energy effective theory for $H_0$ can be written down in terms of a free Dirac Lagrangian\cite{Semenoff1984}, $\mathcal{L}_{\psi,0}=\bar\psi\,\gamma_\mu\partial_\mu\psi$,
where $\partial_\mu=(\partial_\tau, \vec{\nabla})$ and the $\gamma$~matrices are defined in a four-dimensional representation as $\gamma_0=\mathbbm{1}\otimes\sigma_z$, $\gamma_1=\sigma_z\otimes\sigma_y$, and $\gamma_2=\mathbbm{1}\otimes\sigma_x$ in $D=2+1$ dimensions.
The $\sigma_i$ are the conventional Pauli matrices and the conjugate of the Dirac field is given by $\bar\psi=\psi^\dagger\gamma_0$. We consider the Dirac field in Fourier space $\psi(x)=\int d^Dqe^{iqx}\psi(q)$ and express it in terms of Grassmann fields $u,v$ by\cite{herbut2006}
$\psi^\dagger(q)=\left[u^\dagger(K+q),v^\dagger(K+q),u^\dagger(-K+q),v^\dagger(-K+q)\right]$ with momentum vectors $q=(\omega,\vec{q})$ gathering Matsubara frequency $\omega$ and wavevector $\vec{q}$.
The reference frame is chosen such that $q_x=\vec{q}\cdot \vec{K}/|\vec{K}|$ and $q_y=(\vec{K}\times\vec{q})\times \vec{K}/|\vec{K}|^2$.
We define two more $\gamma$ matrices which anticommute with all $\gamma_\mu$, reading $\gamma_3=\sigma_x\otimes\sigma_y$ and $\gamma_5=\sigma_y\otimes\sigma_y$. We define $\gamma_{35}=-i\gamma_3\gamma_5$, which commutes with all $\gamma_\mu$ and anticommutes with $\gamma_3$ and $\gamma_5$.
As a convenient generalization of this model, we introduce an arbitrary number of pairs of Dirac points $\pm\vec{K}_i$ into the spectrum with $i=1,...,N$ and refer to $N$ as the number of fermion flavors of four-component spinors.
We note that the case of spin-1/2 fermions on the honeycomb lattice with spin-rotation invariance as, e.g., relevant for electrons in graphene, can be described by setting the number of fermion flavors to $N=2$ corresponding to an eight-component spinor.

We also introduce a coupling of the fermions to an electromagnetic $U(1)$ gauge field $A_\mu$.
This is formally done by minimal coupling, leading to the Lagrangian\cite{Roy2013c}
\begin{align}
	\mathcal{L}_\psi&=\bar\psi\,\gamma_\mu\left(\partial_\mu-ieA_\mu\right)\psi\,,
\end{align}
with the $U(1)$ charge coupling $e$, which describes the interaction of fermions with the electromagnetic field in terms of a local quantum field theory.

The kinetic term of the gauge field reads
\begin{align}
\mathcal{L}_{\text{em}}=\frac{1}{4}F_{\mu\nu}F_{\mu\nu}\,,
\end{align}
with $F_{\mu\nu}=\partial_\mu A_\nu-\partial_\nu A_\mu$ and we further add a gauge fixing term to the Lagrangian $\mathcal{L}_{\text{gf}}=\frac{1}{2\xi}(\partial_\mu A_\mu)^2$ with gauge fixing parameter $\xi$. Here, $\xi=0$ corresponds to Landau gauge and $\xi=1$ corresponds to Feynman gauge.
In the following, we will discuss two distinct cases\cite{gonzalez1993,Gorbar2001,Herbut2001,Herbut2013}: (1) The electromagnetic field propagates in three-dimensional space, i.e. in 3+1 dimensions whereas the propagation of the fermions is restricted to a lower dimension, e.g., $D=2+1$. (2) The electromagnetic field is defined such that it propagates in the same spacetime dimensions as the fermions, i.e. $D=d+1$.

\subsection{Kekul\'e valence bond solid}

The Kekul\'e VBS corresponds to a particular dimerization pattern of the fermions on the honeycomb lattice that can be captured by adding a bond-dependent modification of the hopping strength~\cite{hou2007} to the simple tight-binding Hamiltonian
\begin{align}\label{eq:Deltah0}
\Delta H_0=-\sum_{\vec{R},i}\left[\Delta t_{\vec{R},i}u^\dagger(\vec{R}) v(\vec{R}+\vec{\delta}_i)+\text{h.c.}\right]
\end{align}
where $\Delta t_{\vec{R},i}=\phi(\vec{R}) e^{i \vec{K}\cdot\vec{\delta}_i}e^{i (\vec{K}-\vec{K}^\prime)\cdot\vec{R}}/3+\text{c.c.}$ and $\phi$ is a complex-valued order parameter possessing $Z_3$ symmetry.
This modified hopping introduces a chiral mixing between the two Dirac points and gives rise to a phase transition from a Dirac semimetal towards a Kekul\'e VBS with a single-particle mass gap.
In terms of the low-energy effective field theory, the coupling between the fermions and the order parameter field reads\cite{roy2010c}
\begin{align}\label{eq:lagb}
	\mathcal{L}_y=\,y\left(\phi_1\bar\psi\, i\gamma_3\psi+\phi_2\bar\psi\, i\gamma_5\psi\right)\,,
\end{align}
with $\phi=\frac{1}{\sqrt{2}}(\phi_1+i\phi_2)$ and $\phi_1,\phi_2$ are real fields.
The dynamics of the bosonic order parameter is described by the Lagrangian
\begin{align}\label{eq:lagb}
	\mathcal{L}_\phi=\phi^\ast\left(\partial_\tau^2-\nabla^2+ m^2\right)\phi+ g\left(\phi^3+\phi^{\ast 3}\right)+\lambda|\phi|^4\,
\end{align}
which due to $Z_3$ symmetry allows the cubic term $\propto g$.
We note that the Kekul\'e order parameter is gauge neutral, and therefore its Lagrangian is not modified by the presence of the $U(1)$ gauge field.
This is in contrast to a superconducting order parameter, which combines two electrons into a Cooper pair and therefore comes with a charge\cite{Roy2013c,Zerf2016} of $2e$.

In sum, we study a $U(1)$ gauge theory for $N$ flavors of Dirac fermions in $D=d+1$ dimensions coupled to a charge-neutral complex order parameter corresponding to a Kekul\'e distortion.
The full Lagrangian $\mathcal{L}$ is composed of a free fermionic part $\mathcal{L}_\psi$, a bosonic part $\mathcal{L}_\phi$, a Yukawa interaction $\mathcal{L}_y$, the Maxwell Lagrangian $\mathcal{L}_{\text{em}}$ and the gauge fixing $\mathcal{L}_{\text{gf}}$, i.e. the action is given by
\begin{align}\label{eq:action}
S=\int d\tau d\vec{x}^{D-1}\left[\mathcal{L}_\psi+\mathcal{L}_\phi+\mathcal{L}_y+\mathcal{L}_{\text{em}}+\mathcal{L}_{\text{gf}}\right]\,.
\end{align}

\section{Renormalization group equations}\label{sec:rge}

A continuous phase transition and the concomitant (quantum) critical behavior is related to the presence of an infrared attractive fixed point in the system's renormalization group equations.
At a RG fixed point the system becomes scale invariant and the free energy exhibits a scaling form giving rise to universal critical exponents\cite{herbutbook}.
To approach the critical point of a phase transition one parameter has to be tuned to its critical value, e.g., the temperature $T\to T_c$ or, in the case of a quantum critical point, the doping level or a coupling strength.
In the RG framework a tuning parameter corresponds to an RG relevant direction, which -- unless precisely tuned -- drives the system away from its fixed point.

Canonical power counting determines the RG relevance of the parameters of a model close to its non-interacting limit, i.e. near the Gaussian fixed point.
At an non-Gaussian fixed point, however, the canonical power counting receives corrections from the interactions which can affect the relevance of the RG directions.
A positive (negative) power counting dimension suggests a relevant (irrelevant) parameter.
Parameters with vanishing power counting dimension are called marginal.
For instance the power counting dimension of the mass parameter in the Lagrangian $\mathcal{L}_\phi$ is $[m^2]=2$, i.e. it is strongly relevant and therefore needs to be fine-tuned to approach the critical point.
The power counting dimension of the quartic coupling $\lambda$ is $[\lambda]=\epsilon$ in $D=4-\epsilon$ dimensions.
It is only slightly positive and is turned negative at a non-Gaussian fixed point even if the interactions are perturbatively small and therefore this coupling automatically approaches its fixed point value toward the infrared.
The same reasoning also holds for the Yukawa coupling $y$ and the electromagnetic coupling $e$.
In most condensed matter models of similar type, this exhausts the available RG parameters leaving only the squared mass as an RG relevant direction, i.e. a single tuning parameter to approach the quantum critical point.

Here, however, the cubic order parameter terms introduce another RG relevant direction with a large power-counting dimension $[g]=1+\epsilon/2$.
In case only one tuning parameter is available, no scaling behavior, as induced by the vicinity to a RG fixed point, will be observed.
On the other hand, the quantum critical behavior of low-dimensional Dirac fermions coupled to order parameters is governed by non-trivial fixed points where the canonical power counting is modified.
When these modifications are large enough, they can flip the sign of the canonically relevant coupling $g$ and therefore restore the continuous transition.
We investigate this scenario and analyze its domain of applicability and control.

For  the  above  action $S$ and  at  zero  temperature we calculate the renormalization group (RG) equations in the Wilsonian scheme by simultaneously integrating out the fermionic as well as the bosonic modes within the narrow momentum shell $\Lambda/b< (\omega^2+\vec{q}^2)<\Lambda$.
At one-loop order in $D=4-\epsilon$ dimensions RG equations for the squared cubic and the quartic bosonic couplings appearing in $\mathcal{L}_\phi$, at $m=0$ read
\begin{align}\label{eq:beta01}
	\beta_{g^2}&=\frac{d g^2}{d \ln b}=(2+\epsilon-3Ny^2-6\lambda-\frac{27}{2}g^2)g^2\,,\\
	\beta_\lambda&=\epsilon\lambda-2Ny^2(\lambda-y^2)-5\lambda^2-162g^4+63\lambda g^2\,.\label{eq:beta02}
\end{align}
Here, we have rescaled the couplings $g^2\to g^2/(4\pi^2 \Lambda^{2+\epsilon})$ and $\lambda\to \lambda/(4\pi^2 \Lambda^{\epsilon})$.
The one-loop $\beta$ functions of the purely bosonic couplings are independent of the charge $e$ as it does not directly couple to the neutral boson field.
The squared Yukawa coupling $y^2$ flows as
\begin{align}\label{eq:betay}
	\beta_{y^2}&=\frac{d y^2}{d \ln b}=\epsilon y^2-(N+1)y^4-\frac{9}{2}g^2y^2+6e^2y^2\,,
\end{align}
where we have rescaled $y^2\to y^2/(4\pi^2 \Lambda^{\epsilon})$.
Here, we obtain contributions to the $\beta$ function from the $U(1)$ field due to a non-trivial renormalization of the Yukawa vertex, as reflected in the last term in Eq.~\eqref{eq:betay}.

Finally, we consider the $\beta$ function of the charge parameter for the case of a gauge field propagating in $3+1$ dimensions whereas fermions and complex bosons propagate on a lower-dimensional brane with $(3-\epsilon) +1$ dimensions~\cite{gonzalez1993,Gorbar2001,Herbut2001,Herbut2013}.
In this case, the bare gauge field propagator of the effective $(3-\epsilon) +1$-dimensional theory is obtained from ``integrating out" the additional spatial direction, which is most conveniently performed in Feynman gauge which leaves the effective gauge parameter invariant when going from 3+1 to $(3-\epsilon) +1$  dimensions. As a result one obtains a bare Maxwell term which is non-analytic in momentum $\propto |q|^{2-\epsilon}$, instead of the conventional, analytic,  $q^2$~dependence\cite{Gorbar2001}. Due to the resulting non-analyticity of the effective gauge field propagator at $q=0$ the charge coupling does not receive renormalization corrections and therefore becomes an exactly marginal coupling;  $\beta_{e^2}=0$. The $\beta$ function for $y^2$, on the other hand, remains unaffected by the above consideration.

We will also consider the $\beta$ function of the charge for the case that the $U(1)$ gauge field propagates in the same dimension as the fermions. In $D=3+1-\epsilon$ and with $e^2\to e^2/(4\pi^2 \Lambda^{\epsilon})$ it is then given by
\begin{align}
	\beta_{e^2}=\frac{d e^2}{d \ln b}=\epsilon e^2-\frac{4}{3}N e^4\,.\label{eq:betae}
\end{align}

\section{Fixed point analysis}\label{sec:FPana}

The set of the four $\beta$ functions, i.e. $\beta_{e^2},\beta_{\lambda},\beta_{y^2}$ and $\beta_{g^2}$, features various fixed points, e.g., the Gaussian (G) fixed point with $e^{2}_\ast =y^{2}_\ast=g^{2}_\ast=\lambda_\ast=0$ and the purely bosonic Wilson-Fisher (WF) fixed point with  $e^{2}_\ast =y^{2}_\ast=g^2_\ast =0\,,\ \lambda_\ast=\epsilon/5$.
In particular, we note that the $\beta$ function for the cubic coupling, Eq.~\eqref{eq:beta01}, always has the trivial fixed-point solution, $g^{2}_\ast=0$, even if the other couplings have non-vanishing fixed-point values.
In case $g^2=0$, a finite value for $g^2$ cannot be generated by the RG flow as the $Z_3$ symmetry of the model is replaced by a larger $U(1)$ symmetry, which is preserved under RG transformations. Here, we will discuss exclusively the Gross-Neveu-Yukawa fixed-points with $g^{2}_\ast=0$ that govern the quantum critical behavior of the model.

\subsection{Exactly marginal charge coupling}

We consider first the scenario where the charge is exactly marginal, and therefore can be taken as a fixed parameter of the model.
The fixed-point equation for the Yukawa coupling then reads
\begin{align}
y^{2}_\ast=\frac{\epsilon}{N+1}+\frac{6\bar e^2}{N+1}
\end{align}
where we have split the expression into a term without the gauge-field coupling $\propto \epsilon$ and a term exclusively originating from a non-vanishing charge.
Accordingly, Eq.~\eqref{eq:beta02} for the four-boson coupling yields a fixed-point value reading
\begin{align}
\lambda_{\ast}=\frac{1-N+R_{N}}{10(N+1)}\epsilon+\frac{6(S_{N}-N)}{5(N+1)}\bar e^2+\mathcal{O}(\epsilon^2,\bar e^4)\,,
\end{align}
to first order in $\epsilon$ with $R_{N}=\sqrt{1+38N+N^2}$ and $S_{N}=\sqrt{N^2+10N}$.
Here, we have also expanded around vanishing $\bar e^2$ as we consider the charge to be a small perturbation.
We note that the limits $\bar e^2\to 0$ and $\epsilon \to 0$ have to be taken with some care, as they do not commute.
We have chosen the order of limits in such a way as to recover best the behavior for $\epsilon \in [0,1]$ and small $\bar e^2 \lesssim 0.12$.
In the limit $\bar e^2 \to0$ our results for $y^2_\ast$ and $\lambda_\ast$ exactly recover the fixed-point values of the neutral model.
The RG scaling related to the interaction parameters $y^2$ and $\lambda$ can be derived from diagonalization of the corresponding block of the stability matrix
\begin{align}\label{eq:stabmat}
	\mathcal{M}=
	\begin{pmatrix}
	\frac{\partial\beta_{y^2}}{\partial y^2} & \frac{\partial\beta_{\lambda}}{\partial y^2}\\
	\frac{\partial\beta_{y^2}}{\partial \lambda} & \frac{\partial\beta_{\lambda}}{\partial \lambda}
	\end{pmatrix}\bigg|_{\bar e^2,y^2_\ast,\lambda_\ast,g_\ast^2}\,,
\end{align}
and we find the two negative eigenvalues
\begin{align}
m_1&=-\epsilon-6 \bar e^2\,,\quad m_2= -\frac{R_N}{N+1}\epsilon-\frac{12S_N}{N+1}\bar e^2 \,.
\end{align}

The RG scaling of the squared mass of the order parameter constitutes the first RG relevant direction of the model and determines the inverse correlation length exponent at the fermionic fixed point.
It is given by the relation
\begin{align}\label{eq:num1}
\theta_{m^2}=\nu^{-1}=2-Ny^{2}_\ast -2\lambda_{\ast}\,.
\end{align}
As the Yukawa coupling can become large for sizable charge this can reduce considerably the value of the RG scaling of the squared mass as compared to its canonical value of 2. We note that $\theta_{m^2}$ remains positive for $\bar e^2 \lesssim 0.12$ for all values of $N$ and $D\in [3,4]$.
Therefore, the bosonic mass parameter serves as a relevant RG direction providing the tuning parameter for the quantum phase transition.

\subsubsection{RG scaling of the cubic coupling}

The $\beta$ function of $g^2$ is proportional to $g^2$ itself
and so the RG scaling due to the cubic term at $g_\ast =0$ is given directly by the derivative of $\beta_{g^2}$ with respect to $g^2$ evaluated at the fermionic fixed point
\begin{align}\label{eq:theta2e2b}
	\theta=\frac{\partial \beta_{g^2}}{\partial g^2}\Big|_{\bar e^2,y^2_\ast,\lambda_\ast,g_\ast^2}\,.
\end{align}
The phase transition to the Kekul\'e VBS state can be rendered continuous when $\theta$ becomes negative. We calculate $\theta=\theta(\bar e^2)$ to order $\epsilon$,
\begin{align}
 \theta&=2+(1-3A_N)\epsilon-18\bar e^2  B_N\label{eq:theta2bare}\,,
\end{align}
where we have defined
\begin{align}
A_N&=\frac{1+4 N+R_N}{5 (N+1)}\,,\quad B_N=\frac{3N+2S_N}{5(N+1)}\,.
\end{align}
The function $A_N$ is positive, and in particular, $A_N\geq 1$ for $N\geq1/2$. To the leading order, the expression in Eq.~\eqref{eq:theta2bare} may become negative then with the increase of $\epsilon$, allowing the transition to become continuous.

The case $\bar e^2=0$, deserves special consideration: For $N<1/2$, we find $A_N<1$ and therefore $\theta(\bar e^2=0)$ is positive for any $\epsilon \leq 1$.
On the other hand, at $N=1/2$, we have $\theta(\bar e^2=0)=2-2\epsilon$, rendering $\theta=0$  for $\epsilon=1$.
Furthermore, at $\epsilon=1$, we find that for $N\geq 1/2$ the exponent $\theta$ is negative.
In the limit of large $N$, on the other hand, $\theta(e^2=0,N\to \infty)=2-2\epsilon$.
Fig.~\ref{fig:theta2e} shows that a negative $\theta(e^2=0)$ generally requires quite a large value of $\epsilon$ for all values of $N$.
For such large values of $\epsilon$ obviously it becomes questionable whether the conclusion about the relevance of the cubic coupling based on the leading order calculation is trustworthy.

For a finite charge one needs to include the function $B(N)$,  which is always positive, with $B_N\to 0$ for $N\to 0$. We find $B_{N=1}\approx 0.96$, $B_{N=2}\approx 1.05$ and $B_{N\to\infty}=1$, for example. A finite charge $\bar e^2\neq 0$ therefore only helps the tendency of rendering the cubic term irrelevant by gapless fermions. This is clearly exhibited in Fig.~\ref{fig:theta2e} for several specific values of $\bar e^2$, where one sees how the quantum critical region can be brought considerably closer to the perturbatively accessible regime of small~$\epsilon$.
A special limit is given in the case $N\to \infty$, where $A_{N\to\infty}=B_{N\to\infty}=1$. Here, for the choice $\bar e^2=1/9$, we have $\theta=-2\epsilon$ and $\theta_{m^2}=4/3-\epsilon$.

\begin{figure}[t!]
\includegraphics[width=1.0\columnwidth]{fig02_Nfcrit.pdf}
\caption{Sign of the dimension of the cubic coupling as a function of the dimensionality $D$ and the number of fermion flavors $N$ in first-order $\epsilon$ expansion with a fixed value of the charge $e^2$. The gray-shaded area shows the region in $(D,N)$ space where $\theta(\bar e^2=0)$ is negative, and the cubic term irrelevant. Quite a large values for $\epsilon$ are evidently needed in this case.
The behavior for finite values of $\bar e^2 \in \{1/100,1/20,1/10,1/9\}$ is indicated by the red lines. To the left of these lines, we have $\theta<0$ and the Kekul\'e transition being continuous. Finite charge $\bar e^2$ substantially decreases the required value of $\epsilon$ that would render $\theta<0$.}
\label{fig:theta2e}
\end{figure}

\subsection{Two-loop analysis for $N\to\infty$ for $\bar e^2=0$}

At two-loop order and at $\bar e^2=0$, the $\beta$ function for the cubic coupling squared  is composed from its one-loop and two-loop contributions, $\beta_{g^2}=\beta_{g^2}^{\text{(1-loop)}}+ \beta_{g^2}^{\text{(2-loop)}}$.
In the large-$N$ limit, we can deduce the exponent $\theta$ at two-loop order without accessing the full set of $\beta$ functions with all numerical coefficients at that order.
Recall that the fixed point value of $g_{\ast}^2=0$ is protected by symmetry to arbitrary order.
Due to the vertex structure of the model, all diagrams contributing to $\beta_{g^2}$ at two-loop order are exactly cubic in the couplings $y^2,g^2, \lambda$, e.g., they are porportional to $(g^2)^3, (g^2)^2 \lambda$, $(g^2) (y^2) \lambda,...\,$.
However, no diagrams that are of order $(g^2)^0$ appear, as these would generate a finite $g^2$ from $g^2=0$.

Diagrams that are of order $\mathcal{O}\left((g^2)^2\right)$ and higher cannot contribute to $\theta$, cf. Eq.~\eqref{eq:theta2e2b}, as their contributions vanish at the fixed point where $g_\ast^2=0$.
Therefore, the only two-loop diagrams in $\beta_{g^2}$ from the class of cubic terms in $y^2,g^2, \lambda$ that can contribute to $\theta$ are the ones that are linear in $g^2$, i.e.,
\begin{align}\label{eq:fast2l}
\hspace{-0.18cm}\Delta\beta_{g^2}^{\text{(2-loop)}}=g^2 \left(a\lambda^2+(b +c N) \lambda y^2+ (d+eN) y^4\right),
\end{align}
with $\{a,b,c,d,e\} \in \mathbbm{R}$.
Here, we have already introduced the appropriate factors of $N$ that can be deduced from simple diagrammatic considerations, i.e. appearing closed fermion loops in the two-loop diagrams, cf.~Fig.~\ref{fig:2loop}.

\begin{figure}[b!]
\includegraphics[width=.7\columnwidth]{fig03_2loop.pdf}
\caption{Diagrams contributing to the two-loop $\beta$ function of the cubic coupling. Here, the dashed lines represent bosons and the solid lines are fermions. The closed fermion loops come with a factor $N$. We do not display two-loop diagrams without $N$ multiplicity.}
\label{fig:2loop}
\end{figure}

The fixed point coordinates of the couplings $\lambda, y^2$ to two-loop order are identical to the ones of the uncharged chiral XY model without the cubic term~\cite{Rosenstein:1993zf}.
The large-$N$ behavior of the fixed point coordinates can be inferred from a canonical consideration of the two-loop contributions to the $\beta$ functions of  $\lambda$ and $y^2$,
\begin{align}
	\Delta\beta_{y^2}^{\text{(2-loop)}}&= a_y\lambda^2 y^2+b_y \lambda y^4+(c_y+N d_y)y^6\,,\\
	\Delta\beta_{\lambda}^{\text{(2-loop)}}&= a_\lambda\lambda^3 +N b_\lambda\lambda^2 y^2+N c_\lambda \lambda y^4+N d_\lambda y^6,\nonumber
\end{align}
which add to the corresponding one-loop expressions.
Again, we have introduced the unspecified coefficients $a_i,b_i,c_i \in\mathbbm{R}$ and appropriate factors of $N$.
With these structural two-loop expressions, the behavior of the order $\epsilon^2$ term of the GNY fixed point values for $y^2$ and $\lambda$ is determined to scale as $y^{2}_\ast(\epsilon^2)\sim\epsilon^2/N^2$ and $\lambda_\ast(\epsilon^2)\sim \epsilon^2/N^2$ for large $N$ which is consistent with Ref.~\onlinecite{Rosenstein:1993zf}.
Therefore, together with the one-loop contributions to $y^2_\ast$ and $\lambda_\ast$, we find that the derivative of Eq.~\eqref{eq:fast2l} w.r.t. $g^2$ scales as $1/N$ and does not provide a non-trivial contribution to $\theta$ in the limit $N\to \infty$.

The only potentially non-trivial term that is left to discuss at order $\epsilon^2$ is the contribution from the term $\sim~\hspace{-0.1cm}3N y^2_\ast$ in Eq.~\eqref{eq:beta01}.
This contribution, however, is also suppressed at order $\epsilon^2$ due to the scaling $y^{2}_\ast(\epsilon^2)\sim\epsilon^2/N^2$, i.e.
\begin{align}
	\theta^{\text{(2-loop)}}(N\to\infty)=\theta^{\text{(1-loop)}}(N\to\infty)=2-2\epsilon\,.
\end{align}
We conjecture this argument to perpetuate to higher loop orders as the introduction of any closed fermion loop coming with a factor $N$ comes with two Yukawa vertices and introduces another factor of $y^{2}_\ast\sim 1/N$ at the fixed point, suggesting $\theta_2^{\text{(exact)}}(N\to\infty)=2-2\epsilon$.

\subsection{Gauge field in $(4-\epsilon)$ dimensions}\label{sec:charge}

In this section, we study the case where the gauge field propagates in the same number of dimensions as the fermions and complex boson fields.
Then, the $\beta$ function for the charge is given by Eq.~\eqref{eq:betae} and charge parameter $e^2$ features a non-trivial fixed point,
\begin{align}\label{eq:east}
e^{2}_\ast =\frac{3}{4 N}\epsilon\,.
\end{align}
Here, the charge corresponds to an RG irrelevant direction, as $(\partial\beta_{e^2}/\partial e^2)|_{e^2=e^{\ast 2}}=-\epsilon$.
We deduce the fixed-point value of the Yukawa coupling, Eqs.~\eqref{eq:betay} and \eqref{eq:east},
\begin{align}\label{eq:fpymode}
	y^{2}_\ast =\frac{1+6 (e^{2}_\ast/\epsilon)}{N+1}\epsilon=\left(1+\frac{9}{2N}\right)\frac{\epsilon}{N+1}\,,
\end{align}
which then allows to calculate the fixed-point value of the quartic coupling
\begin{align}\label{eq:lambdaast}
	\lambda_\ast=\frac{C_N-8-N}{10(1+N)}\epsilon\,,
\end{align}
where $C_N=\sqrt{N^2+56N+424+810/N}$ and $\lambda_\ast>0$ for all $N$.
For large $N$, $\lambda_\ast$ vanishes like $\lambda_\ast\sim 2\epsilon/N$.
The RG scaling related to the interaction parameters $y^2$ and $\lambda$ can be derived from the stability matrix block, Eq.~\eqref{eq:stabmat}, where $\bar e^2$ is replaced by $e_\ast^2$.
This provides the eigenvalues $m_1=-\epsilon(1+9/(2N))$ and $m_2=-\epsilon\, C_N/(1+N)$ which are negative for all $N$.

\begin{figure}[t!]
\includegraphics[width=.9\columnwidth]{fig04a_Nfcrit.pdf}
\includegraphics[width=.9\columnwidth]{fig04b_Nfcrit.pdf}
\caption{Sign of the critical exponent resulting from the presence of the cubic coupling as a function of the dimensionality $D$ and the number of fermion flavors $N_f$ in first-order $\epsilon$ expansion.
The gray-shaded areas show the regions, where $\theta_2<0$. The two shadings distinguish the situations, where the coupling to the gauge field vanishes $e=0$ (darker gray) and where it is at its own fixed point $e=e^\ast$ (lighter gray). Red red dotted line shows Eq.~\eqref{eq:nfcrit}.}
\label{fig:theta2}
\end{figure}

Accordingly, the RG scaling of the cubic term is given by
\begin{align}\label{eq:theta2e2}
	\theta=\frac{\partial \beta_{g^2}}{\partial g^2}\Big|_{e_\ast^2,y^2_\ast,\lambda_\ast,g_\ast^2}=2-\frac{77+14 N +6C_N}{10 (1+N)}\epsilon\,.
\end{align}
In Fig.~\ref{fig:theta2}, the region where $\theta$ becomes negative in this scenario, is indicated by the shaded region which is separated by the dashed line from the region where $\theta>0$.
As two interesting numerical examples, we explicitly note that $\theta(N=2,\epsilon)$ becomes negative for $\epsilon\gtrsim 0.21$ or $D \approx 3.79$ and $\theta(N=1,\epsilon)$ for $\epsilon\approx 0.13$ or $D=3.87$.

In the limit of large $N$, we again find $\theta(N\to\infty)\to 2-2\epsilon$, so a sign change of $\theta$ requires a rather large value for $\epsilon$.
In contrast, for small $N$ we can simplify Eq.~\eqref{eq:theta2e2} and calculate the critical number $N_c$ for a given $\epsilon$ where $\theta$ becomes negative,
\begin{align}\label{eq:nfcrit}
	\theta\approx 2-\frac{27}{5}\sqrt{\frac{10}{N}}\epsilon=0\Rightarrow N_{c}\approx\frac{729}{10}\epsilon^2\,.
\end{align}
For small $N$ the position where $\theta$ changes sign, converges to  $N\to 0$ as $D$ approaches 4, see the red dotted line in Fig.~\ref{fig:theta2}.
This limit, however, has to be considered with some care, as for a given $\epsilon$ the fixed point values can grow large when $N$ decreases.
Therefore, the position where the sign change appears, e.g., in terms of dimensionality $D$, can be associated with large interaction parameters.
On the other hand, for all physically relevant $N$, i.e. $N\in \{1/2,1,2,...\}$ the fixed point values of all couplings remain of order $\mathcal{O}(1)$ in units of $\epsilon$.

\section{Summary \& Conclusions}\label{sec:conc}

We studied the quantum phase transition of $U(1)$ charged Dirac fermions coupled to a Kekul\'e order parameter with $Z_3$ symmetry.
In $4-\epsilon$ dimensions, the canonically relevant cubic term, which is allowed by the symmetry, can be rendered irrelevant through quantum fluctuations. When the charge coupling vanishes this occurs only for large values of $\epsilon$. A finite charge coupling, however, can support strongly the scenario where the renormalization group direction corresponding to the cubic coupling is rendered irrelevant.
In this case, the system gives rise to quantum critical behavior. We have investigated two scenarios:

(1) The charge is an exactly marginal coupling, as the gauge field propagates in four spacetime dimensions.
In that case, the value of the charge parameter $\bar e^2$ may be chosen freely, to control the position of the boundary where the critical exponent, determining the relevance of the cubic term, changes sign. Small values of the charge, $\bar e^2 \lesssim 1/10$, move this boundary from being close to three dimensions to  being very close to four (space-time) dimensions, bringing it under perturbative control.

(2) The charge coupling features a non-trivial fixed point value, when the  gauge field propagates in the same $4-\epsilon$ dimensions as the matter fields.
Here, the presence of a finite charge fixed-point coupling also strongly supports the emerging irrelevance of the cubic term at smaller values of $\epsilon$.
For a small number of fermion flavors, the boundary where the coupling is rendered irrelevant moves close to four dimensions, however, at the cost of the prefactors of the $\sim \epsilon$ fixed point values growing large.

Furthermore, we have carefully discussed the limit of a large number of Dirac fermion flavors in the case of vanishing charge, as well as for finite charge in the scenario~(2). In this limit, the cubic coupling always provides a renormalization group relevant direction unless $\epsilon>1$,  which is again beyond the domain of validity of this approach. Finally, for vanishing charge we have shown that this statement also holds at two-loop order, and provided an argument that suggests that this conclusion continues to hold at arbitrary order in the loop expansion.

\begin{acknowledgments}
The authors are grateful to Fakher Assaad and Igor Boettcher for useful discussions. M.M.S. is supported by DFG Grant No. SCHE 1855/1-1 and I. F. H. is supported by the NSERC of Canada.
\end{acknowledgments}


\section*{Methods}

\subsection{Details of the experiment}

The granular medium was made of a number $N \in [300,2600]$ of spheres
of diameter $d=4$ mm made of non-magnetic steel (mass of each sphere:
$0.267$ g), glass (mass $0.0854$ g), or
delrin\textsuperscript{\textregistered} (mass $0.0462$ g). They were
housed in a plexiglas\textsuperscript{\textregistered} cylinder with a
conical-shaped floor (diameter $90$ mm, minimum height $28,5$ mm,
maximum height $47.5$ mm) in which a plexiglas vane (height $15$ mm,
width $6$ mm, length $35$ mm) was suspended in order to be in contact
with the granular medium and not with the container~\cite{gnoli14}.  The container
was vertically vibrated by an electrodynamic shaker (LDS V450) fed
with an acceleration signal $a(t)$. In most of the experiments $a(t)$
is a white noise with a band-pass filter between $200$ Hz and $400$
Hz, while in the lowest $p_{00}$ case ($p_{00}=78$ Pa) we used a
sinusoidal signal at frequency $53$ Hz. This choice is motivated by
two empirical observations: 1) a lower number of particles (as in the
case of low $p_{00}$) requires a larger energy input to be
homogeneously fluidized and to reach the blade, and this can be
obtained by supplying energy through a sinusoidal signal
at low frequency; 2) in dense cases a sinusoidal signal induces
spurious resonances, while in diluted cases such resonances are never
observed. We have checked that performing the same experiments with
noise signal for $a(t)$ (pushing the shaker to its working limits)
gives flow curves with the same shape.  An accelerometer placed on the
container side measured $a(t)$, allowing us to define
$\Gamma=\sqrt{2\langle a^2(t)\rangle}/g$.
The vane, mounted through its rotation axis to a rotary encoder, was
also connected to a dc motor (typical working voltage $12$ V) as the
source of the driving torque. The motor was directly fed by a dc
voltage supply in the range $0$ to $7$ V. No limit was set for the
maximum current absorbed by the motor that, averaged on the duration
of the experiment, was never higher than $450$ mA. A data acquisition
system collected data for the angular position/velocity of the vane,
the effective motor voltage, the current circulating in the motor and
the root-mean-square vertical acceleration of the container. A
procedure of calibration allowed us to translate average values of
current into average values of applied torque. The same procedure
helped determining the moment of inertia of the rotating block,
$3.2\times 10^2 \textrm{g mm}^2$ (the blade with its axis and the
gears linking it to the motor). The typical experiment, at a given
$\Gamma$ and applied motor voltage, was $3600$ s long, with the
granular always ``reset'' at the beginning of each run for $30$ s at
high shaking intensity ($\Gamma=42$) and motor off. This procedure -
together with periodic replacement of used spheres - guaranteed
reproducible results at a distance of several weeks. Packing fractions
was non-homogeneous (it was larger in regions far from the borders of
the container): its value at rest was estimated to be in the range
$55\% - 70\%$, while it decreased when vibration was switched on. In
the analysis we have identified the shear rate $\dot\gamma$ with the
average of the angular velocity $\omega(t)$ of the rotating blade,
i.e. $\dot\gamma=\langle \omega(t) \rangle$, while the shear stress
$\sigma$ is proportional to the average of the applied torque $T(t)$
through the shear stress constant $\kappa$, i.e. $\sigma=\kappa
\langle T(t) \rangle$ with $\kappa=2 \pi R^2 H$ (with $R$ and $H$ the
blade half-length and height, respectively)~\cite{dzuy83}.

\subsection{Details of data analysis}

Velocity power density spectra (VPDS) are defined as $S(f)=(1/(2
t_{TOT}) |\int_0^{t_{TOT}} \omega(t) e^{i(2 \pi f) t} dt|^2$ with
$t_{TOT}$ the time-length of an experiment ($=3600$ seconds). Some
examples of $S(f)$ curves are shown in the Fig.~1 of Supplementary
Information (SI). In~\cite{scalliet15} VPDS in a similar vibro-fluidized
experimental setup, without applied torque ($\sigma=0$), has been
investigated. In the dilute or gas-like limit, e.g. low number of
spheres at high shaking, the VPDS takes a simple Lorentzian shape
$S(f)=D/[1+(2 \pi f /f_{visc})^2]$, with $D$ the asymptotic (long
time) diffusivity and $f_{visc}$ the effective viscosity due to
granular gas-vane collisions. When the number of particles (density)
is increased and/or the intensity of shaking ($\Gamma$) is reduced,
the system approaches a slow liquid regime and the VPDS develops a
wide bump (or smooth peak) with a maximum near $f \sim 20$ Hz, which
is associated to oscillations of the velocity autocorrelation induced
by liquid cages. At much smaller frequencies the VPDS reaches a
plateau whose height, $\lim_{f \to 0} S(f)$, corresponds to diffusivity
$D$: indeed the blade is not trapped in a cage forever, eventually it
manages to explore a much larger phase space and reaches normal
diffusion. From the low frequencies plateau of VPDS we have extracted
values of $D$ for Figure~\ref{fig:fluct}. We have defined the
cage-exit frequency $f_{cage}$ as the x-position, in the VPDS plot, of the
minimum separating the cage bump from the low-frequency diffusive
plateau (see filled circles in Fig.~1 of SI).

\vspace{.5cm}

{\bf Acknowledgements}

We thank Giorgio Pontuale for useful discussions and suggestions. We
also thank Rachele De Felice and Francisco Vega Reyes for reading
the manuscript.

\vspace{.5cm}

{\bf Corresponding author}

Correspondence to: Andrea Puglisi, andrea.puglisi@roma1.infn.it


\section{Introduction}

The discovery of graphene and other two-dimensional (2D) crystals in 2004 was a major scientific breakthrough in material science \cite{Novoselov2004, Novoselov2005}. In the past decade the family of 2D crystals has widened encompassing a large selection of compositions including almost all the elements of the periodic table \cite{Heine14}. The variety of their electronic properties is extremely broad, including metals and semi-metals (for example graphene \cite{Novoselov2004, CastroNeto09}), insulators (hexagonal boron nitride \cite{Blake2011}) and semiconductors (transition-metal dichalcogenides \cite{Heinz2010, Heinz2014}). These single-layer atomic crystals are stable under ambient conditions and appear continuous on a macroscopic scale  \cite{Novoselov2005}.

Macroscopic continuity is confirmed also by their optical properties \cite{Ramasubramaniam14, Nair2008, Blake2007, Heinz2014, Blake2011, Kravets2010}. In analogy to 3D materials where the susceptibility and the conductivity fix the optical response, single-layer atomic crystals are described as zero-thickness sheets with a surface susceptibility and a surface conductivity \cite{Pershoguba07, Hanson08, Merano15, Merano16}. As for bulk materials, ellipsometry is able to furnish both these parameters, showing that these are the macroscopic physical quantities experimentally accessible from the linear optical response of 2D atomic crystals \cite{Merano16}. Also the nonlinear optical properties are described by modelling them as zero-thickness sheets with a nonlinear surface susceptibility \cite{Zhao13, Heinz13, Merano216}.

A classical model has addressed recently the optical properties of a single-layer 2D crystal at the microscopic level \cite{Luca16}. This approach gave remarkable insight. It explained why the Fresnel coefficients for a single-layer 2D atomic crystal are intrinsically complex quantities even when a null surface conductivity is assumed. The crystal was modelled as atoms with polarizability $\alpha$ distributed on a Bravais lattice. It was shown that (contrary to the 3D case) \cite{Aspnes82, Wolf} all the dipoles contribute to the local electric field and that retardation effects (due to the finite velocity of propagation of the dipolar potential) dephase the local electric field from the electric field incident on the crystal \cite{Luca16}.

Here this classical approach is further extended to treat in a complete way the radiating properties of a 2D crystal. Improvements to ref. \cite{Luca16} must consider two directions that in the end will merge. The first: theory in ref. \cite{Luca16} provides an exact solution for the local electric field when an electromagnetic plane wave is incident on a 2D dipole lattice. This field acting on each dipole, is not equal to the macroscopic field in the 2D crystal. Anyway a complete theory should be able to connect the microscopic description with the macroscopic one. The second: the radiating dipoles in a 2D crystal should experience the radiative reaction force. In \cite{Luca16} only the retarded potentials have been considered while radiative reaction processes are better taken into account by addressing also the role of the advanced potentials \cite{Page24, Dirac38, Wheeler45, Panofsky}.

\section{The mechanism of radiation reaction}

The origin of the force of radiative reaction has puzzled scientists for a long time. Use of action at a distance with field theory as equivalent and complementary tools for the description of nature has for long been prevented by inability of the first point of view fully to account for the mechanism of radiation. Elucidation of this point came with a combined effort of many physicists, among the others: Lorentz \cite{Panofsky}, Page \cite{Page24}, Dirac \cite{Dirac38}, Wheeler and Feynman \cite{Wheeler45}.

For the purpose of this paper it is interesting to refer back to the calculation of radiative reaction made by Lorentz on an extended charge, every part of which exerted a retarded effect upon every other part \cite{Panofsky}. This calculation was further examined in \cite{Page24, Dirac38, Wheeler45} to produce a well-defined and relativistically invariant prescription to compute the magnitude of the radiative damping. If we consider a system of $n$ charged particles, the electromagnetic fields which act on a given particle $j$ arise only from other particles \cite{Wheeler45}. These fields are represented by one-half the retarded plus one-half the advanced Lienard-Wiechert solutions of Maxwell's equations \cite{Dirac38, Wheeler45, Panofsky}. The Lorentz's radiation-reaction can be defined as the difference between one-half the retarded fields and one-half the advanced fields generated by the other particles \cite{Page24, Dirac38, Wheeler45}. This field compensates the half-advanced fields and combines with its half-retarded fields, to produce the full retarded disturbance $\vec{\textbf{\emph{E}}}_{j}$ which is required by experience \cite{Wheeler45}
\begin{eqnarray}
\label{Retarded field}
\vec{\textbf{\emph{E}}}_{j}&=&\frac{1}{2}{\sum_{n}}'\bigg(\vec{\textbf{\emph{E}}}_{ret, n}+\vec{\textbf{\emph{E}}}_{adv, n}\bigg)  \nonumber \\
&+&\frac{1}{2}{\sum_{n}}'\bigg(\vec{\textbf{\emph{E}}}_{ret, n}-\vec{\textbf{\emph{E}}}_{adv, n}\bigg)
\end{eqnarray}
where the dash in the summation indicates that the $j$ particle is excluded from the sum. This procedure is not convenient to apply to cases where relativistic contraction is so strong to modify the form of the crystal  \cite{Wheeler45}. These last cases will not be considered here.

\section{Radiation-reaction electric field in 2D crystals}

A 2D crystal can be treated as an extended charge distribution or better an extended electric dipoles distribution \cite{Luca16}.  I will find an expression for the radiative-reaction electric field $\vec{\textbf{\emph{E}}}_{R}$ using the above prescription that is well-defined and it involves no more than standard electromagnetic theory. An argument based on conservation of energy will further confirm that the proposed expression is the corrected one \cite{Jackson}. Finally, $\vec{\textbf{\emph{E}}}_{R}$ will be connected with the macroscopic electric field.

\begin{figure}
\includegraphics{Fig_1.eps}
\caption{The spherical waves depicted here represent light reemitted by the dipoles (black dots) in a 2D crystal that is exposed to a plane wave entering from the top (solid line). The reflected wave is traveling to the top. The transmitted wave is a coherent superposition of the spherical wavelets and the incident wave.}
\end{figure}

\subsection{Classical model of radiating 2D crystals}

I consider an insulating 2D crystal formed by atoms with polarizability $\alpha$ placed on a 2D Bravais lattice. A linearly polarized electromagnetic plane wave is incident on this 2D crystal (Fig.1). For simplicity normal incidence is assumed. The local electric field $\vec{\textbf{\emph{E}}}_{loc}$ acting on a single dipole (for instance the one in the origin) is given by
\begin{eqnarray}
\label{Local2}
\vec{\textbf{\emph{E}}}_{loc}e^{i(\omega t + \varphi)}=\vec{\textbf{\emph{E}}}_{i}e^{i\omega t}
+{\sum_{(m,n)}}'\vec{\textbf{\emph{E}}}_{n,m}(t)
\end{eqnarray}
where $\vec{\textbf{\emph{E}}}_{loc}$ and $\vec{\textbf{\emph{E}}}_{i}$ are real vectors, $\vec{\textbf{\emph{E}}}_{i}$ is the incident electric field and the last term is the full retarded disturbance which is required by experience \cite{Wheeler45}. This last term was computed in \cite{Luca16} for three different dipole distributions: square, triangular and honeycomb lattices. Formula (\ref{Local2}) evidences that the dipole fields dephase the local field with respect to the incident electric field. On the macroscopic scale this translates to a dephasing between the incident electric field and the surface polarization density $\vec{\textbf{\emph{P}}}(t)$
\begin{eqnarray}
\label{Local1}
\vec{\textbf{\emph{P}}}(t)=N\vec{\textbf{\emph{p}}}(t)=N\alpha \epsilon_{0} \vec{\textbf{\emph{E}}}_{loc}e^{i(\omega t + \varphi)}
\end{eqnarray}
where $\vec{\textbf{\emph{p}}}$ is the induced dipole moment at each reticular point, $N$ is the number of primitive cells per unit area and $\epsilon_{0}$ is the vacuum permittivity.

\subsection{Advanced and retarded electric fields for a dipole}

The theory developed in \cite{Luca16} used only retarded fields. As explained above a complete treatment of the radiating crystal requires also the advanced fields. Following ref. \cite{Feynman} the field generated by an electric dipole is
\begin{eqnarray}
\label{Dipole2}
\vec{\textbf{\emph{E}}}(\vec{\textbf{\emph{r}}}, t)=\frac{1}{4\pi \epsilon_{0}r^3} \biggl(3(\tilde{\vec{\textbf{\emph{p}}}}\cdot\hat{\textbf{\emph{r}}})\hat{\textbf{\emph{r}}}-\tilde{\vec{\textbf{\emph{p}}}} -\frac{( \vec{\textbf{\emph{r}}} \times \ddot{\vec{\textbf{\emph{p}}}})  \times \vec{\textbf{\emph{r}}}}{c^2}  \biggr)
\end{eqnarray}
where
\begin{eqnarray}
\vec{\textbf{\emph{p}}}&=&\vec{\textbf{\emph{p}}}(t \pm \frac{r}{c}) \nonumber \\
\tilde{\vec{\textbf{\emph{p}}}}&=&\vec{\textbf{\emph{p}}}(t \pm \frac{r}{c}) \mp \frac{r}{c}\dot{\vec{\textbf{\emph{p}}}}(t \pm \frac{r}{c})
\end{eqnarray}
here the $\pm$ sign in parentheses means that the dipole moment must be evaluated at the corresponding advanced or retarded time.
Assuming a temporal dependence $e^{i \omega t}$
\begin{eqnarray}
\vec{\textbf{\emph{p}}}=\vec{\textbf{\emph{p}}}_0 e^{i(\omega t \pm kr)}
\end{eqnarray}
I obtain:
\begin{eqnarray}
\label{Dipole3}
\vec{\textbf{\emph{E}}}(\vec{\textbf{\emph{r}}}, t)&=&\frac{e^{i(\omega t \pm kr)}}{4\pi \epsilon_{0}r^3} \biggl((3(\vec{\textbf{\emph{p}}}_0\cdot\hat{\textbf{\emph{r}}})\hat{\textbf{\emph{r}}}-\vec{\textbf{\emph{p}}}_0)(1\mp kr) \nonumber \\
&+&k^2( \vec{\textbf{\emph{r}}}\times \vec{\textbf{\emph{p}}}_0)  \times \vec{\textbf{\emph{r}}}  \biggr)
\end{eqnarray}
As for a plane wave, if the advanced field is
\begin{eqnarray}
\label{Advanced}
\vec{\textbf{\emph{E}}}_{adv}(\vec{\textbf{\emph{r}}}, t)= e^{i \omega t}\vec{\textbf{\emph{f}}}(\vec{\textbf{\emph{r}}})
\end{eqnarray}
the retarded one is
\begin{eqnarray}
\label{Retarded}
\vec{\textbf{\emph{E}}}_{ret}(\vec{\textbf{\emph{r}}}, t)= e^{i \omega t}\vec{\textbf{\emph{f}}}^*(\vec{\textbf{\emph{r}}})
\end{eqnarray}

\subsection{Expression of the radiation-reaction electric field}

I compute the radiation-reaction electric field acting on a dipole as half the difference between the retarded and the advanced electric fields generated by all the other dipoles.

\subsubsection{Square and triangular lattices}

For a square and a triangular lattice, from formula (36) of ref. \cite{Luca16} and formulas (\ref{Advanced}) and (\ref{Retarded}), I obtain
\begin{eqnarray}
\label{Reaction_square}
\vec{\textbf{\emph{E}}}_{R}(t) &= &{\sum_{(m,n)}}' \frac{1}{2}\bigg( \vec{\textbf{\emph{E}}}_{ret, \ m,n}(t)-\vec{\textbf{\emph{E}}}_{adv, \ m,n}(t) \bigg) \qquad \nonumber \\
&=&\frac{i\, \alpha C_1 k}{4\pi a^2} \vec{\textbf{\emph{E}}}_{loc}(t)=- \frac{\eta N}{2}  \dot{\vec{\textbf{\emph{p}}}}(t)= - \frac{\eta}{2}  \dot{\vec{\textbf{\emph{P}}}}(t)
\end{eqnarray}
where $a$ is the lattice spacing, $C_{1}= -2 \pi N a^2$ (formula (50) of ref. \cite{Luca16}), and $\eta$ is the wave impedance of vacuum.

\subsubsection{Honeycomb lattice}

I consider now the case of a honeycomb lattice with two atoms of different polarizability $\alpha_1$ and $\alpha_2$ per primitive cell. From formula (13) and (43) of ref. \cite{Luca16}  and formulas (\ref{Advanced}) and (\ref{Retarded}), I obtain
\begin{eqnarray}
\label{Reaction_honey}
\vec{\textbf{\emph{E}}}_{R}(t) &=& \frac{i\, C_1 k}{4\pi a^2} \bigg(\alpha_1 \vec{\textbf{\emph{E}}}^{(1)}_{loc}(t) + \alpha_2 \vec{\textbf{\emph{E}}}^{(2)}_{loc}(t) \bigg) \nonumber \\
&=&- \frac{\eta N}{2} \bigg( \dot{\vec{\textbf{\emph{p}}}}_1(t)+\dot{\vec{\textbf{\emph{p}}}}_2(t)  \bigg)= - \frac{\eta}{2}  \dot{\vec{\textbf{\emph{P}}}}(t)
\end{eqnarray}
where $\vec{\textbf{\emph{E}}}^{(1)}_{loc}(t)$ and $\vec{\textbf{\emph{E}}}^{(2)}_{loc}(t)$ are the local fields acting on the two atoms of the primitive cell. Formulas (\ref{Reaction_square}) and (\ref{Reaction_honey}) express the radiation-reaction electric field acting on the single dipoles in terms of their dipole moments. Being proportional to the negative time-derivative dipole moments, this force induces a damping and an energy loss due to the emission of radiation.

\section{Reflected and transmitted fields from a radiating 2D crystal}

Formulas (\ref{Reaction_square}) and (\ref{Reaction_honey}) are intriguing because they connect a field acting on a single dipole with a macroscopic quantity: the time derivative of the polarization density.  Wherever the polarization in matter changes with time there is an electric current $\vec{\textbf{J}}_{p}$, a genuine motion of charges. The connection between the radiation-reaction electric field and the current density is
\begin{eqnarray}
\label{current density}
\vec{\textbf{J}}_{p}=\dot{\vec{\textbf{\emph{P}}}}(t)= - \frac{2}{\eta}\vec{\textbf{\emph{E}}}_{R}(t)
\end {eqnarray}
This surface current density is a relevant quantity that determines the boundary conditions for the macroscopic fields and hence the Fresnel coefficients of the 2D crystal. The macroscopic electric field $\vec{\textbf{\emph{E}}}$ in the crystal is proportional to the surface polarization density
\begin{eqnarray}
\label{Polarization}
\vec{\textbf{\emph{P}}}(t)=N\alpha \epsilon_{0} \vec{\textbf{\emph{E}}}_{loc}(t)=\epsilon_{0} \chi \vec{\textbf{\emph{E}}}(t)
\end{eqnarray}
where $\chi$ is the surface susceptibility \cite{Merano16}. With reference to fig. 2, where a plane wave falls onto the 2D material from the top half space, the boundary conditions for the macroscopic fields are
\begin{eqnarray}
\label{Boundary}
\hat{\kappa} \wedge (\vec{\textbf{E}}_{2}-\vec{\textbf{E}}_{1})&=&0 \\
\label{Boundary2}
\hat{\kappa} \wedge (\vec{\textbf{H}}_{2}-\vec{\textbf{H}}_{1})&=&\textbf{J}_{p}
\end{eqnarray}
where $\hat{\kappa}$ is the unit vector along the $z$ axis and the subscripts 1 and 2 refer to the two half spaces separated by the crystal. Without loss of generality I restrict the treatment to an $s$ polarized incident wave.

\begin{figure}
\includegraphics{Fig_2.eps}
\caption{Electric field of an $s$ polarized plane wave normally incident on a single-layer two-dimensional crystal. The two half spaces separated by the atomically thin crystal are indicated as 1 and 2.}
\end{figure}

\subsubsection{Square and triangular lattices}

From equations (\ref{current density}-\ref{Boundary2}) I obtain:
\begin{eqnarray}
\label{Fresnel}
E_{xi}+E_{xr}&=&E_{xt} \\
E_{xi}+E_{xr}&=&\frac{N\alpha}{\chi}E_{loc\ x} \nonumber  \\
H_{yi}-H_{yr}&=&H_{yt}-\frac{2i\alpha C_{1}k}{4\pi a^2\eta}E_{loc\ x}= H_{yt}+\frac{i\alpha N k}{\eta}E_{loc\ x} \nonumber
\end{eqnarray}
where the subscript $i, r, t$ denote the incident, the reflected and the transmitted waves. Defining $r_{s}=E_{r}/E_{i}$ and $t_{s}=E_{t}/E_{i}$ as the reflection and the transmission coefficients, I have
\begin{eqnarray}
\label{Fresnel1}
r_{s}=\frac{-ik\chi}{ik\chi+2}; \qquad t_{s}=r_{s}+1
\end{eqnarray}

\subsubsection{Honeycomb lattice}

For the honeycomb lattice equations (\ref{current density}-\ref{Boundary2}) give:
\begin{eqnarray}
\label{Fresnel2}
E_{xi}+E_{xr}&=&E_{xt} \\
E_{xi}+E_{xr}&=&\frac{N}{\chi}(\alpha_1 E_{loc\ x}^{(1)} +\alpha_2 E_{loc\ x}^{(2)})\nonumber  \\
H_{yi}-H_{yr}&=& H_{yt}+\frac{i N k}{\eta}(\alpha_1 E_{loc\ x}^{(1)} +\alpha_2 E_{loc\ x}^{(2)}) \nonumber
\end{eqnarray}
providing the same Fresnel coefficients of the above case
\begin{eqnarray}
\label{Fresnel3}
r_{s}=\frac{-ik\chi}{ik\chi+2}; \qquad t_{s}=r_{s}+1
\end{eqnarray}
The validity of these expressions has already been proved by experiments \cite{Merano16}.

\section{Work done by the radiation-reaction force}

For the square and the triangular lattice, the time-average work ($W_R$) per unit time and per unit area done by the radiation-reaction electric field is given by $N$ times the one done on a dipole
\begin{eqnarray}
\label{time-average work}
W_R&=&\frac{N}{T}\int_{0}^{T} Re\bigg( \vec{\textbf{\emph{E}}}_{R}(t)\bigg) \cdot Re\bigg( \dot{\vec{\textbf{\emph{p}}}}(t)\bigg) dt \nonumber \\
&=&-\frac{R\vec{\textbf{\emph{E}}}^2_{i}}{\eta}
\end{eqnarray}
where $R$ is the reflectivity of the crystal \cite{Merano16}. For a honeycomb lattice the time-average work per unit time and per unit area done by the radiation-reaction electric field is given by $N$ times the one done on a primitive cell
\begin{eqnarray}
\label{time-average work honeycomb}
W_R&=&\frac{N}{T}\int_{0}^{T} Re\bigg( \vec{\textbf{\emph{E}}}_{R}(t)\bigg) \cdot Re\bigg( \dot{\vec{\textbf{\emph{p}}}}_1 (t)+\dot{\vec{\textbf{\emph{p}}}}_2 (t)\bigg) dt \nonumber \\
&=&-\frac{R\vec{\textbf{\emph{E}}}^2_{i}}{\eta}
\end{eqnarray}
As expected \cite{Jackson}, the expression standing on the right in formulas (\ref{time-average work}) and (\ref{time-average work honeycomb}) is negative and equal in absolute value to the time-average intensity radiated by the 2D crystal. Half of this electromagnetic intensity has an opposite direction to the incident field and it is the reflected field (Fig. 1). The other half travels in the same direction of the incident field and it recombined with it to create the transmitted electromagnetic field. Since I am considering an insulating crystal the  time-average work per unit time and unit area done by the radiation-reaction electric field is opposite to the one done by the incident field
\begin{eqnarray}
\label{time-average work incident}
W_i&=&\frac{N}{T}\int_{0}^{T} Re\bigg( \vec{\textbf{\emph{E}}}_{i}(t)\bigg) \cdot Re\bigg( \dot{\vec{\textbf{\emph{p}}}}(t)\bigg) dt \nonumber \\
&=&\frac{R\vec{\textbf{\emph{E}}}^2_{i}}{\eta}
\end{eqnarray}

\section{Connection between the radiation-reaction electric field and the macroscopic field}

I can now establish an important connection in between $\vec{\textbf{\emph{E}}}_{R}(t)$ and the macroscopic electric field $\vec{\textbf{\emph{E}}}(t)$. From formula (\ref{Polarization})
\begin{eqnarray}
\label{Reaction_reflection}
\vec{\textbf{\emph{E}}}_{R}(t) =  &-& \frac{\eta}{2}  \dot{\vec{\textbf{\emph{P}}}}(t)=\nonumber \\
&-& \frac{i \eta}{2} \omega \chi \epsilon_0  \vec{\textbf{\emph{E}}}(t)=- \frac{i k\chi}{2+ik\chi} \vec{\textbf{\emph{E}}}_{i}(t)
\end{eqnarray}
From formulas (\ref{Fresnel1}) and (\ref{Fresnel3}) it is clear that for an $s$ polarized incident wave
\begin{eqnarray}
\label{Reaction_reflection2}
\vec{\textbf{\emph{E}}}_{R}(t) =r_s \vec{\textbf{\emph{E}}}_{i}(t)
\end{eqnarray}
Formulas (\ref{Reaction_reflection2}) states that $\vec{\textbf{\emph{E}}}_{R}(t)$ is the macroscopic electric field irradiated by the 2D crystal. It gives rise to the reflected field \cite{Merano16} and in superposition with $\vec{\textbf{\emph{E}}}_{i}(t)$ to the transmitted field i.e to the total macroscopic electric field. In formulas
\begin{eqnarray}
\label{Reaction_Macro}
\vec{\textbf{\emph{E}}}(t) =\vec{\textbf{\emph{E}}}_{R}(t)+\vec{\textbf{\emph{E}}}_{i}(t)
\end{eqnarray}
This field is continuous and well-defined also at the crystal location.

\section{Conclusions}

An insulating 2D crystal has been modelled as a 2D Bravais lattice of dipoles. If an electromagnetic plane wave is incident on the crystal, each dipole will scatter it. The effect of the dipoles is obtained by superposition of the scattered wavelets. The scattered radiation will combine coherently to generate the macroscopic transmitted and reflected fields.

The total electric field for this system is the linear superposition of the incident electric field plus the retarded field generated by all the dipoles. The local electric field acting on a single dipole is the linear superposition of the incident electric field plus the retarded fields generated by all the other dipoles. The macroscopic electric field entering in the boundary conditions for the system is provided by the linear superposition of the incident electric field plus the radiation-reaction electric field. This last result is the main contribution of this article.

Wheeler and Feynman, in a very nice paper on the mechanism of the radiative reaction, conclude that radiation is a phenomenon as much of statistical mechanics as of pure electrodynamics \cite{Wheeler45}. Here a fundamental statistical aspect of the radiation-reaction electric field for a 2D crystal has been enlightened. A relation in between the microscopic and the macroscopic description has been found. Remarkably the advanced fields play an important role in this context.

The approach here derived is completely original, it is based on a well-defined and relativistically invariant prescription \cite{Dirac38, Wheeler45}, and it is well supported by experiments \cite{Merano16, Blake2007, Nair2008, Kravets2010}. Starting from a microscopic theory it connects the microscopic and the macroscopic electric fields in the crystal, fixing the bondary conditions for these last ones and in the end the Fresnel coefficients of the 2D crystal.

This paper provides a complete classical theory of radiating two-dimensional crystals. It furnishes for the first time the expression of the radiation-reaction field in these systems, making an essential step in understanding their optical response. This result is strictly valid for single-layer 2D crystals. Anyway the developed method \cite{Luca16} and the results here obtained pave the way to a complete treatment of the optical response of few-layer atomic crystals, hetero-structures \cite{Geim14} and metasurfaces \cite{Kuester03, Capasso14, Arbabi15, Shalaev13, Brongersma14}.

\section{ACKNOWLEDGMENTS}
I acknowledge Luca Dell'Anna for useful discussions.