\section{\label{sec:intro} Introduction }

At their inception in the 30's, Monte Carlo algorithms intended to describe the neutron transport in fissile media \cite{Spanier1969MonteProblems1,Lux1991MonteCalculations} in support of the first applications of nuclear power, might it be for the design of the first nuclear weapons or in support of nuclear energy production. Since then, the nuclear energy industry has used the so-called Monte Carlo criticality codes as reference numerical schemes to solve the linear Boltzmann equation. Nuclear safety demonstrations in particular rely on a Monte Carlo solving of this equation with few, if any, hypotheses, while the use of deterministic codes is preferred for an approximate but fast -and even online- solving of a simplified version of this equation: the two (energy) group critical diffusion equation. In both cases, the aim of such neutronics codes \footnote{neutron transport theory, also known as neutronics, aims at characterizing the properties of the neutron gas in fissile media} is to calculate the neutron spatial distribution in the context of mean-field hypotheses. Indeed, even at the startup of a nuclear reactor (i.e. when the neutron densities are small), such mean-field equations are believed to accurately describe the average behavior of the stochastic regime. Though, in the power regime of the reactor, stabilizing effects come into action. For example, local power excursions and potential sudden rise of local neutron densities are tempered by the Doppler broadening of the neutron cross-sections (giving the probabilities that a neutron at a given energy induce a particular reaction on a nuclei): the quantum resonances of the neutron-nuclei system are broadened by the thermal agitation of the nuclei from which results a modification of the neutron capture cross section. The study of the properties of neutron transport in fissile media subject to feedback mechanisms is called 'reactor physics'. Reactor physics numerical solvers hence couple the mean-field equations of neutronics to thermal or thermal-hydraulics solvers via the definition of intermediate quantities (such as the temperature scalar field) to extract macroscopic measurable quantities characterizing the neutron population (such as local neutron fluxes measured by fission rate chambers for example).

The fluctuations of the neutron population were also largely investigated in the past fifty years \cite{bell_glasstone,williams_book,paszit_pal}, while being assumed to be of marginal concern for nuclear safety, since they vanish as the core power increases. However, a few years ago, this paradigm was questioned \cite{Dumonteil2014ParticleSimulations,Zoia2014ClusteringGeometries,Sutton2017NeutronCalculations} by the study of a neutronics toy-model based on the celebrated branching Brownian motion \cite{harris_book}, that couples the diffusive random walk of neutrons \cite{Zoia2012BranchingStatistics}, mimicking their transport in the Brownian regime, to a Galton-Watson birth-death process, reproducing the induced fission reaction with variable number of outgoing neutrons. Indeed, the first moments of the master equation of this model revealed that the spatial correlations within the neutron population might in some cases diverge, invalidating the use of the mean-field equation. From a physical point of view, this clustering phenomenon (already characterized at the time in the theory of population ecology \cite{Young2001ReproductiveOrganisms, Houchmandzadeh2002ClusteringOrganisms, Houchmandzadeh2008NeutralCommunity}) causes the emergence of spatial patterns in the neutron gas. While always present and increasing with time in 1 and 2 dimensional systems, the spatial correlations however saturate in 3 dimensions. A dedicated experiment taking place at the Rensselaer Poytechnic Institute was designed and took place in 2018 to characterize such spatio-temporal correlations affecting the neutron distributions. This experiment revealed the existence of such correlations and triggered questions relative to their persistence in the power regime of reactor physics where feedbacks occur \cite{Dumonteil2021PatchyReactions}.\newline

In this paper, leveraging on the seminal approach of Janssen and De Dominicis \cite{janssen1976,dominicis1976}, we will build a simplified model of reactor physics to investigate the validity of mean field equations in the power regime of nuclear reactors. In that aim, a field theoretic formulation of neutron physics will be built in Section 2. While being 0 dimensional, this formulation will then be exploited in Section 3 so as to take into account the principal features of neutronics such as the calculation of fluctuations and correlations of the neutron field or the effects of an external neutron source. It will also be extended in order to add the effect of delayed neutrons in Section 4 and to take into account spatial phenomena in Section 5. Along the path, we will highlight the close relationship between  the response functional formalism of Janssen and De Dominics with the Doi and Peliti field theory \cite{doi_1976, peliti_1985} which lend itself to an exact microscopic interpretation, thus strengthening the use of both approaches in neutronics. Finally, in Section 6, we will add a simplified model of thermal feedback, showing that the neutron gas in a nuclear reactor operated at criticality in the power regime can be described as a phase transition belonging to the directed percolation universality class. The calculation of the critical exponents at the percolation threshold using the renormalization group will in return allow to specify the mean field equation, formulated in terms of fractional Laplacian. We will in particular show that in the directed percolation model, the critical diffusion equation shifts to a critical super-diffusion equation, whose stochastic generators are Lévy flights, ultimately dynamically modifying the medium properties in which neutrons are propagating. The analysis of the fundamental mode of the associated fractional Laplacian using an approximate formula will allow to quantify the differences between this formal approach and classical coupling schemes, thus setting an accuracy limit on all actual numerical solvers of reactor physics.

\section{ A Field theoretic formulation of neutronics}
\label{sec:FieldFormulation}
Branching processes are among the simplest models capable of characterizing the phenomenology of neutrons evolving in fissile materials \cite{paszit_pal,williams_book}. In their simplest expression, these models approximate the path of neutrons in matter as a stochastic process, in which they are subject to random extinction or reproduction events, occurring at constant average rates. They are also known to display a second order phase transition between an active and an inactive state \cite{harris_book}. This transition is at the core root of nuclear reactors operations and corresponds, in this context, to the appearance of self-sustained chain reactions: the 'percolation threshold' where neutron chains survive indefinitely is precisely the point at which nuclear reactors do operate.

Field theoretic method have long been used in statistical and condensed matter physics for their efficiency in the study of the intricate behavior of systems undergoing phase transitions. They indeed prove particularly adapted to unravel the universal properties and scaling behaviors in the vicinity of the phase transition's critical point. Along this line of thought, the present paper ought to demonstrate that field theory also provides an adapted framework for neutronics, allowing notably to unravel highly non-trivial behaviors of the neutron population close to \textit{criticality}.

The derivation of a field-theoretic version of neutronics followed in this paper is based on the seminal work of Janssen \cite{janssen1976} and De Dominicis \cite{dominicis1976}. Starting from the mean field equation of the problem, to which is supplemented a random noise source term (as is commonly done to study the stability of nuclear reactors \cite{bell_glasstone}), the neutron population evolving in a nuclear reactor can be decoupled into two components:
\begin{itemize}
    \item a smoothly varying neutron field $N(t)$. It corresponds to the instantaneous number of neutrons that can be measured in a reactor (proportional to the operating power for instance) and it is only sensitive to the global reactivity $\rho$ of the reactor;

    \item a random noise source $\eta(t)$, arising from various perturbations (such as the stochasticity of the induced fission reactions that have a random number of outgoing neutrons, or such as assembly vibrations) which are generally driven by unknown or at least unspecified phenomena.
\end{itemize}

In the simplest model of the branching process associated to neutrons evolving in an idealized, infinite and homogeneous reactor, such a decoupling transforms the mean field equation into a Stochastic Differential Equation (SDE), whose form is given by
\begin{equation}
\label{LangevinEquation}
    \cfrac{d N(t)}{dt} = \rho N(t) + \eta(t).
\end{equation}

Any macroscopic observable quantity $\mathcal{O}[N]$ is insensitive to the rapidly varying and low lying noise term. As such, physical observables must be averaged over all possible configurations of the noise. This is formally achieved by performing a functional integration over the noise probability functional $P[\eta(t)]$,
\begin{equation}
\label{observable}
    \langle \mathcal{O}[N] \rangle \, \propto \, \int \mathcal{D}\eta \, \mathcal{O}[N] \, P[\eta].
\end{equation}

The development of the field-theoretic model starts by noting that the $N$ appearing under the integral sign of Eq.~\ref{observable} is constrained to be a solution of Eq.~\ref{LangevinEquation}. The constraint (valid for all $t$),
\begin{equation}
    \label{constraint_N}
    \mathcal{C}[N] \, = \, \left\{\frac{dN(t)}{dt} \, - \, \rho \, N(t) \,  -  \, \eta(t)\right\} = 0,
\end{equation}
can be enforced with the help of a functional version of the resolution of the identity \cite{cardy_tauber},
\begin{equation}
    1 = \int \mathcal{D}N \, \prod_t  \delta \left( \mathcal{C}[N] \right).
\end{equation}
Performing a (functional) Fourier transform, one can write
\begin{eqnarray}
\label{constrained_field}
    1 & = & \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \, e^{- \int dt \ensuremath{\widetilde{N}}(t) \, \mathcal{C}[N]} \nonumber\\
    & = & \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \, e^{- \int dt \ensuremath{\widetilde{N}}(t) \, \left( \frac{dN(t)}{dt} \, - \, \rho \, N(t) \,  -  \, \eta(t)\right)},
\end{eqnarray}
where the purely imaginary auxiliary field $\ensuremath{\widetilde{N}}(t)$, known as the Martin-Siggia-Rose response field \cite{martin_siggia_rose}, has been introduced. Inserting this Fourier transformed identity back into Eq.~\ref{observable}, one obtains
\begin{eqnarray}
\label{observable2}
    \langle \mathcal{O}[N] \rangle & \propto & \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \left\{ e^{- \int \ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt}  -  \rho   \right)N(t) \, dt} \right. \nonumber \\
    & \times &\mathcal{O}[N]
   \left.  \times  \int \mathcal{D}\eta \, e^{\int \ensuremath{\widetilde{N}}(t) \eta(t) dt} P[\eta] \right\}.
\end{eqnarray}
To go further, one needs to integrate over $\eta$ and thus specify the probability functional to which it is associated. In the following, it will be supposed that at any given moment in time, the noise follows a Gaussian probability distribution
\begin{equation}
\label{noise_proba}
    P[\eta] \propto e^{ -\frac{1}{4}\int \eta(t) \Gamma^{-1} \eta(t) \, dt }
\end{equation}
with zero mean $\langle \eta(t) \rangle = 0$ and a covariance given by,
\begin{equation}
    \langle \eta(t) \eta(t') \rangle = 2 \Gamma(N) \delta(t - t').
\end{equation}
Note that an explicit dependence of the variance $\Gamma(N)$ on the global power level (i.e. $N(t)$) is to be expected on physical ground : the expected noise in a reactor indeed depends on the operating power level. The integral over the noise can now be straightforwardly evaluated as
\begin{equation}
\int \mathcal{D}\eta \, e^{\int \ensuremath{\widetilde{N}}(t) \eta(t) dt} P[\eta] \propto e^{ \Gamma \ensuremath{\widetilde{N}}^2}.
\end{equation}
Inserting this result back into Eq.~\ref{observable2}, one finally ar\-rives at an expression that solely depends on the variance of the noise probability distribution
\begin{equation}
\label{observable_final}
    \langle \mathcal{O}[N] \rangle \, = \mathcal{N}  \,  \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \, \mathcal{O}[N] \,
    e^{- \mathcal{S}[N,\ensuremath{\widetilde{N}}]},
\end{equation}
where both a normalisation factor $\mathcal{N}$ and the so-called \textit{response functional} have been introduced, this last functional being defined through
\begin{equation}
\label{response_functional}
     \mathcal{S}[N,\ensuremath{\widetilde{N}}] = \int \left\{\ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt} - \rho \right)N(t) - \ensuremath{\widetilde{N}}(t) \Gamma \ensuremath{\widetilde{N}}(t)\, \right\} dt.
\end{equation}

It often reveals convenient to include in Eq.~\ref{constraint_N} an initial condition of the type $N(t_0) = N_0$. The initial condition is then included directly into the response functional and it will be taken into account naturally in the further calculation of observables. One can show that the inclusion of an initial condition in the response functional gives
\begin{eqnarray}
\label{response_functional_initial_cond}
     \mathcal{S}[N,\ensuremath{\widetilde{N}}]  & \, & = \int \left\{  \ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt} - \rho \right)N(t)          \right.\nonumber\\
      && \left. -  \ensuremath{\widetilde{N}}(t) \, \delta(t-t_0) \, N_0   -   \ensuremath{\widetilde{N}}(t) \Gamma \ensuremath{\widetilde{N}}(t)\, \right\} dt.
\end{eqnarray}
It can then be demonstrated that the net effect of this supplementary term appearing in the response functional is to multiply each quantity evaluated at the initial time by the factor $N_0$. In the remainder of the paper, only the latter requirement will be retained and the terms corresponding to the initial conditions will be systematically omitted in all of the response functionals that will be considered. \newline

Equations \ref{observable_final}--\ref{response_functional_initial_cond} form the field-theoretic version of the branching process. It presents itself under the form of a (euclidean) path integral. The further calculation of observable quantities can now be addressed thanks to the powerful machinery developed in  Quantum Field Theory (QFT) and statistical mechanics.

\section{Calculation of observables}
\label{sec:observables}

It can generally safely be assumed that in a nuclear reactor, the random power fluctuations driven by random noise are much weaker than the actual operating power level, i.e. the noise term is \textit{small}. This suggests a strategy to extract sensible results from the response functional formalism. One can expand the function $\Gamma(N)$ into the Taylor series
\begin{equation}
\label{NoisePowerSeries}
    \Gamma(N) = \lambda_0 + \lambda_1 N + \lambda_2 N^2 \, + ...,
\end{equation}
with small $\lambda_i$'s. Under this form, the noise term appearing in the response functional can be treated perturbatively, just as one does with an interaction potential in quantum mechanics. This perturbative treatment and the calculation of observables of the type given in Eq.~\ref{observable_final} is most conveniently represented under the form of Feynman diagrams \cite{greiner,zinn}

The observables derived from the response functional will all depend on the parameters $\lambda_i$, which play here the same role as the interaction coupling constants in conventional QFT. The precise meaning of these parameters and in particular their microscopic origin is, however, left unspecified. In this respect, the response functional formalism can be viewed as a purely phenomenological approach: the precise interpretation of the coupling parameters is to be sought either a posteriori or from an external knowledge of the problem.

One might thus rightly be worried that the calculation of observables relies upon an infinite power series expansion, whose structure and convergence depends on unknown parameters. In fact, the inclusion of powers higher than one in Eq.~\ref{NoisePowerSeries} leads in some cases to diverging expressions, rendering the perturbative expansion potentially inoperant (an example of which will be encountered in the last section). This type of divergence is a common problem in QFT and can generally be cured if one insists that all observables quantities truly describe the macro-scale behaviour of the problem, i.e. that all small scale fluctuations are truly being 'integrated over'. In a nutshell, the observables must be derived at a fixed point of a spatio-temporal scaling transformation of the system. By simple dimensional arguments \cite{janssen2004}, one finds that high order terms in the Taylor expansion of the function $\Gamma$ dies off when such a scaling operation is performed. In the Renormalization Group (RG) jargon, they represent \textit{irrelevant} parameters which formally vanish in the vicinity of the RG fixed point. These higher order terms can thus be directly discarded at this stage and $\Gamma$ may thus be simply expressed as

\begin{equation}
    \Gamma(N) =  \lambda_0 +  \lambda_1 N.
\end{equation}

\subsection{ Moments of the neutron population }

Because the noise term has zero mean, it cannot impact the mean number of neutrons $\langle N(t) \rangle = N(t)$. As long as mean values are concerned, one can thus reasonably forget about the noise term in the response functional. Using Eq.~\ref{observable_final}, the mean number of neutrons can be written as,
\begin{equation}
\label{mean_N}
    \langle N \rangle \, = \mathcal{N}  \,  \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \, N \,
    e^{-  \int \left\{\ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt} - \rho\right)N(t) \right\} dt } .
\end{equation}
This expression can be simply evaluated with the help of the functional version of Wick's theorem \cite{zinn}. Taking into account for the initial condition $N(0) = N_0$, one has,
\begin{equation}
    \langle N(t) \rangle = N_0 \times \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}
      \vertex  (a) {\(t\)};
      \vertex [right=1.5cm of a](o1) {\(t_0\)};
    \diagram* {
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture} =
    N_0 \, \mathcal{G}(t-t_0).
\end{equation}
$\mathcal{G}(t-t_0) = \langle N(t) \ensuremath{\widetilde{N}}(t_0) \rangle$ is the retarded Green function, solution of the equation
\begin{equation}
    \left(\frac{d}{dt} - \rho \right) \mathcal{G}(t-t') = \delta(t - t'),
\end{equation}
so that,
\begin{equation}
     \mathcal{G}(t-t') = e^{ \rho (t-t') }\, \theta(t - t'),
\end{equation}
with $\theta(t-t')$ being the Heaviside function, equals to zero whenever $t < t'$ and one otherwise. The mean number of neutrons then reads
\begin{equation}
    \langle N(t) \rangle = N_0 \, e^{ \rho (t-t_0) }\, \theta(t - t_0).
\end{equation}

The second moment $\langle N^2(t) \rangle$ is less trivial, because it depends on the noise variance  $\Gamma(N)$. Forgetting temporarily the $\lambda_0$ term in the series expansion of $\Gamma$, and resorting again to Wick's theorem, one can symbolically write
\begin{eqnarray}
\label{observable_second_moment}
    \langle N^2(t) \rangle \, & = &  2 \lambda_1 \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of a](o1);
    \diagram* {
    (f1) -- (a) -- (f2),
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + 2 \; \;
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above= -0.3cm of a](c) ;
      \vertex [above= 0.4cm of a](b) ;
      \vertex [right=of c](o1) ;
      \vertex [right=of b](o2) ;
    \diagram* {
    (c) -- (o1),
    (b) -- (o2),
    };
    \end{feynman}
  \end{tikzpicture} \nonumber\\ & & \nonumber\\
  & = &  \langle N^2(t) \rangle_c + \langle N(t) \rangle^2,
\end{eqnarray}
where the notation $\langle \rangle_c$ is used to represent the first, so-called \textit{connected}, diagram. The variance can then be simply evaluated as ($\rho \neq 0$)
\begin{eqnarray}
  \langle N^2 \rangle - \langle N \rangle^2 & = & \langle N^2 \rangle_c  =  2 \lambda_1 N_0 \int_{t_0}^t d\tau \mathcal{G}^2(t - \tau)  \mathcal{G}(\tau - t_0) \nonumber\\
  & = &   \cfrac{ 2 \lambda_1 N_0}{\rho}  \,  e^{\rho  (t-t_0)} \left[ e^{\rho  (t-t_0)} -1  \right].
\end{eqnarray}
This expression exactly matches the variance one would derive from a detailed microscopic balance of the neutron population (in an infinite, homogeneous medium) \cite{paszit_pal}, provided that $\lambda_1$ be identified with the rescaled variance $(\bar{\nu}^{(2)} - \rho)/2$ of the offspring distribution of the underlying branching stochastic process. This result is a good illustration of the general outreach of the response functional formalism.  Relying only on the SDE built upon the 'coarse-grained' description of the problem, the formalism is able to take advantage of the phenomenological and microscopic random source term to extract residual information about the structure of the fluctuations around the mean field behaviour of the system.

A strong limitation of this approach is the assumption that the noise term follows a Gaussian profile. Indeed, the calculation of higher moments $\langle N^k \rangle$ for the neutron population would formally require information about the same moments of the offspring distribution. In the response functional formalism, there is no possibility to include this type of information in the calculations. This limitation does not however applies in the present case, at least in the vicinity of criticality. It has indeed been demonstrated in \cite{garcia_millan} that the microscopic branching process under study displays a \textit{universal} behaviour in the vicinity of the critical point $\rho = 0$: it then solely depends on the reactivity $\rho$ and the second moment of its offspring distribution. Thus the Gaussian distribution imposed on the noise source does not limit the extraction of higher moments of the neutron distribution.

The approach proposed in \cite{garcia_millan} also relies on a field-theoretic formulation of the stochastic branching process, pioneered in the work of Doi \cite{doi_1976} and Peliti \cite{peliti_1985}. Starting from a formal identification of the ladder operators $a^\dagger$ and $a$ with the multiplication or derivation by $z$ of the generating function $\psi(z,t) = \sum_n P(n,t) z^n$, one can give a structure to the stochastic process that closely resemble the structure of Fock spaces used in quantum mechanics \cite{baez_2019}. A path integral formulation of the process is then built upon the coherent states of the annihilation operator $a$ \cite{tauber_2005}.

Both the Doi-Peliti and Janssen-De Dominicis approaches allow for the derivation of a field theory, from two different starting points. In this respect, they can be viewed as complementary:
\begin{itemize}
    \item the Doi-Peliti approach relies on the explicit modelling of the underlying microscopic stochastic process. Every term in the theory is given a sound microscopical interpretation from the start. The drawback is the difficulty to include terms that cannot be described as stochastic processes (because their microscopic origin is unclear for instance);

    \item the Janssen-De Dominicis approach is much more phenomenological. Its starting point is the mean field equation associated to the problem. This enables more flexibility in the modelling of phenomena that could not be given a microscopic origin (e.g. the modelling of counter-reaction mechanisms in neutronics for instance). Relying on Gaussian noise distributions and unspecified parameters, this approach is nevertheless more limited in scope when a microscopic model is available.
\end{itemize}

Thanks to the universal behaviour of branching processes close to their critical point, both approaches coincide in the present case: they both lead to the same path integral formulation and the same diagrammatic expansion. As such, the extraction of higher moments for the neutron population in the response functional formalism exactly matches the derivation presented in \cite{garcia_millan}.

\subsection{Addition of an external neutron source}

At low power, external neutron sources, such as the spontaneous fission on some heavy nuclei or $(\alpha,n)$ reactions might have a substantial impact on the moments of the neutron population \cite{williams_book,Dumonteil2021PatchyReactions}. The response functional formalism lends itself easily to describe this phenomenon. Considering an external neutron source, emitting an average of $S$ neutrons per second, the SDE describing the problem now reads
\begin{equation}
    \label{sde_source}
        \cfrac{d N(t)}{dt} = \rho N(t) + S + \eta(t),
\end{equation}
as well as the associated response functional,
\begin{eqnarray}
\label{action_source}
     \mathcal{S}[N,\ensuremath{\widetilde{N}}]  & = & \int \left\{\ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt} -  \rho \right)N(t) - \ensuremath{\widetilde{N}} S   \right.\nonumber\\
     & &  - \left. \ensuremath{\widetilde{N}}(t) \Gamma \ensuremath{\widetilde{N}}(t)\, \right\} dt.
\end{eqnarray}
Considering that the source term also comes with an intrinsic noise, we now can set $\lambda_0 \neq 0$ in the Taylor expansion of $\Gamma(N) = \lambda_0 + \lambda_1 N$.

The moments of the neutron population must now take into account the external source. In the perturbative calculation, this is done in the form of new Feynman diagrams appearing in the calculations.
For the mean number of neutrons at time $t$, one has (setting $t_0 = 0$ in the following),
\begin{eqnarray}
     \langle N(t)\rangle &=  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a](o1);
    \diagram* {
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  +  S \,
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,crossed dot] (o1) { };
    \diagram* {
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture} \nonumber \\
  ~ \nonumber\\
  & = N_0 e^{ \rho t } + S \int_0^t e^{ \rho (t-t') } dt' \nonumber\\
  ~ \nonumber\\
  & = N_0 e^{ \rho t } + \frac{S}{\rho} \left(  e^{ \rho t } - 1 \right) ,
\end{eqnarray}
in agreement with the literature \cite{williams_book}. The second propagator that now appears in the diagram directly account for the creation of neutrons produced at a constant rate $S$ by the source.

As for the variance of the neutron population, one has,
\begin{equation}
     \langle N^2(t)\rangle_c = 2 \lambda_1
     \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of a](o1);
    \diagram* {
    (f1) -- (a) -- (f2),
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
     +  2 \lambda_1 S
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of a, crossed dot](o1) { };
    \diagram* {
    (f1) -- (a) -- (f2),
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + 2 \lambda_0
 \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex [empty dot] (a) { };
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
    \diagram* {
    (f1) -- (a) -- (f2),
    };
    \end{feynman}
  \end{tikzpicture}
\end{equation}
The first diagram has already been calculated in the previous section. The following two diagrams take into account respectively:
\begin{itemize}
    \item the variance of the neutrons' branching process after their creation by the source $S$;
    \item the intrinsic variance $\lambda_0$ of the neutron source.
\end{itemize}
The calculation of the two unknown diagrams leads to the final expression ($\rho \neq 0$):
\begin{eqnarray}
 \langle N^2(t)\rangle_c & = & \frac{ 2 \lambda_1 N_0 }{\rho}  e^{\rho t } \left( e^{\rho t } - 1 \right)  +
 \frac{2 \lambda_1}{\rho} \left( e^{\rho t } - 1 \right)^2 \nonumber\\
 & + & \frac{\lambda_0}{\rho}   \left( e^{  2 \rho t } - 1  \right).
\end{eqnarray}
 The calculation of moments of higher orders now poses a real problem. Although the universal behavior of the neutrons' branching process could, here also, be invoked, one is now limited by the noisy component of the external source term. Only if the source possesses a Gaussian distribution can the calculation be taken further and the moments $\langle N^k\rangle_c$ with $k>2$ be exactly calculated.
\section{Two species model}

Another key aspect in the modelling of the temporal behavior of the neutron population in a nuclear reactor is the role played by the precursors of delayed neutrons $M(t)$. In rare cases indeed, an induced fission does not emit all the outgoing neutrons instantaneously as the fission fragments can also emit neutrons on longer time scale controlled by their electroweak decay, ranging from a second to a few minutes, orders of magnitudes longer than prompt neutrons and occurring in the present model at a rate $\lambda_D$. Thus, while emitting a small fraction $\beta$of the overall number of neutrons, precursors control the kinetics of nuclear reactor which could not be operated on human time scales otherwise.

Including the coupling between the neutron and precursor populations, the point-like model transforms to (with $r = \rho - \beta)$
\begin{equation}
\label{pointKinetics}
\begin{cases}
        \cfrac{d N(t)}{dt} & = r N(t) + \lambda_D M(t) + \eta_1(t), \\
         \\
    \cfrac{d M(t)}{dt} & = \beta N(t) - \lambda_D M(t) + \eta_2(t) ,
\end{cases}
\end{equation}
where external sources have been neglected, for simpli\-city. As before, the two random source terms are supposed to follow Gaussian probability distributions with zero mean and a covariance given by the relation ($a,b = 1,2$),
\begin{equation}
    \langle \eta_a(t) \eta_b(t') \rangle = 2  \, \Gamma_{ab}(N, M) \,  \delta(t - t'),
\end{equation}
The derivation of the response functional associated to this coupled system of equations still follows the same recipe as the one described in section \ref{sec:FieldFormulation}. Four independent fields $\ensuremath{\widetilde{N}}, \ensuremath{\widetilde{M}}, M$ and $N$ are now introduced into the theory and the response functional can be written $\mathcal{S} = \mathcal{S}_{\text{neut.}} + \mathcal{S}_{\text{prec.}} + \mathcal{S}_{\Gamma}$, with,
\begin{equation}
\label{action_neutron_2species}
   \mathcal{S}_{\text{neut.}}  =  \int \left\{\ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt}-r\right)N(t) -  \lambda_D \ensuremath{\widetilde{N}}(t) M(t)  \right\} dt,
\end{equation}

\begin{equation}
\label{action_prec_2species}
     \mathcal{S}_{\text{prec.}} =  \int \left\{ \,\ensuremath{\widetilde{M}}(t)  \left( \frac{d}{dt} + \lambda_D \right)M(t) - \beta \ensuremath{\widetilde{M}}(t) N(t) \right\} dt,
\end{equation}
and, in its most general form,
\begin{eqnarray}
       & & \mathcal{S}_{\Gamma} =  \int \Gamma_{11} \ensuremath{\widetilde{N}}^2 +  \Gamma_{22} \ensuremath{\widetilde{M}}^2 + 2 \Gamma_{12} \ensuremath{\widetilde{N}} \ensuremath{\widetilde{M}}  \nonumber\\
 & =&  \int   \left\{ \lambda_1 \ensuremath{\widetilde{N}}^2 N
  +  \lambda_2 \ensuremath{\widetilde{N}}^2 M + \lambda_3 \ensuremath{\widetilde{M}}^2 N \right.\nonumber\\
   &\hspace{0.3 cm}&  + \left. \lambda_4 \ensuremath{\widetilde{M}}^2 M  + 2 \lambda_5 \ensuremath{\widetilde{N}} \ensuremath{\widetilde{M}}  N + 2 \lambda_6 \ensuremath{\widetilde{N}} \ensuremath{\widetilde{M}} M\right\} dt.
\end{eqnarray}
Two Green functions are to be associated to the free parts of the neutron and precursor functionals,
\begin{equation}
   \mathcal{G}_N(t - t') =  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}
      \vertex  (a) {\(t\)};
      \vertex [right=1.5cm of a](o1) {\(t'\)};
    \diagram* {
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture} =
    e^{ r(t-t') } \theta(t - t')
\end{equation}
and
\begin{equation}
   \mathcal{G}_M(t - t') =  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}
      \vertex  (a) {\(t\)};
      \vertex [right=1.5cm of a](o1) {\(t'\)};
    \diagram* {
    (a) -- [scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture} =
    e^{ -\lambda_D (t-t') } \theta(t - t').
\end{equation}

Because of the couplings between the neutron and precursor populations appearing in Eq.~\ref{action_neutron_2species} and Eq.~\ref{action_prec_2species}, observables must however be ultimately  constructed from a generalized set of propagators, which solve the full coupled set of equations Eq.~\ref{pointKinetics}. These generalized propagators may be symbolically written,

\begin{widetext}
\begin{eqnarray}
 &  & \Sigma_{NN}(t,t') = \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
    (a) --  (b) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture} =
     \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a](o1);
    \diagram* {
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.5 cm of c](o1);
    \diagram* {
    (a) -- (b) ,
    (b)--[scalar] (c),
    (c) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \left( \beta \, \lambda_D  \right)^2
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.7 cm of c,dot](d) {};
      \vertex [right=0.7 cm of d,dot](e) {};
      \vertex [right=0.5 cm of e](o1);
    \diagram* {
    (a) -- (b) ,
    (b)--[scalar] (c),
    (c) -- (d) ,
    (d)--[scalar] (e),
    (e) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + ... \, , \, \label{sig_nn} \\ & & \nonumber \\
    & & \Sigma_{MM}(t,t')  =
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
    (a) -- [scalar] (b) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture} =
     \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a](o1);
    \diagram* {
    (a) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.5 cm of c](o1);
    \diagram* {
    (a) --[scalar] (b) ,
    (b)-- (c),
    (c) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \left( \beta \, \lambda_D  \right)^2
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.7 cm of c,dot](d) {};
      \vertex [right=0.7 cm of d,dot](e) {};
      \vertex [right=0.5 cm of e](o1);
    \diagram* {
    (a) --[scalar] (b) ,
    (b)-- (c),
    (c) --[scalar] (d) ,
    (d)-- (e),
    (e) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + ... \, , \, \\ & & \nonumber \\
& & \Sigma_{NM}(t,t')  =  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
    (a) -- [scalar] (b) --(o1)
    };
    \end{feynman}
  \end{tikzpicture} = \beta \left\{ \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b](o1);
    \diagram* {
    (a) --[scalar] (b),
    (b) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \, \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b,dot](c) {};
      \vertex [right=of c,dot](d) {};
      \vertex [right=of d](o1);
    \diagram* {
    (a) --[scalar] (b) ,
    (b)-- (c),
    (c)-- [scalar] (d),
    (d) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + ...\right\} , \\ & & \nonumber\\
 & &  \Sigma_{MN}(t,t')  =  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.8 cm of b](o1);
    \diagram* {
    (a) --  (b) -- [scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture} = \lambda_D \left\{ \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b](o1);
    \diagram* {
    (a) -- (b),
    (b) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \, \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b,dot](c) {};
      \vertex [right=of c,dot](d) {};
      \vertex [right=of d](o1);
    \diagram* {
    (a) -- (b) ,
    (b)-- [scalar](c),
    (c)-- [] (d),
    (d) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
 + ...\right\}.
\end{eqnarray}
\end{widetext}
The calculation of each one of the $\Sigma_{IJ}$ involves an infinite series of nested integrals. Fortunately, these infinite series can be summed over and it is possible to obtain analytical results in all four cases (see annex \ref{app_a}).
A much more convenient but indirect way to extract an analytical expression out of these propagators is to realize that they all directly correspond to the mean value of either the neutron or precursor populations for different initial conditions. Indeed, setting $( N(t_0), M(t_0)) =(N_0, 0 )$, one can show that,
    \begin{equation}
        \langle N(t) \rangle_{(N_0, 0 )} = N_0 \times \Sigma_{NN}(t,t_0),
    \end{equation}
    \begin{equation}
        \langle M(t) \rangle_{(N_0, 0 )} = N_0 \times \Sigma_{NM}(t,t_0),
    \end{equation}
whereas, for $( N(t_0), M(t_0)) =(0, M_0 )$,
    \begin{equation}
        \langle N(t) \rangle_{(0, M_0 )} = M_0 \times \Sigma_{MN}(t,t_0),
    \end{equation}
    \begin{equation}
        \langle M(t) \rangle_{(0, M_0 )} = M_0 \times \Sigma_{MM}(t,t_0),
    \end{equation}
so that each of these propagators is the solution of the system of equations Eq.~\ref{pointKinetics} for a specific set of initial conditions (and with the stochastic terms set to zero). This system of equations can be straightforwardly integrated and solutions extracted. The explicit form of each generalized propagator is provided in annex \ref{app_a}.

Equipped with an analytical expression for the propagators of the theory, the variance of the process may then be evaluated. Imposing the initial condition $( N(0), M(0)) =(N_0, 0 )$  , the variance of the neutron and precursor populations may be diagrammatically written,
\begin{widetext}
\begin{eqnarray}
\label{var_N}
     \langle N^2(t) \rangle_c & = & 2 N_0 \lambda_1 \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- (f1)  -- (a) --  (f2) -- (ff2),
    (a) -- (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + 2 N_0 \lambda_2
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- (f1) -- (a) -- (f2) -- (ff2),
    (a) -- [scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + 2 N_0 \lambda_3
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- (f1) --[scalar] (a) --[scalar]  (f2) -- (ff2),
    (a) -- (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture} \nonumber\\
   & & + 2 N_0 \lambda_4
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- (f1) --[scalar] (a) --[scalar]  (f2) -- (ff2),
    (a) --[scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + 2 N_0 \lambda_5
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- (f1) --[scalar] (a) --  (f2) -- (ff2),
    (a) -- (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
    + 2 N_0 \lambda_6
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- (f1) --[scalar] (a) --  (f2) -- (ff2),
    (a) -- [scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture} ,
\end{eqnarray}
\begin{eqnarray}
\label{var_M}
     \langle M^2(t) \rangle_c &=& 2 N_0 \lambda_1 \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) --[scalar] (f1)  -- (a) --  (f2) --[scalar] (ff2),
    (a) -- (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + 2 N_0 \lambda_2
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) --[scalar] (f1) -- (a) -- (f2) --[scalar] (ff2),
    (a) -- [scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + 2 N_0 \lambda_3
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) --[scalar] (f1) --[scalar] (a) --[scalar]  (f2) --[scalar] (ff2),
    (a) -- (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture} \nonumber \\
   & & + 2 N_0 \lambda_4
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) --[scalar] (f1) --[scalar] (a) --[scalar]  (f2) --[scalar] (ff2),
    (a) --[scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
    + 2 N_0 \lambda_5
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) --[scalar] (f1) --[scalar] (a) --  (f2) --[scalar](ff2),
    (a) -- (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
    + 2 N_0 \lambda_6
   \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) --[scalar] (f1) --[scalar] (a) --  (f2) --[scalar] (ff2),
    (a) -- [scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture} .
\end{eqnarray}
\end{widetext}

Just as it was the case in the last section, both expressions depend on the  set of parameters $\lambda_i$, which are not specified by the theory.
With a total of 6 unknown parameters for only two observables, it may seem at first sight difficult to give a precise meaning to each term. However, the structure of the diagrammatic expansion and the results of the one species model derived earlier allows to do so, at least partially.
In light of the results obtained in the previous section, it might indeed  be inferred that the diagrams associated to $\lambda_1$ appearing in both expressions represent the variance associated to the neutron’s offspring distribution. There is furthermore no reason to expect that the value for $\lambda_1$ obtained in the case of the one species model should be any different in the present situation,  except for the replacement of the reactivity $\rho$  with the prompt reactivity $r = \rho - \beta$.

On the same model, the diagrams associated with $\lambda_3$ are most certainly to be  identified with the variance of the precursor emitted by fission.

The diagrams associated to $\lambda_2$ and $\lambda_4$ share the same output structure as for $\lambda_1$ and $\lambda_3$ but possess a precursor at the origin of the branching term. This suggest that these diagrams encode for a source of variance intrinsic to the precursor population. The only process responsible for such kind of fluctuation in the precursor population is their decay process which, being a Poisson process, has a variance $\lambda_D$.\newline

At this point, one can thus postulate,
\begin{equation}
    \begin{cases}
            \lambda_1 = (\bar{\nu}_n^{(2)} - \rho + \beta)/2, \\
            \lambda_3 = (\bar{\nu}_m^{(2)} + \beta)/2, \\
            \lambda_2 = \lambda_4 = \lambda_D/2,
    \end{cases}
\end{equation}
so that only two unknown parameters, $\lambda_5$ and $\lambda_6$, remain to be determined for the two variances. Comparing equations Eq.~\ref{var_N}-\ref{var_M} with a microscopic balance of the neutron and precursor populations \cite{paszit_pal} for different values of $\rho$, $\beta$ and $\lambda_D$, one gets (the offspring distributions of neutrons and precursors were supposed independent),
\begin{equation}
    \begin{cases}
            \lambda_5 =  \overline{\nu_m ( \nu_n -1)} = \bar{\nu}_n^{(2)} \frac{\beta}{1-\beta}, \\
            \lambda_6 = - \lambda_D.
    \end{cases}
\end{equation}
A comparison of the diagrammatic expansion with the results of the literature is presented in Figure \ref{Var_NM}. The temporal evolution of the variance of both the neutron and precursor populations are displayed, in all three regimes: sub-critical ($\rho<0)$, critical $(\rho=0)$ and super-critical $(\rho>0)$. An excellent agreement is found in all three cases.

\begin{figure}[!ht]
\centering
\hspace{-1 cm}
\includegraphics[scale=0.75]{Figures/VNM.pdf}
\caption{Comparison of the variance of the neutron and precursor populations extracted from the solution of the stochastic process (plain lines) and from the calculation of Feynman diagrams (dotted lines), for $\rho =$ -0.2 (red), 0 (green), 0.2 (blue). The calculations are made assuming the following set of parameters: $\beta = 0.8$, $\lambda_D = 0.03$, $N_0 = 1000$, $\bar{\nu}_n^{(2)} = 2$, $\bar{\nu}_m^{(2)} = 1.2$. The two offspring distributions were supposed independent so that $ \overline{\nu_m(\nu_n - 1)} = \frac{\beta}{1-\beta} \bar{\nu}_n^{(2)} $. }
\label{Var_NM}
\end{figure}

\section{Stochastic diffusion equation}

Now equipped with the main phenomena of 0 dimensional neutronics, our field theoretic model can transition to proper neutron transport theory. A convenient and easy way to do so, yet realistic \footnote{the diffusion equation accurately reproduces the transport equation for isotropic media and far from the boundaries, see \cite{Duderstadt1979TransportTheory} for more details} (it is used in the operational context of reactor core monitoring), is to lean on the Brownian motion and its associated diffusion equation. The diffusion equation, which can also formally be built as a coarse grained description of a microscopic, stochastic process. When supplemented by an external random noise term, the (one energy group) equation can be written as
\begin{equation}
    \label{stoch_diffusion}
    \frac{\partial }{ \partial t} N(\vec{x},t) - D ( \Delta + r)  N(\vec{x},t) + \eta(\vec{x},t) = 0,
\end{equation}
where the rescaled parameter $r = \rho/D$ as been defined for convenience and where now $N(\vec{x},t)$ should be interpreted as the neutron density field. The transcription of this $d+1$ diffusion equation into a response functional follows, here also, exactly the same line as in section \ref{sec:FieldFormulation}. One thus have (in $d$ spatial dimensions),
\begin{equation}
    \mathcal{S} = \int d^d \vec{x} \int dt \ensuremath{\widetilde{N}}(\vec{x},t) \left( \partial_t - D ( \Delta + r)  \right) N(\vec{x},t) - \ensuremath{\widetilde{N}} \Gamma \ensuremath{\widetilde{N}}.
\end{equation}
The propagator of the theory is most conveniently written in the Fourier transformed spatial domain ($\vec{p}$,$t$),
 \begin{equation}
    \mathcal{G}(\vec{p},t) = e^{ -D(\vec{p}^{\,2} - r)t }.
\end{equation}
The structure of the interaction term of the response functional is not affected by the dimensionality of the problem. As a consequence, the diagrammatic expansion of the theory is the same as the one obtained in section \ref{sec:observables}. It can for instance be immediately concluded that provided with the initial condition $ N(\vec{x},t) = N_0 \, \delta(t) \delta^{(d)}(\vec{x})$, one finds that
\begin{equation}
    \langle N(\vec{x},t) \rangle = N_0  \mathcal{G}(\vec{x},t) = N_0 \, \sqrt{\frac{\pi}{4Dt}} \,  e^{ r D t - \frac{|\vec{x}|^2}{4 D t} }.
\end{equation}
A substantial amount of work has been done in the past, revolving around the possibility that neutrons evolving in fissile matter might gather into clusters and that the neutron density map should consequently present patchiness. This neutron clustering phenomenon results from the very nature of the underlying branching process they are subjected to \cite{Dumonteil2014ParticleSimulations, Zoia2014ClusteringGeometries,  DeMulatier2015TheRevisited,Agarwal2001, Sutton2017NeutronCalculations,Dumonteil2017ClusteringMedia,Dumonteil2021PatchyReactions}, and more specifically from the asymmetry between their uniform capture in the medium and their production only at spatial sites correlated to their parents \. This phenomenon can be quantified through the equal-time spatial correlation function $\langle N(\vec{x},t) N(\vec{y},t) \rangle$ which, in the response functional formal\-ism can be written as
\begin{widetext}
\begin{align}
\label{spatial_corr_func}
    & \langle N(\vec{x},t) N(\vec{y},t) \rangle_c  = 2 \lambda_1
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of a](o1);
    \diagram* {
     (f1)  -- [anti fermion] (a) -- [fermion] (f2),
    (a) -- [anti fermion] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
    = 2 \lambda_1 \int d^d\vec{w}  \int_0^t dt_1  \mathcal{G}(\vec{w},t_1) \mathcal{G}(\vec{x}-\vec{w},t-t_1) \mathcal{G}(\vec{y}-\vec{w},t'-t_1) \nonumber\\

    & = 2 \lambda_1 \int_0^t dt_1 \iiint  \frac{d^d\vec{p_1}}{(2\pi)^d} \frac{d^d\vec{p_2}}{(2\pi)^d} \frac{d^d\vec{p_3}}{(2\pi)^d} \mathcal{G}(\vec{p_1},t_1) \mathcal{G}(\vec{p_2},t-t_1) \mathcal{G}(\vec{p_3},t-t_1) e^{-i \vec{p_2} \cdot \vec{x}} e^{-i \vec{p_3} \cdot \vec{y}} \int d^d\vec{w} e^{-i \vec{w} \cdot (\vec{p_1} - \vec{p_2} - \vec{p_3})} \nonumber\\
    & = 2 \lambda_1 \int_0^t dt_1 \iiint  \frac{d^d\vec{p_1}}{(2\pi)^d} \frac{d^d\vec{p_2}}{(2\pi)^d} \frac{d^d\vec{p_3}}{(2\pi)^d} \mathcal{G}(\vec{p_1},t_1) \mathcal{G}(\vec{p_2},t-t_1) \mathcal{G}(\vec{p_3},t-t_1) e^{-i \vec{p_2} \cdot \vec{x}} e^{-i \vec{p_3} \cdot \vec{y}} (2\pi)^d \delta^{(d)}( \vec{p_1} - \vec{p_2} - \vec{p_3} ) \nonumber\\
    & = 2 \lambda_1 \int_0^t dt_1 \iint  \frac{d^d\vec{p_2}}{(2\pi)^d} \frac{d^d\vec{p_3}}{(2\pi)^d} \mathcal{G}(\vec{p_2} + \vec{p_3},t_1) \mathcal{G}(\vec{p_2},t-t_1) \mathcal{G}(\vec{p_3},t-t_1) e^{-i \vec{p_2} \cdot \vec{x}} e^{-i \vec{p_3} \cdot \vec{y}}.
\end{align}
 \end{widetext}
Note that because of the appearance of Dirac delta functions enforcing momentum conservation at each vertex, it becomes useful to indicate on the Feynman graphs the directions of the flowing momenta.

To proceed further with the calculation, one has to pay special attention to the initial condition chosen for the neutron population. In the case of neutrons evolving in an homogeneous and infinite system, and setting a uniform initial neutron density $N_0$, one can use the translation invariance of the problem to set $\vec{p_1} = \vec{0}$. Momentum conservation then implies that Eq.~\ref{spatial_corr_func} be transformed into
\begin{eqnarray}
    \langle N(\vec{x},t) N(\vec{y},t) & &\rangle_c  =  2  \lambda_1 N_0
     \int_0^t dt_1 \int  \frac{d^d\vec{p}}{(2\pi)^d} \mathcal{G}(\vec{0},t_1) \times \nonumber\\
    & &\mathcal{G}(\vec{p},t-t_1) \mathcal{G}(-\vec{p},t-t_1) e^{-i \vec{p} \cdot (\vec{x} - \vec{y})}.
\end{eqnarray}
Performing the integration over the momentum $\vec{p}$ in $d$-dimension, one obtains,
\begin{equation}
    \langle N(\vec{x},t) N(\vec{y},t) \rangle_c  = 2 \lambda_1 N_0 e^{D r t} \int_0^t  \frac{e^{Drt'}}{(8\pi D t')^{d/2}} \, e^{ - \frac{|\vec{x} - \vec{y}|^2}{ 8 D t'} } dt' ,
\end{equation}
in agreement with the literature \cite{mulatier2015}.

\section{From neutronics to reactor physics - addition of a simple feedback mechanism }
\label{sec:directed_percolation}

Reactor physics is neutron transport theory supplemented by feedback mechanisms, which happen whenever the reactor is operated at high power.
These feedback mechanisms are a central aspect of the stability of nuclear reactors but their detailed modelling is generally difficult to implement in a neutron balance, principally because they cannot be associated to a microscopic stochastic process: for example, the Doppler effect is proportional to local neutron densities, since neutrons heat up their surrounding media. In the response functional formalism, these effect can however be directly incorporated at the macroscopic level of description, i.e. in the mean field equations.

In the following, this advantage is therefore exploited to explore a simplified but realistic model of a reactor in the presence of counter-reactions. It will be shown that the feedback mechanism imposed on the neutron population has a strong impact on the structure of the fluctuations close to the percolation transition. The field theoretic viewpoint adopted here will in particular allow to establish that the neutron gas undergoes, whenever the reactor is exactly critical, a phase transition whose universality class will be identified.

\subsection{A simplified model with Doppler broadening}

One of the most paradigmatic feedback effect affecting the neutron population in a reactor core is the Doppler broadening of the capture resonances of $^{238}$U, which couples the neutron population to the temperature of the fuel. As the temperature increases, more and more neutrons are captured by the medium and disappear: the overall neutron population is pulled back towards a reference operating power. This feedback is quite difficult to include in a detailed microscopic model because it involves the reactor power or the temperature of the medium which are macroscopic quantities. It can however be straightforwardly implemented under the form of a restoring force of the type $\sigma N(N - N_{\text{ref.}})$ in the diffusion equation \cite{bell_glasstone}. The behavior of the neutron population evolving in a nuclear reactor operating at nominal power (and thus subject to Doppler broadening) can consequently be described by the equation
\begin{equation}
    \label{stoch_diffusion_doppler}
    \frac{\partial }{ \partial t} N(\vec{x},t) - D ( \Delta + r)  N(\vec{x},t)  +  \sigma N^2(\vec{x},t)  = 0,
\end{equation}
where the shifted parameter $ r = (\rho + \sigma N_{\text{ref.}})/D$ has been introduced, for convenience.

The response functional and the resulting theory to which this equation is associated has been extensively studied. It corresponds to the functional of the \textit{Directed Percolation} (DP) universality class \cite{janssen2004,tauber_2005}:
\begin{equation}
    \mathcal{S} = \int d^d \vec{x} \int dt \ensuremath{\widetilde{N}} \left( \partial_t - D ( \Delta + r)  \right) N + \sigma \ensuremath{\widetilde{N}} N^2 - \lambda \ensuremath{\widetilde{N}}^2 N  .
\end{equation}

The inclusion of the term $\sigma N^2$ has a strong impact on the structure of the theory, for now the two following processes are allowed
\begin{equation*}
    2 \lambda \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of a](o1);
    \diagram* {
    (f1) --[anti fermion] (a) --[fermion] (f2),
    (a) --[anti fermion] (o1)
    };
    \end{feynman}
  \end{tikzpicture} \hspace{2 cm} 2 \sigma
      \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above right=of a](f1);
      \vertex [below right=of a](f2);
      \vertex [left=of a](o1);
    \diagram* {
    (f1) --[fermion] (a) --[anti fermion] (f2),
    (a) --[fermion] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
\end{equation*}
as well as their combination, i.e. the loop diagram
\begin{equation}
    4 \sigma \lambda \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (z);
      \vertex [right=of z](a);
      \vertex [right=of a] (b);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of b](o1);
    \diagram* {
     (z) --[anti fermion] (a)
     -- [fermion,half left] (b)
     -- [fermion,half left] (a),
    (b) --[anti fermion] (o1)
    };
    \end{feynman}
  \end{tikzpicture} \, ,
\end{equation}
whose contribution leads to badly divergent behaviors of two different types \cite{tauber_2005}:
\begin{itemize}
    \item a divergence is observed in dimension $d \geq 2$, because the momentum  associated to the loop is not constrained by any conservation rule. Its integration thus runs over theoretically infinite momenta and include small distances, short time, unbounded fluctuations of the neutron population. This is a paradigmatic example of an Ultra-Violet (UV) divergence;

    \item another and more serious divergence occurs in dimension $d \leq 4$. An Infrared (IR) divergence is observed in the long time, long distance limit. The highest dimension $d_c = 4$ for which the divergence is observed is called the \textit{critical dimension} of the problem.
\end{itemize}

Both types of divergences can be eliminated through the well-known renormalization procedure \cite{Delamotte_2004}, which, in essence, consists in the application of successive spatio-temporal scaling transformations on the problem, so that the small distance and UV divergent fluctuations be truly integrated over. The true macro-scale description of the problem is obtained at a fixed point of the iterated scaling process which generally has the induced benefit of taming the residual IR divergence of the problem.
\begin{comment}
\begin{itemize}
    \item The asymptotic behavior of the theory does not depend on the fine grain details and, in particular, on the precise origin of the phenomena at play. The resulting theory is \textit{universal}.

    \item Because the theory resulting from the renormalization procedure is derived at a fixed point of a spatio-temporal scaling transformation, the observables extracted from it will typically present power law behaviors. An imprint of the system's fluctuations will generally subsist, through characteristic \textit{critical exponents} that are unique to each universality class.
\end{itemize}

\end{comment}

Although the starting point of the model are the mean field equations of the problem, the renormalized observables one ultimately extract from the theory present asymptotic behaviors (i.e. critical exponents) that differ from the predictions of the mean field equations, at least in the vicinity of the phase transition. Close to the critical point, mean field theory breaks down and only the renormalized field theoretic description of the system can correctly describe the evolution of observables. \newline

\subsection{The directed percolation transition}

The largest deviation from the mean field predictions is obtained in the 0D model. Indeed, setting $D =0$ in Eq.~\ref{stoch_diffusion_doppler}, the point-like DP model can be defined as
\begin{equation}
    \frac{dN}{dt} = r N - \sigma N^2 .
\end{equation}
As long as $r>0$, this equation possesses a positive stationary solution $N_{\infty} \propto r^\beta$ ($\beta = 1$ is the first critical exponent of directed percolation). A detailed analysis with the help of the RG leads to a strikingly different prediction \cite{Canet2006,Benitez2013}: in zero dimension, whatever the values chosen for $r$ and $\sigma$, the neutron population always asymptotically vanish. This finding is in sharp contrast with the solution of the mean field equations but has been tested and validated against Monte-Carlo simulations.
Translated into the domain of nuclear reactor physics, this means that adding a feedback term of the form $\sigma N^2$ into the point kinetics equation leads to a theory that do not possess a stable, non-trivial, equilibrium solution:  every point-like reactors subject to counter-reactions are sub-critical~!

The triviality of this point-like model can be easily circumvented though, by considering instead the case of a tri-dimensional, infinite and homogeneous reactor. The system then do admit $N_{\infty} \propto r^\beta$ as a solution, but the RG predicts a value for the exponent $\beta$ close to 0.8 in 3$d$ \cite{tauber_2005, janssen2004} (as compared to 1 for mean field). Although conceptually important, this deviation in the scaling relation between the power level and the reactivity of the reactor would be difficult to observe in practice.  \newline

Even though the feedback mechanism proposed in this paper might appear as a simplistic approximation of the true processes at play in an operating reactor, the asymptotic scaling laws and the overall behaviors predicted by the RG analysis of the directed percolation process might still persist in a more refined treatment. The \textit{DP conjecture} \cite{Hinrichsen2000,janssen2004} indeed states that the DP universality class is the most general universality class observed in transitions between an active to an absorbing state and persist under very general conditions. The modelling of critical nuclear reactors might thus fall into the general class of the DP percolation transitions, even with a more elaborate modelling of counter-reaction mechanisms.

Observable consequences of the DP transition and direct violation of the mean field equation predictions might however be difficult to observe in practice. Apart from the exponent $\beta$, two other critical exponents characterize the DP universality class. They are associated to the correlation length $\xi$ and the diffusion time scale $\tau$ that naturally appear in equation \ref{stoch_diffusion_doppler}:
\begin{equation}
\label{scalinglaws}
    \xi \sim |r|^{-\nu} \; \; \text{ and } \; \; \tau \sim \xi^z/D \sim |r|^{-z \nu}.
\end{equation}
The DP class is entirely characterized by these exponents and any observable consequence of the DP transition can be traced back to a deviation of these exponents from their mean field predictions. Deviations from the mean field theory are expected to be larger for system in lower dimensions, in models farther away from the critical dimension $d_c=4$. Table \ref{critical_exponent_DP} list the values of the exponents of directed percolation, extracted from numerical simulations, for systems in 1, 2 or 3 dimensions.
\begin{center}
\begin{table}[h!]
\caption{Critical exponents of the DP phase transition \cite{Hinrichsen2008}. \label{critical_exponent_DP} }
\begin{tabular}{|c | c | c| c | c|}
 \hline
 exponent & d = 1  & d = 2 & d = 3 & mean field \\
 \hline\hline
 $\beta$ & 0.276486(8) & 0.5834(30)  & 0.813(9) & 1  \\
 \hline
 $\nu$ & 1.096854(4) & 0.7333(75) & 0.584(5) & 0.5  \\
 \hline
 $z$ & 1.580745(10) & 1.7660(16) & 1.901(5) & 2  \\
 \hline
\end{tabular}
\end{table}
\end{center}

A detailed account of all of the possible experimental evidences of the DP transition in nuclear reactors would by far exceed the scope of this paper. As a matter of illustration, the following section will describe what observable and concrete consequence the dynamic exponent $z$ might have in a real reactor.

\begin{comment}
\sout{The practical observation of deviations from mean field theory in a nuclear reactor experiment would be challenging because of the typical experimental resolution would not be precise enough to discriminate with certainty between the mean field or the directed percolation prediction. A detailed account of all of the possible experimental evidences of the DP transition in nuclear reactors would by far exceed the scope of this paper (CONCLURE: faites pas chier). \newline

It might also appear important to mention some other effects that could change the conclusions of this section and might thus require some careful and dedicated treatment.
One very important ingredient towards a more realistic modeling of a nuclear reactors would be to take into account and incorporate into the model a finite size associated to the reactor. Such type of models have been partly studied in the case of the DP transition (REF finite size effect) but their manipulation is generally uneasy. }
Energetics might be quite easily implemented, at least for a multi-group treatment, and would most certainly fall into the category of multi-species directed percolation (REF) transitions. In this respect, a model with precursors might be more interesting to study, simply because precursors are insensitive to the restoring force of the model, but remain tightly linked to the neutron population.
\end{comment}

\subsection{A Consequence on reactor physics}

Reactor physics aims at accurately calculating (with a precision of $\sim 1\%$) the spatial distribution of the neutron flux, from which diverse measurable local quantities can be extracted, such as fission rates, neutron counts, etc. Those quantities are measured either by in-core detection devices (e.g. the aeroball measurement system) or by ex-core detection devices (using source-range, intermediate-range or power-range neutron flux instrumentations). Any significant deviation from the flux measured by these detectors to 'expected' values should trigger various safety warnings notified by the automatic protection system, ultimately leading to the core tripping procedure. Noticeably these expected values of the flux are set by an online solving of a simplified -yet robust- neutron transport model.

During the intermediate and power range operation of the reactor, the thermal feedback effects have to be taken into account as they modify the cross-sections dictating the laws of propagation of the neutrons. As seen in the previous section, this modification has another consequence: the macroscopic properties of the medium emerge from an intimate coupling at the mesoscopic scale between the random walk of the neutrons and the feedback reactions depending on the neutrons local densities. Usual coupling schemes relying on mean field equations are therefore not valid anymore since the mean field equations themselves are wrong, and the properties of the system can only be obtained through the critical exponents of the directed percolation defined in the previous section. In the following, since we are interested in a simplified model of reactor physics, we will restrict ourselves to $d=3$. In this case, the value of $z$, the dynamic critical exponent, is about $1.9$. This exponent describes how fast any local perturbation spreads \cite{Tauber2012CriticalBehavior} and can be related to a modification of the diffusive ($z=2$) random walk which becomes slightly super-diffusive for $d=3$ (indeed $z$ can be interpreted as an anisotropy parameter). The scaling of the spatial correlation length versus the time correlation length follows $\xi^z \sim D \tau$ (see Eq.~\ref{scalinglaws}, with $z=1.901$ in the following). These Lévy flights (see for instance \cite{Lischke2020WhatResults} for a review on fractional Laplacians and their connection to Lévy flights) call for a more quantitative modelling of the mean-field Eq.~\ref{stoch_diffusion_doppler} which is hence modified according to:
\begin{equation}
    \label{stoch_diffusion_doppler_z}
     D (-(-\Delta)^{z/2} + r)  N(x,t)  -  \sigma N^2(x,t)  = \frac{\partial }{ \partial t} N(x,t),
\end{equation}
where $-(-\Delta)^{z/2}$ is the fractional Laplacian defined through $-(-\Delta)^{z/2}=\frac{\partial^z}{\partial |x|^z}$ ($\frac{\partial^z}{\partial |x|^z}$ being the Riesz-Feller derivative of fractional order $z$).
Yet, whenever the reactor is in a stationary state, the Doppler effect does hardly comes into action: it mainly stabilizes the core upon local perturbations or control long-term statistical fluctuations \cite{Dumonteil2021PatchyReactions}. Since the critical exponents obtained in the previous section did not explicitly depend on $\sigma$, representing the strength of the coupling between the neutron scalar field and the feedback mechanism, this equation can further be simplified by considering the parametric regime where this coupling is weak. The feedback term $\sigma N^2$ in the mean field equation can consequently be safely neglected. In this case the stationary behavior of Eq.~\ref{stoch_diffusion_doppler_z} boils down to
\begin{equation}
    \label{stoch_diffusion_doppler_nofeedback}
     D(-\Delta)^{z/2} N(x,t) = r  N(x,t).
\end{equation}
This critical diffusion equation essentially tells that, in the frame of the so-called one-group diffusion equation, moving from neutron transport theory (without feedback) to reactor physics theory (i.e. neutronics with feedback) mainly consists in replacing the standard Laplacian ($z=2$) by the fractional Laplacian with $z=1.901$. The equilibrium spatial distribution of the neutron field can therefore be obtained through the analysis of the first eigenmode of this fractional Laplacian. Since the reactor has a finite spatial extent, it is necessary to estimate this eigenmode in a bounded domain, keeping in mind that the critical exponents obtained through the use of the renormalization group are valid only for an infinite medium. To avoid the intricate analysis of critical exponents of the directed percolation in finite domain, we will anyway use this (infinite medium) dynamic critical exponent while keeping in mind that the results of the analysis will be only valid within the bulk of the medium (i.e. far from the boundaries). The precise topic of the properties of the fractional Laplacian in bounded domains has been investigated in the seminal paper of Zoia et al \cite{Zoia2007FractionalDomains}. This paper proposes different way to numerically tackle such equations. One possible -and elegant- way being to perform a Monte Carlo simulation of the neutron transport where the standard diffusion kernel is replaced by Lévy flights. Given the different hypotheses of our reactor model (diffusion instead of transport, no energy dependence, uniform medium, restriction to the bulk of the domain,...) we are only interested in approximate solutions, and therefore another possibility is to have recourse to approximate analytical expression as the one given by \cite{Kwasnicki2012EigenvaluesInterval} for example. Figure \ref{compa_z} reports the comparison between the exact solution of the critical diffusion equation that models neutron transport and such an approximate solution of our critical super-diffusion equation that models reactor physics. This approximate solution has the form
\begin{equation}
    \label{1DLévyFlight}
     N(x) = A \cos \Bigl( \left( \frac{\pi}{2}-(2-z)\frac{\pi}{8}\right) \frac{x}{L} \Bigr)
\end{equation}
with A being a normalization constant, and is obtained for a $1d$ uniform medium by imposing Dirichlet boundary conditions ($N(\pm L) = 0$) at the boundaries $L$ of the domain for $z=2$. For $z<2$, the neutron density at the boundaries does not vanish: this can be interpreted as long-distance flight bringing some neutrons to cross the boundary instead of being absorbed at it. But in the present case we are not interested in the properties of the medium close to the boundaries (where the dynamic critical coefficient is not valid anymore), hence the proper boundaries conditions for $z=1.901$ are not discussed in this paper. It is however worth noting that the topic of boundaries conditions for Lévy flights and their physics interpretation is currently under active investigation (see \cite{Constantin} for instance).
Equation Eq.~\ref{1DLévyFlight} can also be extended to $3d$ medium by multiplication of the modes along each axis since  $N(x,y,z)=N(x)N(y)N(z)$ (this separation of variables stands for homogeneous partial differential equations in general \cite{Duderstadt1979TransportTheory}), making the study of a $1d$ medium meaningful. Also, for $z$ values close to $2$ the approximated solution comes close to being exact (the $L^2$ norm between the exact and the approximate solution is controled by $(2-z)$ \cite{Kwasnicki2012EigenvaluesInterval}). When $z=2$, this approximate solution becomes exact since one retrieves the usual solution of the critical diffusion equation $N(x) \sim \cos \Bigl( \frac{\pi}{2} \frac{x}{L} \Bigr)$. The $1d$ comparison between the spatial distributions of $z=2$ and $z=1.901$ (top plot of Figure \ref{compa_z}) finally reveals that, as could be qualitatively guessed, the super-diffusion slightly flatten the spatial distribution compared to diffusion. This macroscopic point of view emerges from the fact that the far-distance events characteristics of Lévy flights indeed tend to spatially redistribute the neutrons. In the frame of all the aforementioned hypotheses, the differences between the diffusion and the $z=1.901$ super-diffusion on local neutron densities are of the order of $1-2\%$ (bottom plot of Figure \ref{compa_z}), which is precisely the accuracy targeted by high precision reactor physics codes and methods, thus making these relative differences non-negligible.

\begin{figure}[!ht]
\centering
\hspace{-1 cm}
\includegraphics[scale=0.4]{Figures/compa_z.png}
\caption{ Top plot: comparison between the first normalized eigenfunction of the linear diffusion equation ($z=2$ in plain line) with Dirichlet boundary conditions in $-1$ and $1$, and an approximation of the first normalized eigenfunction of the fractional Laplacian (as given by \cite{Kwasnicki2012EigenvaluesInterval}) with $z=1.901$ (dashed line). The approximation is valid far from the boundaries, i.e. within the bulk. Bottom plot: relative differences (in \%) between both eigenfunctions.}
\label{compa_z}
\end{figure}

\section{Conclusion}
\label{sec:conclusion}

The response functional formalism of Janssen and De Dominicis has proven to be a versatile and powerful companion for exploring the physics of branching processes. Applied to the well known equations of neutron transport in fissile media (neutronics), we have been able to recover with relative ease the main recent results of the field : moments of the neutron population within the reactor, effect of external neutron sources on the chain reaction, taking into account of delayed neutrons, and characterization of spatial and temporal correlations that may develop at low power. Along this path, the close connection between the statistical field theory of Janssen and De Dominics and  the coherent path integral method of Doi and Peliti has been underlined. This last approach is based on establishing a second quantized version of the neutron's branching process, and therefore relies on the solid-ground microscopic description of branching spatial Brownian motion that has been build in the past few decades. But the generality of the response functional approach lend itself to the inclusion of very general ad-hoc models whose microscopic origin would be difficult to implement. The restoring force that has been used in section \ref{sec:directed_percolation} is an example of such a heuristic approach. It allowed the study of a simplistic model of nuclear reactor, that couples neutronics to a prototype feedback model (the Doppler effect), which is usually not possible via the traditional master equations approach. Showing that the phase transition of the neutron population at criticality falls into the directed percolation universality class permitted to retrieve the critical exponents of the phase transition obtained thanks to the well established results of the renormalization group technics. \newline

The main consequence of this work in the realm of nuclear reactor physics were evaluated by focussing on the corrected mean field equation of our simplified reactor physics model, using in particular the dynamic critical exponent $z$ previously obtained. In the same manner that this coefficient transforms the diffusion law $z=2$ into a super-diffusive Lévy flights with $z=1.901$ for $3d$ medium, the usual critical diffusion equation used to grasp the main features of neutronics was shown to be replaced, in this prototype model of reactor physics, by a critical super-diffusion equation. In the low feedback regime, and far from the boundaries, the main eigenmode of this critical super-diffusion equation was approximated using recent developments on the spectrum of the fractional Laplacian in bounded domains. In this frame, the solution of the standard mean field equation for the neutron local distributions might show deviations up to $1-2\%$, which invalidates many formal approaches and reactor physics analytical models (see \cite{Dellacherie, Tomatis} among many others) sometimes used to benchmark reactor physics codes, while questioning the accuracy targeted by the current high precision reactor physics codes themselves, which aim at calculating neutron fluxes with an accuracy below this precision (see \cite{Demaziere} for instance). \newline

Finally, to further assess the accuracy of mean field approaches (and hence reactor physics codes), it would be necessary to investigate the properties of the directed percolation in finite medium, while studying the effects of neutron energy (in the frame of a 2-group critical super-diffusion equation for instance) and improving the modeling of the feedback mechanism. A prominent role might also be played by the precursor of delayed neutrons. Indeed, these precursors being not subject to any feedback mechanism, they might inflect the system's dynamic into a different universality class.
Ultimately, reactor physics experiments aiming at measuring and detecting distortions of the neutron field could definitely sign this phenomenon of neutron percolation in reactor cores.
\section{Analytical calculation of the generalized propagators $\Sigma_{ab}$ ($a,b$ = $N,M$) in the two species model}
\label{app_a}

Focusing on the evaluation of $\Sigma_{NN}(t,t')$, and using the shorthand notation $a = r + \lambda_D$, the diagram of Eq.~\ref{sig_nn} translates into the infinite series,
\begin{widetext}
\begin{equation}
    \Sigma_{NN}(t,t') = e^{ r (t - t')} \left\{  1 + \beta \lambda_D \int_{t'}^{t} dt_1 e^{-a t_1} \int_{t'}^{t_1} dt_2 e^{a t_2} +
    (\beta \lambda_D)^2 \int_{t'}^{t} dt_1 e^{-a t_1} \int_{t'}^{t_1} dt_2 e^{a t_2} \int_{t'}^{t_2} dt_3 e^{-a t_3} \int_{t'}^{t_3} dt_4 e^{a t_4} + ... \right\},
\end{equation}
or,
\end{widetext}
\begin{equation}
   \Sigma_{NN}(t,t') = e^{ r (t - t')} \left( 1 + \sum_{n \geq 1} (\beta \lambda_D)^{n}  A_{n}(t,t') \right),
\end{equation}
with,
\begin{eqnarray}
   & & A_n(t, t') =\nonumber\\
   & & \hspace{-0.2 cm} \int_{t'}^t\int_{t'}^{t_1} ... \int_{t'}^{t_{2n-1}} \hspace{-0.1 cm} dt_1 ...  dt_{2n} e^{-a(t_1 -t_2 + ... +t_{2n-1} - t_{2n})}
\end{eqnarray}
One can first easily show that $A_n(t,t') = A_n(t-t',0)$, so that only $A_n(t) := A_n(t,0)$ needs to be evaluated. Then, the form of the exponential calls for a change of variable of the form $t_1, ..., t_{2n} \rightarrow t_1 -t_2, ..., t_{2n-1} - t_{2n}, t_{2n} = u_{1}, ..., u_{2n}$. Using this new set of variables, $A_n(t)$ can be written,
\begin{equation}
     A_n(t) = \idotsint\limits_{ \substack{0 < u_1, \hdots, u_{2n}\\  u_1 + .. + u_{2n} < t}} du_1 \hdots  du_{2n}  e^{-a(u_1 + \hdots + u_{2n-1})},
\end{equation}
where the exponential now only depends on the odd terms $u_{2i+1}$. The constraint  $\sum_{i=1}^{2n}u_i < t$ can be enforced through the introduction of an indicator function $\mathbbm{1}_{\sum_{i=1}^{2n}u_i < t}$, so that,
\begin{eqnarray}
\label{appA:int_to_split}
    & & A_n(t) =\nonumber\\
    & & \hspace{-0.2 cm}\idotsint\limits_{0 < u_1, \hdots, u_{2n}} du_1 \hdots  du_{2n}  e^{-a(u_1 + \hdots + u_{2n-1})} \mathbbm{1}_{\sum_{i=1}^{2n}u_i < t}
\end{eqnarray}
The integration over even indices corresponds now to a volume integral over an $n$-dimensional simplex,
\begin{equation}
\label{vol_simplex}
     \idotsint\limits_{0 < x_1, \hdots, x_p} dx_1 \hdots  dx_p  \mathbbm{1}_{\sum_{i=1}^{p}x_i < \alpha} = \frac{\alpha^p}{p!}.
\end{equation}
Splitting Eq.~\ref{appA:int_to_split} over even/odd indices, one can then show that,
\begin{widetext}
\begin{equation}
\begin{split}
A_n(t) & = \idotsint\limits_{0 < u_1, \hdots, u_{2n-1}} du_1 \hdots du_{2n-1}  \mathbbm{1}_{\sum_{i=0}^{n-1}u_{2i+1} < t}  e^{-a(u_1 + \hdots + u_{2n-1})}
\idotsint\limits_{0 < u_2, \hdots, u_{2n}} du_2 \hdots du_{2n}  \mathbbm{1}_{\sum_{i=1}^{n}u_{2i} < t - \sum_{i=0}^{n-1}u_{2i+1}} \\
& = \idotsint\limits_{0 < u_1, \hdots, u_{2n-1}} du_1 \hdots du_{2n-1}  \mathbbm{1}_{\sum_{i=0}^{n-1}u_{2i+1} < t}  e^{-a(u_1 + \hdots + u_{2n-1})} \frac{\left(t- u_1 - \hdots u_{2n-1}\right)^n}{n!}
\end{split}
\end{equation}
\end{widetext}
The integrand in the last expression only depends on the sum $u_1 + \hdots + u_{2n-1}$. Thus, defining $x = \sum\limits_{i=0}^{n-1}u_{2i+1}$ and changing the integration variables accordingly, one ends up with,
\begin{equation}
A_n(t) =  \int\limits_{0}^t  dx  e^{-ax} \frac{\left(t-x\right)^n}{n!} \frac{x^{n-1}}{(n-1)!}.
\end{equation}
$\Sigma_{NN}$ can thus be written,
\begin{eqnarray}
\Sigma_{NN}(t,0) & = &  e^{r t}
\Bigg\{ 1 +  \nonumber\\
  & & \hspace{-2 cm}  \left. \sum\limits_{n\geq0} (\beta \lambda_D)^{n+1} \int_0^t dx ~ e^{-ax}  \frac{\left(t-x\right)^{n+1}  x^{n} }{(n+1)! ~ n!} \right\}.
\end{eqnarray}
The latter expression involves the series development of the Modified Bessel function of the first kind $I_1$. Reinstating the lower bound $t'$ into the equation, one then has,
\begin{eqnarray}
\Sigma_{NN}(t, t') & & =  e^{r (t-t')} \Bigg\{ 1 + \nonumber\\
& &  \int_0^{t-t'} d\tau  e^{-r\tau} \sqrt{\frac{\beta \lambda_D(t-t'-\tau)}{\tau}} \nonumber \\
\times & & I_1\left(2\sqrt{\beta \lambda_D \tau (t-t'-\tau)}\right)
\Bigg\}.
\end{eqnarray}
One the same model, one can show that,
\begin{eqnarray}
 \Sigma_{MM}(t,t') & & = e^{-\lambda_D (t-t')} \Bigg\{1 + \nonumber\\
 & &\int_{0}^{t-t'} d\tau  e^{a \tau} \sqrt{\frac{\beta \lambda_D(t-t'-\tau)}{\tau}} \nonumber\\
 & & \times I_1\left(2\sqrt{\beta \lambda_D \tau (t-t'-\tau)}\right) \Bigg\},
\end{eqnarray}
\begin{eqnarray}
     \Sigma_{NM}(t,t') & & = \beta  e^{-\lambda_d(t-t')} \nonumber\\
     & & \hspace{-1 cm}\times \int_{0}^{t-t'}  d\tau ~ e^{a\tau} I_0\left(2\sqrt{\beta \lambda_D \tau (t-t'-\tau)}\right),
\end{eqnarray}
and
\begin{eqnarray}
 \Sigma_{MN}(t,t') & & = \lambda_D ~e^{r(t-t')} \nonumber\\
 & & \hspace{-1 cm}\times \int_{0}^{t-t'} d\tau  e^{-a\tau} I_0\left(2\sqrt{\beta \lambda_D \tau (t-t'-\tau)}\right).
\end{eqnarray}
All these expressions contain integrals that can be further evaluated. After successive change of variables, and setting $\Omega = 4 \beta \lambda_D + \lambda_D^2 + 2 \lambda_D r + r^2$, one can show \cite{gradshteyn} that,
\begin{eqnarray}
    \Sigma_{NN}(t,t') &=& \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \left\{  \sqrt{\Omega} \cosh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right) \right. \nonumber\\
    & & \left.  + (r + \lambda_D) \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right)  \right\},
\end{eqnarray}

\begin{eqnarray}
    \Sigma_{MM}(t,t') &=& \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \left\{  \sqrt{\Omega} \cosh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right) \right. \nonumber\\
    & & \left.  - (r + \lambda_D) \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right)  \right\},
\end{eqnarray}

\begin{equation}
        \Sigma_{NM}(t,t') = 2 \beta \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right),
\end{equation}

\begin{equation}
    \Sigma_{MN}(t,t') = 2 \lambda_D \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right) .
\end{equation}

\section{\label{sec:level1}First-level heading:\protect\\ The line
break was forced \lowercase{via} \textbackslash\textbackslash}

This sample document demonstrates proper use of REV\TeX~4.2 (and
\LaTeXe) in manuscripts prepared for submission to AIP
journals. Further information can be found in the documentation included in the distribution or available at
\url{http://authors.aip.org} and in the documentation for
REV\TeX~4.2 itself.

When commands are referred to in this example file, they are always
shown with their required arguments, using normal \TeX{} format. In
this format, \verb+#1+, \verb+#2+, etc. stand for required
author-supplied arguments to commands. For example, in
\verb+\section{#1}+ the \verb+#1+ stands for the title text of the
author's section heading, and in \verb+\title{#1}+ the \verb+#1+
stands for the title text of the paper.

Line breaks in section headings at all levels can be introduced using
\textbackslash\textbackslash. A blank input line tells \TeX\ that the
paragraph has ended.

\subsection{\label{sec:level2}Second-level heading: Formatting}

This file may be formatted in both the \texttt{preprint} (the default) and
\texttt{reprint} styles; the latter format may be used to
mimic final journal output. Either format may be used for submission
purposes; however, for peer review and production, AIP will format the
article using the \texttt{preprint} class option. Hence, it is
essential that authors check that their manuscripts format acceptably
under \texttt{preprint}. Manuscripts submitted to AIP that do not
format correctly under the \texttt{preprint} option may be delayed in
both the editorial and production processes.

The \texttt{widetext} environment will make the text the width of the
full page, as on page~\pageref{eq:wideeq}. (Note the use the
\verb+\pageref{#1}+ to get the page number right automatically.) The
width-changing commands only take effect in \texttt{twocolumn}
formatting. It has no effect if \texttt{preprint} formatting is chosen
instead.

\subsubsection{\label{sec:level3}Third-level heading: Citations and Footnotes}

Citations in text refer to entries in the Bibliography;
they use the commands \verb+\cite{#1}+ or \verb+\onlinecite{#1}+.
Because REV\TeX\ uses the \verb+natbib+ package of Patrick Daly,
its entire repertoire of commands are available in your document;
see the \verb+natbib+ documentation for further details.
The argument of \verb+\cite+ is a comma-separated list of \emph{keys};
a key may consist of letters and numerals.

By default, citations are numerical; \cite{feyn54} author-year citations are an option.
To give a textual citation, use \verb+\onlinecite{#1}+: (Refs.~\onlinecite{witten2001,epr,Bire82}).
REV\TeX\ ``collapses'' lists of consecutive numerical citations when appropriate.
REV\TeX\ provides the ability to properly punctuate textual citations in author-year style;
this facility works correctly with numerical citations only with \texttt{natbib}'s compress option turned off.
To illustrate, we cite several together \cite{feyn54,witten2001,epr,Berman1983},
and once again (Refs.~\onlinecite{epr,feyn54,Bire82,Berman1983}).
Note that, when numerical citations are used, the references were sorted into the same order they appear in the bibliography.

A reference within the bibliography is specified with a \verb+\bibitem{#1}+ command,
where the argument is the citation key mentioned above.
\verb+\bibitem{#1}+ commands may be crafted by hand or, preferably,
generated by using Bib\TeX.
The AIP styles for REV\TeX~4 include Bib\TeX\ style files
\verb+aipnum.bst+ and \verb+aipauth.bst+, appropriate for
numbered and author-year bibliographies,
respectively.
REV\TeX~4 will automatically choose the style appropriate for
the document's selected class options: the default is numerical, and
you obtain the author-year style by specifying a class option of \verb+author-year+.

This sample file demonstrates a simple use of Bib\TeX\
via a \verb+\bibliography+ command referencing the \verb+aipsamp.bib+ file.
Running Bib\TeX\ (in this case \texttt{bibtex
aipsamp}) after the first pass of \LaTeX\ produces the file
\verb+aipsamp.bbl+ which contains the automatically formatted
\verb+\bibitem+ commands (including extra markup information via
\verb+\bibinfo+ commands). If not using Bib\TeX, the
\verb+thebibiliography+ environment should be used instead.

\paragraph{Fourth-level heading is run in.}%
Footnotes are produced using the \verb+\footnote{#1}+ command.
Numerical style citations put footnotes into the
bibliography\footnote{Automatically placing footnotes into the bibliography requires using BibTeX to compile the bibliography.}.
Author-year and numerical author-year citation styles (each for its own reason) cannot use this method.
Note: due to the method used to place footnotes in the bibliography, \emph{you
must re-run BibTeX every time you change any of your document's
footnotes}.

\section{Math and Equations}
Inline math may be typeset using the \verb+$+ delimiters. Bold math
symbols may be achieved using the \verb+bm+ package and the
\verb+\bm{#1}+ command it supplies. For instance, a bold $\alpha$ can
be typeset as \verb+$\bm{\alpha}$+ giving $\bm{\alpha}$. Fraktur and
Blackboard (or open face or double struck) characters should be
typeset using the \verb+\mathfrak{#1}+ and \verb+\mathbb{#1}+ commands
respectively. Both are supplied by the \texttt{amssymb} package. For
example, \verb+$\mathbb{R}$+ gives $\mathbb{R}$ and
\verb+$\mathfrak{G}$+ gives $\mathfrak{G}$

In \LaTeX\ there are many different ways to display equations, and a
few preferred ways are noted below. Displayed math will center by
default. Use the class option \verb+fleqn+ to flush equations left.

Below we have numbered single-line equations, the most common kind:
\begin{eqnarray}
\chi_+(p)\alt{\bf [}2|{\bf p}|(|{\bf p}|+p_z){\bf ]}^{-1/2}
\left(
\begin{array}{c}
|{\bf p}|+p_z\\
px+ip_y
\end{array}\right)\;,
\\
\left\{%
 \openone234567890abc123\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}%
\label{eq:one}.
\end{eqnarray}
Note the open one in Eq.~(\ref{eq:one}).

Not all numbered equations will fit within a narrow column this
way. The equation number will move down automatically if it cannot fit
on the same line with a one-line equation:
\begin{equation}
\left\{
 ab12345678abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}.
\end{equation}

When the \verb+\label{#1}+ command is used [cf. input for
Eq.~(\ref{eq:one})], the equation can be referred to in text without
knowing the equation number that \TeX\ will assign to it. Just
use \verb+\ref{#1}+, where \verb+#1+ is the same name that used in
the \verb+\label{#1}+ command.

Unnumbered single-line equations can be typeset
using the \verb+\[+, \verb+\]+ format:
\[g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \]

\subsection{Multiline equations}

Multiline equations are obtained by using the \verb+eqnarray+
environment.  Use the \verb+\nonumber+ command at the end of each line
to avoid assigning a number:
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
\delta_{\sigma_1,-\sigma_2}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_jl_i\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1),
\end{eqnarray}
\begin{eqnarray}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\nonumber \\
 & &\times \left( \sum_{i<j}\right)
  \sum_{\text{perm}}
 \frac{1}{S_{12}}
 \frac{1}{S_{12}}
 \sum_\tau c^f_\tau~.
\end{eqnarray}
\textbf{Note:} Do not use \verb+\label{#1}+ on a line of a multiline
equation if \verb+\nonumber+ is also used on that line. Incorrect
cross-referencing will result. Notice the use \verb+\text{#1}+ for
using a Roman font within a math environment.

To set a multiline equation without \emph{any} equation
numbers, use the \verb+\begin{eqnarray*}+,
\verb+\end{eqnarray*}+ format:
\begin{eqnarray*}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\\
 & &\times \left( \sum_{i<j}\right)
 \left(
  \sum_{\text{perm}}\frac{1}{S_{12}S_{23}S_{n1}}
 \right)
 \frac{1}{S_{12}}~.
\end{eqnarray*}
To obtain numbers not normally produced by the automatic numbering,
use the \verb+\tag{#1}+ command, where \verb+#1+ is the desired
equation number. For example, to get an equation number of
(\ref{eq:mynum}),
\begin{equation}
g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \tag{2.6$'$}\label{eq:mynum}
\end{equation}

A few notes on \verb=\tag{#1}=. \verb+\tag{#1}+ requires
\texttt{amsmath}. The \verb+\tag{#1}+ must come before the
\verb+\label{#1}+, if any. The numbering set with \verb+\tag{#1}+ is
\textit{transparent} to the automatic numbering in REV\TeX{};
therefore, the number must be known ahead of time, and it must be
manually adjusted if other equations are added. \verb+\tag{#1}+ works
with both single-line and multiline equations. \verb+\tag{#1}+ should
only be used in exceptional case - do not use it to number all
equations in a paper.

Enclosing single-line and multiline equations in
\verb+\begin{subequations}+ and \verb+\end{subequations}+ will produce
a set of equations that are ``numbered'' with letters, as shown in
Eqs.~(\ref{subeq:1}) and (\ref{subeq:2}) below:
\begin{subequations}
\label{eq:whole}
\begin{equation}
\left\{
 abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}
\right\},\label{subeq:1}
\end{equation}
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1).\label{subeq:2}
\end{eqnarray}
\end{subequations}
Putting a \verb+\label{#1}+ command right after the
\verb+\begin{subequations}+, allows one to
reference all the equations in a subequations environment. For
example, the equations in the preceding subequations environment were
Eqs.~(\ref{eq:whole}).

\subsubsection{Wide equations}
The equation that follows is set in a wide format, i.e., it spans
across the full page. The wide format is reserved for long equations
that cannot be easily broken into four lines or less:
\begin{widetext}
\begin{equation}
{\cal R}^{(\text{d})}=
 g_{\sigma_2}^e
 \left(
   \frac{[\Gamma^Z(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^Z(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)
 + x_WQ_e
 \left(
   \frac{[\Gamma^\gamma(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^\gamma(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)\;. \label{eq:wideeq}
\end{equation}
\end{widetext}
This is typed to show the output is in wide format.
(Since there is no input line between \verb+\equation+ and
this paragraph, there is no paragraph indent for this paragraph.)
\section{Cross-referencing}
REV\TeX{} will automatically number sections, equations, figure
captions, and tables. In order to reference them in text, use the
\verb+\label{#1}+ and \verb+\ref{#1}+ commands. To reference a
particular page, use the \verb+\pageref{#1}+ command.

The \verb+\label{#1}+ should appear in a section heading, within an
equation, or in a table or figure caption. The \verb+\ref{#1}+ command
is used in the text where the citation is to be displayed.  Some
examples: Section~\ref{sec:level1} on page~\pageref{sec:level1},
Table~\ref{tab:table1},%
\begin{table}
\caption{\label{tab:table1}This is a narrow table which fits into a
text column when using \texttt{twocolumn} formatting. Note that
REV\TeX~4 adjusts the intercolumn spacing so that the table fills the
entire width of the column. Table captions are numbered
automatically. This table illustrates left-aligned, centered, and
right-aligned columns.  }
\begin{ruledtabular}
\begin{tabular}{lcr}
Left\footnote{Note a.}&Centered\footnote{Note b.}&Right\\
\hline
1 & 2 & 3\\
10 & 20 & 30\\
100 & 200 & 300\\
\end{tabular}
\end{ruledtabular}
\end{table}
and Fig.~\ref{fig:epsart}.

\section{Figures and Tables}
Figures and tables are typically ``floats''; \LaTeX\ determines their
final position via placement rules.
\LaTeX\ isn't always successful in automatically placing floats where you wish them.

Figures are marked up with the \texttt{figure} environment, the content of which
imports the image (\verb+\includegraphics+) followed by the figure caption (\verb+\caption+).
The argument of the latter command should itself contain a \verb+\label+ command if you
wish to refer to your figure with \verb+\ref+.

Import your image using either the \texttt{graphics} or
\texttt{graphix} packages. These packages both define the
\verb+\includegraphics{#1}+ command, but they differ in the optional
arguments for specifying the orientation, scaling, and translation of the figure.
Fig.~\ref{fig:epsart}%
\begin{figure}
\includegraphics{fig_1
\caption{\label{fig:epsart} A figure caption. The figure captions are
automatically numbered.}
\end{figure}
is small enough to fit in a single column, while
Fig.~\ref{fig:wide}%
\begin{figure*}
\includegraphics{fig_2
\caption{\label{fig:wide}Use the \texttt{figure*} environment to get a wide
figure, spanning the page in \texttt{twocolumn} formatting.}
\end{figure*}
is too wide for a single column,
so instead the \texttt{figure*} environment has been used.

The analog of the \texttt{figure} environment is \texttt{table}, which uses
the same \verb+\caption+ command.
However, you should type your caption command first within the \texttt{table},
instead of last as you did for \texttt{figure}.

The heart of any table is the \texttt{tabular} environment,
which represents the table content as a (vertical) sequence of table rows,
each containing a (horizontal) sequence of table cells.
Cells are separated by the \verb+&+ character;
the row terminates with \verb+\\+.
The required argument for the \texttt{tabular} environment
specifies how data are displayed in each of the columns.
For instance, a column
may be centered (\verb+c+), left-justified (\verb+l+), right-justified (\verb+r+),
or aligned on a decimal point (\verb+d+).
(Table~\ref{tab:table4}%
\begin{table}
\caption{\label{tab:table4}Numbers in columns Three--Five have been
aligned by using the ``d'' column specifier (requires the
\texttt{dcolumn} package).
Non-numeric entries (those entries without
a ``.'') in a ``d'' column are aligned on the decimal point.
Use the
``D'' specifier for more complex layouts. }
\begin{ruledtabular}
\begin{tabular}{ccddd}
One&Two&\mbox{Three}&\mbox{Four}&\mbox{Five}\\
\hline
one&two&\mbox{three}&\mbox{four}&\mbox{five}\\
He&2& 2.77234 & 45672. & 0.69 \\
C\footnote{Some tables require footnotes.}
  &C\footnote{Some tables need more than one footnote.}
  & 12537.64 & 37.66345 & 86.37 \\
\end{tabular}
\end{ruledtabular}
\end{table}
illustrates the use of decimal column alignment.)

Extra column-spacing may be be specified as well, although
REV\TeX~4 sets this spacing so that the columns fill the width of the
table.
Horizontal rules are typeset using the \verb+\hline+
command.
The doubled (or Scotch) rules that appear at the top and
bottom of a table can be achieved by enclosing the \texttt{tabular}
environment within a \texttt{ruledtabular} environment.
Rows whose columns span multiple columns can be typeset using \LaTeX's
\verb+\multicolumn{#1}{#2}{#3}+ command
(for example, see the first row of Table~\ref{tab:table3}).%
\begin{table*}
\caption{\label{tab:table3}This is a wide table that spans the page
width in \texttt{twocolumn} mode. It is formatted using the
\texttt{table*} environment. It also demonstrates the use of
\textbackslash\texttt{multicolumn} in rows with entries that span
more than one column.}
\begin{ruledtabular}
\begin{tabular}{ccccc}
 &\multicolumn{2}{c}{$D_{4h}^1$}&\multicolumn{2}{c}{$D_{4h}^5$}\\
 Ion&1st alternative&2nd alternative&lst alternative
&2nd alternative\\ \hline
 K&$(2e)+(2f)$&$(4i)$ &$(2c)+(2d)$&$(4f)$ \\
 Mn&$(2g)$\footnote{The $z$ parameter of these positions is $z\sim\frac{1}{4}$.}
 &$(a)+(b)+(c)+(d)$&$(4e)$&$(2a)+(2b)$\\
 Cl&$(a)+(b)+(c)+(d)$&$(2g)$\footnote{This is a footnote in a table that spans the full page
width in \texttt{twocolumn} mode. It is supposed to set on the full width of the page, just as the caption does. }
 &$(4e)^{\text{a}}$\\
 He&$(8r)^{\text{a}}$&$(4j)^{\text{a}}$&$(4g)^{\text{a}}$\\
 Ag& &$(4k)^{\text{a}}$& &$(4h)^{\text{a}}$\\
\end{tabular}
\end{ruledtabular}
\end{table*}

The tables in this document illustrate various effects.
Tables that fit in a narrow column are contained in a \texttt{table}
environment.
Table~\ref{tab:table3} is a wide table, therefore set with the
\texttt{table*} environment.
Lengthy tables may need to break across pages.
A simple way to allow this is to specify
the \verb+[H]+ float placement on the \texttt{table} or
\texttt{table*} environment.
Alternatively, using the standard \LaTeXe\ package \texttt{longtable}
gives more control over how tables break and allows headers and footers
to be specified for each page of the table.
An example of the use of \texttt{longtable} can be found
in the file \texttt{summary.tex} that is included with the REV\TeX~4
distribution.

There are two methods for setting footnotes within a table (these
footnotes will be displayed directly below the table rather than at
the bottom of the page or in the bibliography).
The easiest
and preferred method is just to use the \verb+\footnote{#1}+
command. This will automatically enumerate the footnotes with
lowercase roman letters.
However, it is sometimes necessary to have
multiple entries in the table share the same footnote.
In this case,
create the footnotes using
\verb+\footnotemark[#1]+ and \verb+\footnotetext[#1]{#2}+.
\texttt{\#1} is a numeric value.
Each time the same value for \texttt{\#1} is used,
the same mark is produced in the table.
The \verb+\footnotetext[#1]{#2}+ commands are placed after the \texttt{tabular}
environment.
Examine the \LaTeX\ source and output for Tables~\ref{tab:table1} and
\ref{tab:table2}%
\begin{table}
\caption{\label{tab:table2}A table with more columns still fits
properly in a column. Note that several entries share the same
footnote. Inspect the \LaTeX\ input for this table to see
exactly how it is done.}
\begin{ruledtabular}
\begin{tabular}{cccccccc}
 &$r_c$ (\AA)&$r_0$ (\AA)&$\kappa r_0$&
 &$r_c$ (\AA) &$r_0$ (\AA)&$\kappa r_0$\\
\hline
Cu& 0.800 & 14.10 & 2.550 &Sn\footnotemark[1]
& 0.680 & 1.870 & 3.700 \\
Ag& 0.990 & 15.90 & 2.710 &Pb\footnotemark[2]
& 0.450 & 1.930 & 3.760 \\
Au& 1.150 & 15.90 & 2.710 &Ca\footnotemark[3]
& 0.750 & 2.170 & 3.560 \\
Mg& 0.490 & 17.60 & 3.200 &Sr\footnotemark[4]
& 0.900 & 2.370 & 3.720 \\
Zn& 0.300 & 15.20 & 2.970 &Li\footnotemark[2]
& 0.380 & 1.730 & 2.830 \\
Cd& 0.530 & 17.10 & 3.160 &Na\footnotemark[5]
& 0.760 & 2.110 & 3.120 \\
Hg& 0.550 & 17.80 & 3.220 &K\footnotemark[5]
&  1.120 & 2.620 & 3.480 \\
Al& 0.230 & 15.80 & 3.240 &Rb\footnotemark[3]
& 1.330 & 2.800 & 3.590 \\
Ga& 0.310 & 16.70 & 3.330 &Cs\footnotemark[4]
& 1.420 & 3.030 & 3.740 \\
In& 0.460 & 18.40 & 3.500 &Ba\footnotemark[5]
& 0.960 & 2.460 & 3.780 \\
Tl& 0.480 & 18.90 & 3.550 & & & & \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{Here's the first, from Ref.~\onlinecite{feyn54}.}
\footnotetext[2]{Here's the second.}
\footnotetext[3]{Here's the third.}
\footnotetext[4]{Here's the fourth.}
\footnotetext[5]{And etc.}
\end{table}
for an illustration.

All AIP journals require that the initial citation of
figures or tables be in numerical order.
\LaTeX's automatic numbering of floats is your friend here:
just put each \texttt{figure} environment immediately following
its first reference (\verb+\ref+), as we have done in this example file.

\begin{acknowledgments}
We wish to acknowledge the support of the author community in using
REV\TeX{}, offering suggestions and encouragement, testing new versions,
\dots.
\end{acknowledgments}

\section{\label{sec:level1}First-level heading:\protect\\ The line
break was forced \lowercase{via} \textbackslash\textbackslash}

This sample document demonstrates proper use of REV\TeX~4.2 (and
\LaTeXe) in manuscripts prepared for submission to AAPM
journals. Further information can be found in the documentation included in the distribution or available at
\url{http://www.aapm.org} and in the documentation for
REV\TeX~4.2 itself.

When commands are referred to in this example file, they are always
shown with their required arguments, using normal \TeX{} format. In
this format, \verb+#1+, \verb+#2+, etc. stand for required
author-supplied arguments to commands. For example, in
\verb+\section{#1}+ the \verb+#1+ stands for the title text of the
author's section heading, and in \verb+\title{#1}+ the \verb+#1+
stands for the title text of the paper.

Line breaks in section headings at all levels can be introduced using
\textbackslash\textbackslash. A blank input line tells \TeX\ that the
paragraph has ended.

\subsection{\label{sec:level2}Second-level heading: Formatting}

This file may be formatted in both the \texttt{preprint} (the default) and
\texttt{reprint} styles; the latter format may be used to
mimic final journal output. In addition, there is another
option available, \texttt{lengthcheck}, which formats the document as closely
as possible to an actual journal article, to facilitate the author's
performance of a length check. Either format may be used for submission
purposes; however, for peer review and production, AAPM will format the
article using the \texttt{preprint} class option. Hence, it is
essential that authors check that their manuscripts format acceptably
under \texttt{preprint}. Manuscripts submitted to AAPM that do not
format correctly under the \texttt{preprint} option may be delayed in
both the editorial and production processes.

The \texttt{widetext} environment will make the text the width of the
full page, as on page~\pageref{eq:wideeq}. (Note the use the
\verb+\pageref{#1}+ to get the page number right automatically.) The
width-changing commands only take effect in \texttt{twocolumn}
formatting. It has no effect if \texttt{preprint} formatting is chosen
instead.

\subsubsection{\label{sec:level3}Third-level heading: Citations and Footnotes}

Citations in text refer to entries in the Bibliography;
they use the commands \verb+\cite{#1}+ or \verb+\onlinecite{#1}+.
Because REV\TeX\ uses the \verb+natbib+ package of Patrick Daly,
its entire repertoire of commands are available in your document;
see the \verb+natbib+ documentation for further details.
The argument of \verb+\cite+ is a comma-separated list of \emph{keys};
a key may consist of letters and numerals.

By default, AAPM citations are numerical; \cite{feyn54}
to give a textual citation, use \verb+\onlinecite{#1}+: (Refs.~\onlinecite{witten2001,epr,Bire82}).
REV\TeX\ ``collapses'' lists of consecutive numerical citations when appropriate.
To illustrate, we cite several together \cite{feyn54,witten2001,epr,Berman1983},
and once again (Refs.~\onlinecite{epr,feyn54,Bire82,Berman1983}).
Note that, when numerical citations are used, the references were sorted into the same order they appear in the bibliography.

A reference within the bibliography is specified with a \verb+\bibitem{#1}+ command,
where the argument is the citation key mentioned above.
\verb+\bibitem{#1}+ commands may be crafted by hand or, preferably,
generated by using Bib\TeX.
The AAPM styles for REV\TeX~4 include Bib\TeX\ style file
\verb+aapmrev4-2.bst+, appropriate for
numbered bibliography.
REV\TeX~4 will automatically choose the style appropriate for
the document's selected class options: the default is numerical.

This sample file demonstrates a simple use of Bib\TeX\
via a \verb+\bibliography+ command referencing the \verb+aapmsamp.bib+ file.
Running Bib\TeX\ (in this case \texttt{bibtex
aapmsamp}) after the first pass of \LaTeX\ produces the file
\verb+aapmsamp.bbl+ which contains the automatically formatted
\verb+\bibitem+ commands (including extra markup information via
\verb+\bibinfo+ commands). If not using Bib\TeX, the
\verb+thebibiliography+ environment should be used instead.

\paragraph{Fourth-level heading is run in.}%
Footnotes are produced using the \verb+\footnote{#1}+ command.
Numerical style citations put footnotes into the
bibliography\footnote{Automatically placing footnotes into the bibliography requires using BibTeX to compile the bibliography.}.
Note: due to the method used to place footnotes in the bibliography, \emph{you
must re-run BibTeX every time you change any of your document's
footnotes}.

\section{Math and Equations}
Inline math may be typeset using the \verb+$+ delimiters. Bold math
symbols may be achieved using the \verb+bm+ package and the
\verb+\bm{#1}+ command it supplies. For instance, a bold $\alpha$ can
be typeset as \verb+$\bm{\alpha}$+ giving $\bm{\alpha}$. Fraktur and
Blackboard (or open face or double struck) characters should be
typeset using the \verb+\mathfrak{#1}+ and \verb+\mathbb{#1}+ commands
respectively. Both are supplied by the \texttt{amssymb} package. For
example, \verb+$\mathbb{R}$+ gives $\mathbb{R}$ and
\verb+$\mathfrak{G}$+ gives $\mathfrak{G}$

In \LaTeX\ there are many different ways to display equations, and a
few preferred ways are noted below. Displayed math will flush left by
default.

Below we have numbered single-line equations, the most common kind:
\begin{eqnarray}
\chi_+(p)\alt{\bf [}2|{\bf p}|(|{\bf p}|+p_z){\bf ]}^{-1/2}
\left(
\begin{array}{c}
|{\bf p}|+p_z\\
px+ip_y
\end{array}\right)\;,
\\
\left\{%
 \openone234567890abc123\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}%
\label{eq:one}.
\end{eqnarray}
Note the open one in Eq.~(\ref{eq:one}).

Not all numbered equations will fit within a narrow column this
way. The equation number will move down automatically if it cannot fit
on the same line with a one-line equation:
\begin{equation}
\left\{
 ab12345678abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}.
\end{equation}

When the \verb+\label{#1}+ command is used [cf. input for
Eq.~(\ref{eq:one})], the equation can be referred to in text without
knowing the equation number that \TeX\ will assign to it. Just
use \verb+\ref{#1}+, where \verb+#1+ is the same name that used in
the \verb+\label{#1}+ command.

Unnumbered single-line equations can be typeset
using the \verb+\[+, \verb+\]+ format:
\[g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \]

\subsection{Multiline equations}

Multiline equations are obtained by using the \verb+eqnarray+
environment.  Use the \verb+\nonumber+ command at the end of each line
to avoid assigning a number:
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
\delta_{\sigma_1,-\sigma_2}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_jl_i\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1),
\end{eqnarray}
\begin{eqnarray}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\nonumber \\
 & &\times \left( \sum_{i<j}\right)
  \sum_{\text{perm}}
 \frac{1}{S_{12}}
 \frac{1}{S_{12}}
 \sum_\tau c^f_\tau~.
\end{eqnarray}
\textbf{Note:} Do not use \verb+\label{#1}+ on a line of a multiline
equation if \verb+\nonumber+ is also used on that line. Incorrect
cross-referencing will result. Notice the use \verb+\text{#1}+ for
using a Roman font within a math environment.

To set a multiline equation without \emph{any} equation
numbers, use the \verb+\begin{eqnarray*}+,
\verb+\end{eqnarray*}+ format:
\begin{eqnarray*}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\\
 & &\times \left( \sum_{i<j}\right)
 \left(
  \sum_{\text{perm}}\frac{1}{S_{12}S_{23}S_{n1}}
 \right)
 \frac{1}{S_{12}}~.
\end{eqnarray*}
To obtain numbers not normally produced by the automatic numbering,
use the \verb+\tag{#1}+ command, where \verb+#1+ is the desired
equation number. For example, to get an equation number of
(\ref{eq:mynum}),
\begin{equation}
g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \tag{2.6$'$}\label{eq:mynum}
\end{equation}

A few notes on \verb=\tag{#1}=. \verb+\tag{#1}+ requires
\texttt{amsmath}. The \verb+\tag{#1}+ must come before the
\verb+\label{#1}+, if any. The numbering set with \verb+\tag{#1}+ is
\textit{transparent} to the automatic numbering in REV\TeX{};
therefore, the number must be known ahead of time, and it must be
manually adjusted if other equations are added. \verb+\tag{#1}+ works
with both single-line and multiline equations. \verb+\tag{#1}+ should
only be used in exceptional case - do not use it to number all
equations in a paper.

Note the equation number gets reset again:
\begin{equation}
g^+g^+g^+ \rightarrow g^+g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~.
\end{equation}

Enclosing single-line and multiline equations in
\verb+\begin{subequations}+ and \verb+\end{subequations}+ will produce
a set of equations that are ``numbered'' with letters, as shown in
Eqs.~(\ref{subeq:1}) and (\ref{subeq:2}) below:
\begin{subequations}
\label{eq:whole}
\begin{equation}
\left\{
 abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}
\right\},\label{subeq:1}
\end{equation}
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1).\label{subeq:2}
\end{eqnarray}
\end{subequations}
Putting a \verb+\label{#1}+ command right after the
\verb+\begin{subequations}+, allows one to
reference all the equations in a subequations environment. For
example, the equations in the preceding subequations environment were
Eqs.~(\ref{eq:whole}).

\subsubsection{Wide equations}
The equation that follows is set in a wide format, i.e., it spans
across the full page. The wide format is reserved for long equations
that cannot be easily broken into four lines or less:
\begin{widetext}
\begin{equation}
{\cal R}^{(\text{d})}=
 g_{\sigma_2}^e
 \left(
   \frac{[\Gamma^Z(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^Z(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)
 + x_WQ_e
 \left(
   \frac{[\Gamma^\gamma(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^\gamma(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)\;. \label{eq:wideeq}
\end{equation}
\end{widetext}
This is typed to show the output is in wide format.
(Since there is no input line between \verb+\equation+ and
this paragraph, there is no paragraph indent for this paragraph.)
\section{Cross-referencing}
REV\TeX{} will automatically number sections, equations, figure
captions, and tables. In order to reference them in text, use the
\verb+\label{#1}+ and \verb+\ref{#1}+ commands. To reference a
particular page, use the \verb+\pageref{#1}+ command.

The \verb+\label{#1}+ should appear in a section heading, within an
equation, or in a table or figure caption. The \verb+\ref{#1}+ command
is used in the text where the citation is to be displayed.  Some
examples: Section~\ref{sec:level1} on page~\pageref{sec:level1},
Table~\ref{tab:table1},%
\begin{table}
\caption{\label{tab:table1}This is a narrow table which fits into a
text column when using \texttt{twocolumn} formatting. Note that
REV\TeX~4 adjusts the intercolumn spacing so that the table fills the
entire width of the column. Table captions are numbered
automatically. This table illustrates left-aligned, centered, and
right-aligned columns.  }
\begin{ruledtabular}
\begin{tabular}{lcr}
Left\footnote{Note a.}&Centered\footnote{Note b.}&Right\\
\hline
1 & 2 & 3\\
10 & 20 & 30\\
100 & 200 & 300\\
\end{tabular}
\end{ruledtabular}
\end{table}
and Fig.~\ref{fig:epsart}.

\section{Figures and Tables}
Figures and tables are typically ``floats''; \LaTeX\ determines their
final position via placement rules.
\LaTeX\ isn't always successful in automatically placing floats where you wish them.

Figures are marked up with the \texttt{figure} environment, the content of which
imports the image (\verb+\includegraphics+) followed by the figure caption (\verb+\caption+).
The argument of the latter command should itself contain a \verb+\label+ command if you
wish to refer to your figure with \verb+\ref+.

Import your image using either the \texttt{graphics} or
\texttt{graphix} packages. These packages both define the
\verb+\includegraphics{#1}+ command, but they differ in the optional
arguments for specifying the orientation, scaling, and translation of the figure.
Fig.~\ref{fig:epsart}%
\begin{figure}
\includegraphics{fig_1
\caption{\label{fig:epsart} A figure caption. The figure captions are
automatically numbered.}
\end{figure}
is small enough to fit in a single column, while
Fig.~\ref{fig:wide}%
\begin{figure*}
\includegraphics{fig_2
\caption{\label{fig:wide}Use the \texttt{figure*} environment to get a wide
figure, spanning the page in \texttt{twocolumn} formatting.}
\end{figure*}
is too wide for a single column,
so instead the \texttt{figure*} environment has been used.

The analog of the \texttt{figure} environment is \texttt{table}, which uses
the same \verb+\caption+ command.
However, you should type your caption command first within the \texttt{table},
instead of last as you did for \texttt{figure}.

The heart of any table is the \texttt{tabular} environment,
which represents the table content as a (vertical) sequence of table rows,
each containing a (horizontal) sequence of table cells.
Cells are separated by the \verb+&+ character;
the row terminates with \verb+\\+.
The required argument for the \texttt{tabular} environment
specifies how data are displayed in each of the columns.
For instance, a column
may be centered (\verb+c+), left-justified (\verb+l+), right-justified (\verb+r+),
or aligned on a decimal point (\verb+d+).
(Table~\ref{tab:table4}%
\begin{table}
\caption{\label{tab:table4}Numbers in columns Three--Five have been
aligned by using the ``d'' column specifier (requires the
\texttt{dcolumn} package).
Non-numeric entries (those entries without
a ``.'') in a ``d'' column are aligned on the decimal point.
Use the
``D'' specifier for more complex layouts. }
\begin{ruledtabular}
\begin{tabular}{ccddd}
One&Two&\mbox{Three}&\mbox{Four}&\mbox{Five}\\
\hline
one&two&\mbox{three}&\mbox{four}&\mbox{five}\\
He&2& 2.77234 & 45672. & 0.69 \\
C\footnote{Some tables require footnotes.}
  &C\footnote{Some tables need more than one footnote.}
  & 12537.64 & 37.66345 & 86.37 \\
\end{tabular}
\end{ruledtabular}
\end{table}
illustrates the use of decimal column alignment.)

Extra column-spacing may be be specified as well, although
REV\TeX~4 sets this spacing so that the columns fill the width of the
table.
Horizontal rules are typeset using the \verb+\hline+
command.
The doubled (or Scotch) rules that appear at the top and
bottom of a table can be achieved by enclosing the \texttt{tabular}
environment within a \texttt{ruledtabular} environment.
Rows whose columns span multiple columns can be typeset using \LaTeX's
\verb+\multicolumn{#1}{#2}{#3}+ command
(for example, see the first row of Table~\ref{tab:table3}).%
\begin{table*}
\caption{\label{tab:table3}This is a wide table that spans the page
width in \texttt{twocolumn} mode. It is formatted using the
\texttt{table*} environment. It also demonstrates the use of
\textbackslash\texttt{multicolumn} in rows with entries that span
more than one column.}
\begin{ruledtabular}
\begin{tabular}{ccccc}
 &\multicolumn{2}{c}{$D_{4h}^1$}&\multicolumn{2}{c}{$D_{4h}^5$}\\
 Ion&1st alternative&2nd alternative&lst alternative
&2nd alternative\\ \hline
 K&$(2e)+(2f)$&$(4i)$ &$(2c)+(2d)$&$(4f)$ \\
 Mn&$(2g)$\footnote{The $z$ parameter of these positions is $z\sim\frac{1}{4}$.}
 &$(a)+(b)+(c)+(d)$&$(4e)$&$(2a)+(2b)$\\
 Cl&$(a)+(b)+(c)+(d)$&$(2g)$\footnote{This is a footnote in a table that spans the full page
width in \texttt{twocolumn} mode. It is supposed to set on the full width of the page, just as the caption does. }
 &$(4e)^{\text{a}}$\\
 He&$(8r)^{\text{a}}$&$(4j)^{\text{a}}$&$(4g)^{\text{a}}$\\
 Ag& &$(4k)^{\text{a}}$& &$(4h)^{\text{a}}$\\
\end{tabular}
\end{ruledtabular}
\end{table*}

The tables in this document illustrate various effects.
Tables that fit in a narrow column are contained in a \texttt{table}
environment.
Table~\ref{tab:table3} is a wide table, therefore set with the
\texttt{table*} environment.
Lengthy tables may need to break across pages.
A simple way to allow this is to specify
the \verb+[H]+ float placement on the \texttt{table} or
\texttt{table*} environment.
Alternatively, using the standard \LaTeXe\ package \texttt{longtable}
gives more control over how tables break and allows headers and footers
to be specified for each page of the table.
An example of the use of \texttt{longtable} can be found
in the file \texttt{summary.tex} that is included with the REV\TeX~4
distribution.

There are two methods for setting footnotes within a table (these
footnotes will be displayed directly below the table rather than at
the bottom of the page or in the bibliography).
The easiest
and preferred method is just to use the \verb+\footnote{#1}+
command. This will automatically enumerate the footnotes with
lowercase roman letters.
However, it is sometimes necessary to have
multiple entries in the table share the same footnote.
In this case,
create the footnotes using
\verb+\footnotemark[#1]+ and \verb+\footnotetext[#1]{#2}+.
\texttt{\#1} is a numeric value.
Each time the same value for \texttt{\#1} is used,
the same mark is produced in the table.
The \verb+\footnotetext[#1]{#2}+ commands are placed after the \texttt{tabular}
environment.
Examine the \LaTeX\ source and output for Tables~\ref{tab:table1} and
\ref{tab:table2}%
\begin{table}
\caption{\label{tab:table2}A table with more columns still fits
properly in a column. Note that several entries share the same
footnote. Inspect the \LaTeX\ input for this table to see
exactly how it is done.}
\begin{ruledtabular}
\begin{tabular}{cccccccc}
 &$r_c$ (\AA)&$r_0$ (\AA)&$\kappa r_0$&
 &$r_c$ (\AA) &$r_0$ (\AA)&$\kappa r_0$\\
\hline
Cu& 0.800 & 14.10 & 2.550 &Sn\footnotemark[1]
& 0.680 & 1.870 & 3.700 \\
Ag& 0.990 & 15.90 & 2.710 &Pb\footnotemark[2]
& 0.450 & 1.930 & 3.760 \\
Au& 1.150 & 15.90 & 2.710 &Ca\footnotemark[3]
& 0.750 & 2.170 & 3.560 \\
Mg& 0.490 & 17.60 & 3.200 &Sr\footnotemark[4]
& 0.900 & 2.370 & 3.720 \\
Zn& 0.300 & 15.20 & 2.970 &Li\footnotemark[2]
& 0.380 & 1.730 & 2.830 \\
Cd& 0.530 & 17.10 & 3.160 &Na\footnotemark[5]
& 0.760 & 2.110 & 3.120 \\
Hg& 0.550 & 17.80 & 3.220 &K\footnotemark[5]
&  1.120 & 2.620 & 3.480 \\
Al& 0.230 & 15.80 & 3.240 &Rb\footnotemark[3]
& 1.330 & 2.800 & 3.590 \\
Ga& 0.310 & 16.70 & 3.330 &Cs\footnotemark[4]
& 1.420 & 3.030 & 3.740 \\
In& 0.460 & 18.40 & 3.500 &Ba\footnotemark[5]
& 0.960 & 2.460 & 3.780 \\
Tl& 0.480 & 18.90 & 3.550 & & & & \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{Here's the first, from Ref.~\onlinecite{feyn54}.}
\footnotetext[2]{Here's the second.}
\footnotetext[3]{Here's the third.}
\footnotetext[4]{Here's the fourth.}
\footnotetext[5]{And etc.}
\end{table}
for an illustration.

All AAPM journals require that the initial citation of
figures or tables be in numerical order.
\LaTeX's automatic numbering of floats is your friend here:
just put each \texttt{figure} environment immediately following
its first reference (\verb+\ref+), as we have done in this example file.

\begin{acknowledgments}
We wish to acknowledge the support of the author community in using
REV\TeX{}, offering suggestions and encouragement, testing new versions,
\dots.
\end{acknowledgments}

\section{\label{sec:level1}First-level heading:\protect\\ The line
break was forced \lowercase{via} \textbackslash\textbackslash}

This sample document demonstrates proper use of REV\TeX~4.2 (and
\LaTeXe) in manuscripts prepared for submission to AIP
journals. Further information can be found in the documentation included in the distribution or available at
\url{http://authors.aip.org} and in the documentation for
REV\TeX~4.2 itself.

When commands are referred to in this example file, they are always
shown with their required arguments, using normal \TeX{} format. In
this format, \verb+#1+, \verb+#2+, etc. stand for required
author-supplied arguments to commands. For example, in
\verb+\section{#1}+ the \verb+#1+ stands for the title text of the
author's section heading, and in \verb+\title{#1}+ the \verb+#1+
stands for the title text of the paper.

Line breaks in section headings at all levels can be introduced using
\textbackslash\textbackslash. A blank input line tells \TeX\ that the
paragraph has ended.

\subsection{\label{sec:level2}Second-level heading: Formatting}

This file may be formatted in both the \texttt{preprint} (the default) and
\texttt{reprint} styles; the latter format may be used to
mimic final journal output. Either format may be used for submission
purposes; however, for peer review and production, AIP will format the
article using the \texttt{preprint} class option. Hence, it is
essential that authors check that their manuscripts format acceptably
under \texttt{preprint}. Manuscripts submitted to AIP that do not
format correctly under the \texttt{preprint} option may be delayed in
both the editorial and production processes.

The \texttt{widetext} environment will make the text the width of the
full page, as on page~\pageref{eq:wideeq}. (Note the use the
\verb+\pageref{#1}+ to get the page number right automatically.) The
width-changing commands only take effect in \texttt{twocolumn}
formatting. It has no effect if \texttt{preprint} formatting is chosen
instead.

\subsubsection{\label{sec:level3}Third-level heading: Citations and Footnotes}

Citations in text refer to entries in the Bibliography;
they use the commands \verb+\cite{#1}+ or \verb+\onlinecite{#1}+.
Because REV\TeX\ uses the \verb+natbib+ package of Patrick Daly,
its entire repertoire of commands are available in your document;
see the \verb+natbib+ documentation for further details.
The argument of \verb+\cite+ is a comma-separated list of \emph{keys};
a key may consist of letters and numerals.

By default, citations are numerical; \cite{feyn54} author-year citations are an option.
To give a textual citation, use \verb+\onlinecite{#1}+: (Refs.~\onlinecite{witten2001,epr,Bire82}).
REV\TeX\ ``collapses'' lists of consecutive numerical citations when appropriate.
REV\TeX\ provides the ability to properly punctuate textual citations in author-year style;
this facility works correctly with numerical citations only with \texttt{natbib}'s compress option turned off.
To illustrate, we cite several together \cite{feyn54,witten2001,epr,Berman1983},
and once again (Refs.~\onlinecite{epr,feyn54,Bire82,Berman1983}).
Note that, when numerical citations are used, the references were sorted into the same order they appear in the bibliography.

A reference within the bibliography is specified with a \verb+\bibitem{#1}+ command,
where the argument is the citation key mentioned above.
\verb+\bibitem{#1}+ commands may be crafted by hand or, preferably,
generated by using Bib\TeX.
The AIP styles for REV\TeX~4 include Bib\TeX\ style files
\verb+aipnum.bst+ and \verb+aipauth.bst+, appropriate for
numbered and author-year bibliographies,
respectively.
REV\TeX~4 will automatically choose the style appropriate for
the document's selected class options: the default is numerical, and
you obtain the author-year style by specifying a class option of \verb+author-year+.

This sample file demonstrates a simple use of Bib\TeX\
via a \verb+\bibliography+ command referencing the \verb+sorsamp.bib+ file.
Running Bib\TeX\ (in this case \texttt{bibtex
sorsamp}) after the first pass of \LaTeX\ produces the file
\verb+sorsamp.bbl+ which contains the automatically formatted
\verb+\bibitem+ commands (including extra markup information via
\verb+\bibinfo+ commands). If not using Bib\TeX, the
\verb+thebibiliography+ environment should be used instead.

\paragraph{Fourth-level heading is run in.}%
Footnotes are produced using the \verb+\footnote{#1}+ command.
Numerical style citations put footnotes into the
bibliography\footnote{Automatically placing footnotes into the bibliography requires using BibTeX to compile the bibliography.}.
Author-year and numerical author-year citation styles (each for its own reason) cannot use this method.
Note: due to the method used to place footnotes in the bibliography, \emph{you
must re-run BibTeX every time you change any of your document's
footnotes}.

\section{Math and Equations}
Inline math may be typeset using the \verb+$+ delimiters. Bold math
symbols may be achieved using the \verb+bm+ package and the
\verb+\bm{#1}+ command it supplies. For instance, a bold $\alpha$ can
be typeset as \verb+$\bm{\alpha}$+ giving $\bm{\alpha}$. Fraktur and
Blackboard (or open face or double struck) characters should be
typeset using the \verb+\mathfrak{#1}+ and \verb+\mathbb{#1}+ commands
respectively. Both are supplied by the \texttt{amssymb} package. For
example, \verb+$\mathbb{R}$+ gives $\mathbb{R}$ and
\verb+$\mathfrak{G}$+ gives $\mathfrak{G}$

In \LaTeX\ there are many different ways to display equations, and a
few preferred ways are noted below. Displayed math will center by
default. Use the class option \verb+fleqn+ to flush equations left.

Below we have numbered single-line equations, the most common kind:
\begin{eqnarray}
\chi_+(p)\alt{\bf [}2|{\bf p}|(|{\bf p}|+p_z){\bf ]}^{-1/2}
\left(
\begin{array}{c}
|{\bf p}|+p_z\\
px+ip_y
\end{array}\right)\;,
\\
\left\{%
 \openone234567890abc123\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}%
\label{eq:one}.
\end{eqnarray}
Note the open one in Eq.~(\ref{eq:one}).

Not all numbered equations will fit within a narrow column this
way. The equation number will move down automatically if it cannot fit
on the same line with a one-line equation:
\begin{equation}
\left\{
 ab12345678abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}.
\end{equation}

When the \verb+\label{#1}+ command is used [cf. input for
Eq.~(\ref{eq:one})], the equation can be referred to in text without
knowing the equation number that \TeX\ will assign to it. Just
use \verb+\ref{#1}+, where \verb+#1+ is the same name that used in
the \verb+\label{#1}+ command.

Unnumbered single-line equations can be typeset
using the \verb+\[+, \verb+\]+ format:
\[g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \]

\subsection{Multiline equations}

Multiline equations are obtained by using the \verb+eqnarray+
environment.  Use the \verb+\nonumber+ command at the end of each line
to avoid assigning a number:
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
\delta_{\sigma_1,-\sigma_2}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_jl_i\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1),
\end{eqnarray}
\begin{eqnarray}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\nonumber \\
 & &\times \left( \sum_{i<j}\right)
  \sum_{\text{perm}}
 \frac{1}{S_{12}}
 \frac{1}{S_{12}}
 \sum_\tau c^f_\tau~.
\end{eqnarray}
\textbf{Note:} Do not use \verb+\label{#1}+ on a line of a multiline
equation if \verb+\nonumber+ is also used on that line. Incorrect
cross-referencing will result. Notice the use \verb+\text{#1}+ for
using a Roman font within a math environment.

To set a multiline equation without \emph{any} equation
numbers, use the \verb+\begin{eqnarray*}+,
\verb+\end{eqnarray*}+ format:
\begin{eqnarray*}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\\
 & &\times \left( \sum_{i<j}\right)
 \left(
  \sum_{\text{perm}}\frac{1}{S_{12}S_{23}S_{n1}}
 \right)
 \frac{1}{S_{12}}~.
\end{eqnarray*}
To obtain numbers not normally produced by the automatic numbering,
use the \verb+\tag{#1}+ command, where \verb+#1+ is the desired
equation number. For example, to get an equation number of
(\ref{eq:mynum}),
\begin{equation}
g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \tag{2.6$'$}\label{eq:mynum}
\end{equation}

A few notes on \verb=\tag{#1}=. \verb+\tag{#1}+ requires
\texttt{amsmath}. The \verb+\tag{#1}+ must come before the
\verb+\label{#1}+, if any. The numbering set with \verb+\tag{#1}+ is
\textit{transparent} to the automatic numbering in REV\TeX{};
therefore, the number must be known ahead of time, and it must be
manually adjusted if other equations are added. \verb+\tag{#1}+ works
with both single-line and multiline equations. \verb+\tag{#1}+ should
only be used in exceptional case - do not use it to number all
equations in a paper.

Enclosing single-line and multiline equations in
\verb+\begin{subequations}+ and \verb+\end{subequations}+ will produce
a set of equations that are ``numbered'' with letters, as shown in
Eqs.~(\ref{subeq:1}) and (\ref{subeq:2}) below:
\begin{subequations}
\label{eq:whole}
\begin{equation}
\left\{
 abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}
\right\},\label{subeq:1}
\end{equation}
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1).\label{subeq:2}
\end{eqnarray}
\end{subequations}
Putting a \verb+\label{#1}+ command right after the
\verb+\begin{subequations}+, allows one to
reference all the equations in a subequations environment. For
example, the equations in the preceding subequations environment were
Eqs.~(\ref{eq:whole}).

\subsubsection{Wide equations}
The equation that follows is set in a wide format, i.e., it spans
across the full page. The wide format is reserved for long equations
that cannot be easily broken into four lines or less:
\begin{widetext}
\begin{equation}
{\cal R}^{(\text{d})}=
 g_{\sigma_2}^e
 \left(
   \frac{[\Gamma^Z(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^Z(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)
 + x_WQ_e
 \left(
   \frac{[\Gamma^\gamma(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^\gamma(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)\;. \label{eq:wideeq}
\end{equation}
\end{widetext}
This is typed to show the output is in wide format.
(Since there is no input line between \verb+\equation+ and
this paragraph, there is no paragraph indent for this paragraph.)
\section{Cross-referencing}
REV\TeX{} will automatically number sections, equations, figure
captions, and tables. In order to reference them in text, use the
\verb+\label{#1}+ and \verb+\ref{#1}+ commands. To reference a
particular page, use the \verb+\pageref{#1}+ command.

The \verb+\label{#1}+ should appear in a section heading, within an
equation, or in a table or figure caption. The \verb+\ref{#1}+ command
is used in the text where the citation is to be displayed.  Some
examples: Section~\ref{sec:level1} on page~\pageref{sec:level1},
Table~\ref{tab:table1},%
\begin{table}
\caption{\label{tab:table1}This is a narrow table which fits into a
text column when using \texttt{twocolumn} formatting. Note that
REV\TeX~4 adjusts the intercolumn spacing so that the table fills the
entire width of the column. Table captions are numbered
automatically. This table illustrates left-aligned, centered, and
right-aligned columns.  }
\begin{ruledtabular}
\begin{tabular}{lcr}
Left\footnote{Note a.}&Centered\footnote{Note b.}&Right\\
\hline
1 & 2 & 3\\
10 & 20 & 30\\
100 & 200 & 300\\
\end{tabular}
\end{ruledtabular}
\end{table}
and Fig.~\ref{fig:epsart}.

\section{Figures and Tables}
Figures and tables are typically ``floats''; \LaTeX\ determines their
final position via placement rules.
\LaTeX\ isn't always successful in automatically placing floats where you wish them.

Figures are marked up with the \texttt{figure} environment, the content of which
imports the image (\verb+\includegraphics+) followed by the figure caption (\verb+\caption+).
The argument of the latter command should itself contain a \verb+\label+ command if you
wish to refer to your figure with \verb+\ref+.

Import your image using either the \texttt{graphics} or
\texttt{graphix} packages. These packages both define the
\verb+\includegraphics{#1}+ command, but they differ in the optional
arguments for specifying the orientation, scaling, and translation of the figure.
Fig.~\ref{fig:epsart}%
\begin{figure}
\includegraphics{fig_1
\caption{\label{fig:epsart} A figure caption. The figure captions are
automatically numbered.}
\end{figure}
is small enough to fit in a single column, while
Fig.~\ref{fig:wide}%
\begin{figure*}
\includegraphics{fig_2
\caption{\label{fig:wide}Use the \texttt{figure*} environment to get a wide
figure, spanning the page in \texttt{twocolumn} formatting.}
\end{figure*}
is too wide for a single column,
so instead the \texttt{figure*} environment has been used.

The analog of the \texttt{figure} environment is \texttt{table}, which uses
the same \verb+\caption+ command.
However, you should type your caption command first within the \texttt{table},
instead of last as you did for \texttt{figure}.

The heart of any table is the \texttt{tabular} environment,
which represents the table content as a (vertical) sequence of table rows,
each containing a (horizontal) sequence of table cells.
Cells are separated by the \verb+&+ character;
the row terminates with \verb+\\+.
The required argument for the \texttt{tabular} environment
specifies how data are displayed in each of the columns.
For instance, a column
may be centered (\verb+c+), left-justified (\verb+l+), right-justified (\verb+r+),
or aligned on a decimal point (\verb+d+).
(Table~\ref{tab:table4}%
\begin{table}
\caption{\label{tab:table4}Numbers in columns Three--Five have been
aligned by using the ``d'' column specifier (requires the
\texttt{dcolumn} package).
Non-numeric entries (those entries without
a ``.'') in a ``d'' column are aligned on the decimal point.
Use the
``D'' specifier for more complex layouts. }
\begin{ruledtabular}
\begin{tabular}{ccddd}
One&Two&\mbox{Three}&\mbox{Four}&\mbox{Five}\\
\hline
one&two&\mbox{three}&\mbox{four}&\mbox{five}\\
He&2& 2.77234 & 45672. & 0.69 \\
C\footnote{Some tables require footnotes.}
  &C\footnote{Some tables need more than one footnote.}
  & 12537.64 & 37.66345 & 86.37 \\
\end{tabular}
\end{ruledtabular}
\end{table}
illustrates the use of decimal column alignment.)

Extra column-spacing may be be specified as well, although
REV\TeX~4 sets this spacing so that the columns fill the width of the
table.
Horizontal rules are typeset using the \verb+\hline+
command.
The doubled (or Scotch) rules that appear at the top and
bottom of a table can be achieved by enclosing the \texttt{tabular}
environment within a \texttt{ruledtabular} environment.
Rows whose columns span multiple columns can be typeset using \LaTeX's
\verb+\multicolumn{#1}{#2}{#3}+ command
(for example, see the first row of Table~\ref{tab:table3}).%
\begin{table*}
\caption{\label{tab:table3}This is a wide table that spans the page
width in \texttt{twocolumn} mode. It is formatted using the
\texttt{table*} environment. It also demonstrates the use of
\textbackslash\texttt{multicolumn} in rows with entries that span
more than one column.}
\begin{ruledtabular}
\begin{tabular}{ccccc}
 &\multicolumn{2}{c}{$D_{4h}^1$}&\multicolumn{2}{c}{$D_{4h}^5$}\\
 Ion&1st alternative&2nd alternative&lst alternative
&2nd alternative\\ \hline
 K&$(2e)+(2f)$&$(4i)$ &$(2c)+(2d)$&$(4f)$ \\
 Mn&$(2g)$\footnote{The $z$ parameter of these positions is $z\sim\frac{1}{4}$.}
 &$(a)+(b)+(c)+(d)$&$(4e)$&$(2a)+(2b)$\\
 Cl&$(a)+(b)+(c)+(d)$&$(2g)$\footnote{This is a footnote in a table that spans the full page
width in \texttt{twocolumn} mode. It is supposed to set on the full width of the page, just as the caption does. }
 &$(4e)^{\text{a}}$\\
 He&$(8r)^{\text{a}}$&$(4j)^{\text{a}}$&$(4g)^{\text{a}}$\\
 Ag& &$(4k)^{\text{a}}$& &$(4h)^{\text{a}}$\\
\end{tabular}
\end{ruledtabular}
\end{table*}

The tables in this document illustrate various effects.
Tables that fit in a narrow column are contained in a \texttt{table}
environment.
Table~\ref{tab:table3} is a wide table, therefore set with the
\texttt{table*} environment.
Lengthy tables may need to break across pages.
A simple way to allow this is to specify
the \verb+[H]+ float placement on the \texttt{table} or
\texttt{table*} environment.
Alternatively, using the standard \LaTeXe\ package \texttt{longtable}
gives more control over how tables break and allows headers and footers
to be specified for each page of the table.
An example of the use of \texttt{longtable} can be found
in the file \texttt{summary.tex} that is included with the REV\TeX~4
distribution.

There are two methods for setting footnotes within a table (these
footnotes will be displayed directly below the table rather than at
the bottom of the page or in the bibliography).
The easiest
and preferred method is just to use the \verb+\footnote{#1}+
command. This will automatically enumerate the footnotes with
lowercase roman letters.
However, it is sometimes necessary to have
multiple entries in the table share the same footnote.
In this case,
create the footnotes using
\verb+\footnotemark[#1]+ and \verb+\footnotetext[#1]{#2}+.
\texttt{\#1} is a numeric value.
Each time the same value for \texttt{\#1} is used,
the same mark is produced in the table.
The \verb+\footnotetext[#1]{#2}+ commands are placed after the \texttt{tabular}
environment.
Examine the \LaTeX\ source and output for Tables~\ref{tab:table1} and
\ref{tab:table2}%
\begin{table}
\caption{\label{tab:table2}A table with more columns still fits
properly in a column. Note that several entries share the same
footnote. Inspect the \LaTeX\ input for this table to see
exactly how it is done.}
\begin{ruledtabular}
\begin{tabular}{cccccccc}
 &$r_c$ (\AA)&$r_0$ (\AA)&$\kappa r_0$&
 &$r_c$ (\AA) &$r_0$ (\AA)&$\kappa r_0$\\
\hline
Cu& 0.800 & 14.10 & 2.550 &Sn\footnotemark[1]
& 0.680 & 1.870 & 3.700 \\
Ag& 0.990 & 15.90 & 2.710 &Pb\footnotemark[2]
& 0.450 & 1.930 & 3.760 \\
Au& 1.150 & 15.90 & 2.710 &Ca\footnotemark[3]
& 0.750 & 2.170 & 3.560 \\
Mg& 0.490 & 17.60 & 3.200 &Sr\footnotemark[4]
& 0.900 & 2.370 & 3.720 \\
Zn& 0.300 & 15.20 & 2.970 &Li\footnotemark[2]
& 0.380 & 1.730 & 2.830 \\
Cd& 0.530 & 17.10 & 3.160 &Na\footnotemark[5]
& 0.760 & 2.110 & 3.120 \\
Hg& 0.550 & 17.80 & 3.220 &K\footnotemark[5]
&  1.120 & 2.620 & 3.480 \\
Al& 0.230 & 15.80 & 3.240 &Rb\footnotemark[3]
& 1.330 & 2.800 & 3.590 \\
Ga& 0.310 & 16.70 & 3.330 &Cs\footnotemark[4]
& 1.420 & 3.030 & 3.740 \\
In& 0.460 & 18.40 & 3.500 &Ba\footnotemark[5]
& 0.960 & 2.460 & 3.780 \\
Tl& 0.480 & 18.90 & 3.550 & & & & \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{Here's the first, from Ref.~\onlinecite{feyn54}.}
\footnotetext[2]{Here's the second.}
\footnotetext[3]{Here's the third.}
\footnotetext[4]{Here's the fourth.}
\footnotetext[5]{And etc.}
\end{table}
for an illustration.

All AIP journals require that the initial citation of
figures or tables be in numerical order.
\LaTeX's automatic numbering of floats is your friend here:
just put each \texttt{figure} environment immediately following
its first reference (\verb+\ref+), as we have done in this example file.

\begin{acknowledgments}
We wish to acknowledge the support of the author community in using
REV\TeX{}, offering suggestions and encouragement, testing new versions,
\dots.
\end{acknowledgments}

\section{\label{intro} Introduction }

This sample document demonstrates proper use of REV\TeX~4.2 (and
\LaTeXe) in mansucripts prepared for submission to APS
journals. Further information can be found in the REV\TeX~4.2
documentation included in the distribution or available at
\url{http://journals.aps.org/revtex/}.

When commands are referred to in this example file, they are always
shown with their required arguments, using normal \TeX{} format. In
this format, \verb+#1+, \verb+#2+, etc. stand for required
author-supplied arguments to commands. For example, in
\verb+\section{#1}+ the \verb+#1+ stands for the title text of the
author's section heading, and in \verb+\title{#1}+ the \verb+#1+
stands for the title text of the paper.

Line breaks in section headings at all levels can be introduced using
\textbackslash\textbackslash. A blank input line tells \TeX\ that the
paragraph has ended. Note that top-level section headings are
automatically uppercased. If a specific letter or word should appear in
lowercase instead, you must escape it using \verb+\lowercase{#1}+ as
in the word ``via'' above.

\subsection{\label{sec:level2}Second-level heading: Formatting}

This file may be formatted in either the \texttt{preprint} or
\texttt{reprint} style. \texttt{reprint} format mimics final journal output.
Either format may be used for submission purposes. \texttt{letter} sized paper should
be used when submitting to APS journals.

\subsubsection{Wide text (A level-3 head)}
The \texttt{widetext} environment will make the text the width of the
full page, as on page~\pageref{eq:wideeq}. (Note the use the
\verb+\pageref{#1}+ command to refer to the page number.)
\paragraph{Note (Fourth-level head is run in)}
The width-changing commands only take effect in two-column formatting.
There is no effect if text is in a single column.

\subsection{\label{sec:citeref}Citations and References}
A citation in text uses the command \verb+\cite{#1}+ or
\verb+\onlinecite{#1}+ and refers to an entry in the bibliography.
An entry in the bibliography is a reference to another document.

\subsubsection{Citations}
Because REV\TeX\ uses the \verb+natbib+ package of Patrick Daly,
the entire repertoire of commands in that package are available for your document;
see the \verb+natbib+ documentation for further details. Please note that
REV\TeX\ requires version 8.31a or later of \verb+natbib+.

\paragraph{Syntax}
The argument of \verb+\cite+ may be a single \emph{key},
or may consist of a comma-separated list of keys.
The citation \emph{key} may contain
letters, numbers, the dash (-) character, or the period (.) character.
New with natbib 8.3 is an extension to the syntax that allows for
a star (*) form and two optional arguments on the citation key itself.
The syntax of the \verb+\cite+ command is thus (informally stated)
\begin{quotation}\flushleft\leftskip1em
\verb+\cite+ \verb+{+ \emph{key} \verb+}+, or\\
\verb+\cite+ \verb+{+ \emph{optarg+key} \verb+}+, or\\
\verb+\cite+ \verb+{+ \emph{optarg+key} \verb+,+ \emph{optarg+key}\ldots \verb+}+,
\end{quotation}\noindent
where \emph{optarg+key} signifies
\begin{quotation}\flushleft\leftskip1em
\emph{key}, or\\
\texttt{*}\emph{key}, or\\
\texttt{[}\emph{pre}\texttt{]}\emph{key}, or\\
\texttt{[}\emph{pre}\texttt{]}\texttt{[}\emph{post}\texttt{]}\emph{key}, or even\\
\texttt{*}\texttt{[}\emph{pre}\texttt{]}\texttt{[}\emph{post}\texttt{]}\emph{key}.
\end{quotation}\noindent
where \emph{pre} and \emph{post} is whatever text you wish to place
at the beginning and end, respectively, of the bibliographic reference
(see Ref.~[\onlinecite{witten2001}] and the two under Ref.~[\onlinecite{feyn54}]).
(Keep in mind that no automatic space or punctuation is applied.)
It is highly recommended that you put the entire \emph{pre} or \emph{post} portion
within its own set of braces, for example:
\verb+\cite+ \verb+{+ \texttt{[} \verb+{+\emph{text}\verb+}+\texttt{]}\emph{key}\verb+}+.
The extra set of braces will keep \LaTeX\ out of trouble if your \emph{text} contains the comma (,) character.

The star (*) modifier to the \emph{key} signifies that the reference is to be
merged with the previous reference into a single bibliographic entry,
a common idiom in APS and AIP articles (see below, Ref.~[\onlinecite{epr}]).
When references are merged in this way, they are separated by a semicolon instead of
the period (full stop) that would otherwise appear.

\paragraph{Eliding repeated information}
When a reference is merged, some of its fields may be elided: for example,
when the author matches that of the previous reference, it is omitted.
If both author and journal match, both are omitted.
If the journal matches, but the author does not, the journal is replaced by \emph{ibid.},
as exemplified by Ref.~[\onlinecite{epr}].
These rules embody common editorial practice in APS and AIP journals and will only
be in effect if the markup features of the APS and AIP Bib\TeX\ styles is employed.

\paragraph{The options of the cite command itself}
Please note that optional arguments to the \emph{key} change the reference in the bibliography,
not the citation in the body of the document.
For the latter, use the optional arguments of the \verb+\cite+ command itself:
\verb+\cite+ \texttt{*}\allowbreak
\texttt{[}\emph{pre-cite}\texttt{]}\allowbreak
\texttt{[}\emph{post-cite}\texttt{]}\allowbreak
\verb+{+\emph{key-list}\verb+}+.

\subsubsection{Example citations}
By default, citations are numerical\cite{Beutler1994}.
Author-year citations are used when the journal is RMP.
To give a textual citation, use \verb+\onlinecite{#1}+:
Refs.~\onlinecite{[][{, and references therein}]witten2001,Bire82}.
By default, the \texttt{natbib} package automatically sorts your citations into numerical order and ``compresses'' runs of three or more consecutive numerical citations.
REV\TeX\ provides the ability to automatically change the punctuation when switching between journal styles that provide citations in square brackets and those that use a superscript style instead. This is done through the \texttt{citeautoscript} option. For instance, the journal style \texttt{prb} automatically invokes this option because \textit{Physical
Review B} uses superscript-style citations. The effect is to move the punctuation, which normally comes after a citation in square brackets, to its proper position before the superscript.
To illustrate, we cite several together
\cite{[See the explanation of time travel in ]feyn54,*[The classical relativistic treatment of ][ is a relative classic]epr,witten2001,Berman1983,Davies1998,Bire82},
and once again in different order (Refs.~\cite{epr,feyn54,Bire82,Berman1983,witten2001,Davies1998}).
Note that the citations were both compressed and sorted. Futhermore, running this sample file under the \texttt{prb} option will move the punctuation to the correct place.

When the \verb+prb+ class option is used, the \verb+\cite{#1}+ command
displays the reference's number as a superscript rather than in
square brackets. Note that the location of the \verb+\cite{#1}+
command should be adjusted for the reference style: the superscript
references in \verb+prb+ style must appear after punctuation;
otherwise the reference must appear before any punctuation. This
sample was written for the regular (non-\texttt{prb}) citation style.
The command \verb+\onlinecite{#1}+ in the \texttt{prb} style also
displays the reference on the baseline.

\subsubsection{References}
A reference in the bibliography is specified by a \verb+\bibitem{#1}+ command
with the same argument as the \verb+\cite{#1}+ command.
\verb+\bibitem{#1}+ commands may be crafted by hand or, preferably,
generated by Bib\TeX.
REV\TeX~4.2 includes Bib\TeX\ style files
\verb+apsrev4-2.bst+, \verb+apsrmp4-2.bst+ appropriate for
\textit{Physical Review} and \textit{Reviews of Modern Physics},
respectively.

\subsubsection{Example references}
This sample file employs the \verb+\bibliography+ command,
which formats the \texttt{\jobname .bbl} file
and specifies which bibliographic databases are to be used by Bib\TeX\
(one of these should be by arXiv convention \texttt{\jobname .bib}).
Running Bib\TeX\ (via \texttt{bibtex \jobname})
after the first pass of \LaTeX\ produces the file
\texttt{\jobname .bbl} which contains the automatically formatted
\verb+\bibitem+ commands (including extra markup information via
\verb+\bibinfo+ and \verb+\bibfield+ commands).
If not using Bib\TeX, you will have to create the \verb+thebibiliography+ environment
and its \verb+\bibitem+ commands by hand.

Numerous examples of the use of the APS bibliographic entry types appear in the bibliography of this sample document.
You can refer to the \texttt{\jobname .bib} file,
and compare its information to the formatted bibliography itself.

\subsection{Footnotes}%
Footnotes, produced using the \verb+\footnote{#1}+ command,
usually integrated into the bibliography alongside the other entries.
Numerical citation styles do this%
\footnote{Automatically placing footnotes into the bibliography requires using BibTeX to compile the bibliography.};
author-year citation styles place the footnote at the bottom of the text column.
Note: due to the method used to place footnotes in the bibliography,
\emph{you must re-run Bib\TeX\ every time you change any of your document's footnotes}.

\section{Math and Equations}
Inline math may be typeset using the \verb+$+ delimiters. Bold math
symbols may be achieved using the \verb+bm+ package and the
\verb+\bm{#1}+ command it supplies. For instance, a bold $\alpha$ can
be typeset as \verb+$\bm{\alpha}$+ giving $\bm{\alpha}$. Fraktur and
Blackboard (or open face or double struck) characters should be
typeset using the \verb+\mathfrak{#1}+ and \verb+\mathbb{#1}+ commands
respectively. Both are supplied by the \texttt{amssymb} package. For
example, \verb+$\mathbb{R}$+ gives $\mathbb{R}$ and
\verb+$\mathfrak{G}$+ gives $\mathfrak{G}$

In \LaTeX\ there are many different ways to display equations, and a
few preferred ways are noted below. Displayed math will center by
default. Use the class option \verb+fleqn+ to flush equations left.

Below we have numbered single-line equations; this is the most common
type of equation in \textit{Physical Review}:
\begin{eqnarray}
\chi_+(p)\alt{\bf [}2|{\bf p}|(|{\bf p}|+p_z){\bf ]}^{-1/2}
\left(
\begin{array}{c}
|{\bf p}|+p_z\\
px+ip_y
\end{array}\right)\;,
\\
\left\{%
 \openone234567890abc123\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}%
\label{eq:one}.
\end{eqnarray}
Note the open one in Eq.~(\ref{eq:one}).

Not all numbered equations will fit within a narrow column this
way. The equation number will move down automatically if it cannot fit
on the same line with a one-line equation:
\begin{equation}
\left\{
 ab12345678abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}%
\right\}.
\end{equation}

When the \verb+\label{#1}+ command is used [cf. input for
Eq.~(\ref{eq:one})], the equation can be referred to in text without
knowing the equation number that \TeX\ will assign to it. Just
use \verb+\ref{#1}+, where \verb+#1+ is the same name that used in
the \verb+\label{#1}+ command.

Unnumbered single-line equations can be typeset
using the \verb+\[+, \verb+\]+ format:
\[g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \]

\subsection{Multiline equations}

Multiline equations are obtained by using the \verb+eqnarray+
environment.  Use the \verb+\nonumber+ command at the end of each line
to avoid assigning a number:
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
\delta_{\sigma_1,-\sigma_2}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_jl_i\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1),
\end{eqnarray}
\begin{eqnarray}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\nonumber \\
 & &\times \left( \sum_{i<j}\right)
  \sum_{\text{perm}}
 \frac{1}{S_{12}}
 \frac{1}{S_{12}}
 \sum_\tau c^f_\tau~.
\end{eqnarray}
\textbf{Note:} Do not use \verb+\label{#1}+ on a line of a multiline
equation if \verb+\nonumber+ is also used on that line. Incorrect
cross-referencing will result. Notice the use \verb+\text{#1}+ for
using a Roman font within a math environment.

To set a multiline equation without \emph{any} equation
numbers, use the \verb+\begin{eqnarray*}+,
\verb+\end{eqnarray*}+ format:
\begin{eqnarray*}
\sum \vert M^{\text{viol}}_g \vert ^2&=&g^{2n-4}_S(Q^2)~N^{n-2}
        (N^2-1)\\
 & &\times \left( \sum_{i<j}\right)
 \left(
  \sum_{\text{perm}}\frac{1}{S_{12}S_{23}S_{n1}}
 \right)
 \frac{1}{S_{12}}~.
\end{eqnarray*}

To obtain numbers not normally produced by the automatic numbering,
use the \verb+\tag{#1}+ command, where \verb+#1+ is the desired
equation number. For example, to get an equation number of
(\ref{eq:mynum}),
\begin{equation}
g^+g^+ \rightarrow g^+g^+g^+g^+ \dots ~,~~q^+q^+\rightarrow
q^+g^+g^+ \dots ~. \tag{2.6$'$}\label{eq:mynum}
\end{equation}

\paragraph{A few notes on \texttt{tag}s}
\verb+\tag{#1}+ requires the \texttt{amsmath} package.
Place the \verb+\tag{#1}+ command before the \verb+\label{#1}+, if any.
The numbering produced by \verb+\tag{#1}+ \textit{does not affect}
the automatic numbering in REV\TeX;
therefore, the number must be known ahead of time,
and it must be manually adjusted if other equations are added.
\verb+\tag{#1}+ works with both single-line and multiline equations.
\verb+\tag{#1}+ should only be used in exceptional cases---%
do not use it to number many equations in your paper.
Please note that this feature of the \texttt{amsmath} package
is \emph{not} compatible with the \texttt{hyperref} (6.77u) package.

Enclosing display math within
\verb+\begin{subequations}+ and \verb+\end{subequations}+ will produce
a set of equations that are labeled with letters, as shown in
Eqs.~(\ref{subeq:1}) and (\ref{subeq:2}) below.
You may include any number of single-line and multiline equations,
although it is probably not a good idea to follow one display math
directly after another.
\begin{subequations}
\label{eq:whole}
\begin{eqnarray}
{\cal M}=&&ig_Z^2(4E_1E_2)^{1/2}(l_i^2)^{-1}
(g_{\sigma_2}^e)^2\chi_{-\sigma_2}(p_2)\nonumber\\
&&\times
[\epsilon_i]_{\sigma_1}\chi_{\sigma_1}(p_1).\label{subeq:2}
\end{eqnarray}
\begin{equation}
\left\{
 abc123456abcdef\alpha\beta\gamma\delta1234556\alpha\beta
 \frac{1\sum^{a}_{b}}{A^2}
\right\},\label{subeq:1}
\end{equation}
\end{subequations}
Giving a \verb+\label{#1}+ command directly after the \verb+\begin{subequations}+,
allows you to reference all the equations in the \texttt{subequations} environment.
For example, the equations in the preceding subequations environment were
Eqs.~(\ref{eq:whole}).

\subsubsection{Wide equations}
The equation that follows is set in a wide format, i.e., it spans the full page.
The wide format is reserved for long equations
that cannot easily be set in a single column:
\begin{widetext}
\begin{equation}
{\cal R}^{(\text{d})}=
 g_{\sigma_2}^e
 \left(
   \frac{[\Gamma^Z(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^Z(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)
 + x_WQ_e
 \left(
   \frac{[\Gamma^\gamma(3,21)]_{\sigma_1}}{Q_{12}^2-M_W^2}
  +\frac{[\Gamma^\gamma(13,2)]_{\sigma_1}}{Q_{13}^2-M_W^2}
 \right)\;.
 \label{eq:wideeq}
\end{equation}
\end{widetext}
This is typed to show how the output appears in wide format.
(Incidentally, since there is no blank line between the \texttt{equation} environment above
and the start of this paragraph, this paragraph is not indented.)

\section{Cross-referencing}
REV\TeX{} will automatically number such things as
sections, footnotes, equations, figure captions, and table captions.
In order to reference them in text, use the
\verb+\label{#1}+ and \verb+\ref{#1}+ commands.
To reference a particular page, use the \verb+\pageref{#1}+ command.

The \verb+\label{#1}+ should appear
within the section heading,
within the footnote text,
within the equation, or
within the table or figure caption.
The \verb+\ref{#1}+ command
is used in text at the point where the reference is to be displayed.
Some examples: Section~\ref{sec:level1} on page~\pageref{sec:level1},
Table~\ref{tab:table1},%
\begin{table}[b
\caption{\label{tab:table1}%
A table that fits into a single column of a two-column layout.
Note that REV\TeX~4 adjusts the intercolumn spacing so that the table fills the
entire width of the column. Table captions are numbered
automatically.
This table illustrates left-, center-, decimal- and right-aligned columns,
along with the use of the \texttt{ruledtabular} environment which sets the
Scotch (double) rules above and below the alignment, per APS style.
}
\begin{ruledtabular}
\begin{tabular}{lcdr}
\textrm{Left\footnote{Note a.}}&
\textrm{Centered\footnote{Note b.}}&
\multicolumn{1}{c}{\textrm{Decimal}}&
\textrm{Right}\\
\colrule
1 & 2 & 3.001 & 4\\
10 & 20 & 30 & 40\\
100 & 200 & 300.0 & 400\\
\end{tabular}
\end{ruledtabular}
\end{table}
and Fig.~\ref{fig:epsart}.%
\begin{figure}[b]
\includegraphics{fig_1
\caption{\label{fig:epsart} A figure caption. The figure captions are
automatically numbered.}
\end{figure}

\section{Floats: Figures, Tables, Videos, etc.}
Figures and tables are usually allowed to ``float'', which means that their
placement is determined by \LaTeX, while the document is being typeset.

Use the \texttt{figure} environment for a figure, the \texttt{table} environment for a table.
In each case, use the \verb+\caption+ command within to give the text of the
figure or table caption along with the \verb+\label+ command to provide
a key for referring to this figure or table.
The typical content of a figure is an image of some kind;
that of a table is an alignment.%
\begin{figure*}
\includegraphics{fig_2
\caption{\label{fig:wide}Use the figure* environment to get a wide
figure that spans the page in \texttt{twocolumn} formatting.}
\end{figure*}
\begin{table*}
\caption{\label{tab:table3}This is a wide table that spans the full page
width in a two-column layout. It is formatted using the
\texttt{table*} environment. It also demonstates the use of
\textbackslash\texttt{multicolumn} in rows with entries that span
more than one column.}
\begin{ruledtabular}
\begin{tabular}{ccccc}
 &\multicolumn{2}{c}{$D_{4h}^1$}&\multicolumn{2}{c}{$D_{4h}^5$}\\
 Ion&1st alternative&2nd alternative&lst alternative
&2nd alternative\\ \hline
 K&$(2e)+(2f)$&$(4i)$ &$(2c)+(2d)$&$(4f)$ \\
 Mn&$(2g)$\footnote{The $z$ parameter of these positions is $z\sim\frac{1}{4}$.}
 &$(a)+(b)+(c)+(d)$&$(4e)$&$(2a)+(2b)$\\
 Cl&$(a)+(b)+(c)+(d)$&$(2g)$\footnotemark[1]
 &$(4e)^{\text{a}}$\\
 He&$(8r)^{\text{a}}$&$(4j)^{\text{a}}$&$(4g)^{\text{a}}$\\
 Ag& &$(4k)^{\text{a}}$& &$(4h)^{\text{a}}$\\
\end{tabular}
\end{ruledtabular}
\end{table*}

Insert an image using either the \texttt{graphics} or
\texttt{graphix} packages, which define the \verb+\includegraphics{#1}+ command.
(The two packages differ in respect of the optional arguments
used to specify the orientation, scaling, and translation of the image.)
To create an alignment, use the \texttt{tabular} environment.

The best place to locate the \texttt{figure} or \texttt{table} environment
is immediately following its first reference in text; this sample document
illustrates this practice for Fig.~\ref{fig:epsart}, which
shows a figure that is small enough to fit in a single column.

In exceptional cases, you will need to move the float earlier in the document, as was done
with Table~\ref{tab:table3}: \LaTeX's float placement algorithms need to know
about a full-page-width float earlier.

Fig.~\ref{fig:wide}
has content that is too wide for a single column,
so the \texttt{figure*} environment has been used.%
\begin{table}[b]
\caption{\label{tab:table4}%
Numbers in columns Three--Five are aligned with the ``d'' column specifier
(requires the \texttt{dcolumn} package).
Non-numeric entries (those entries without a ``.'') in a ``d'' column are aligned on the decimal point.
Use the ``D'' specifier for more complex layouts. }
\begin{ruledtabular}
\begin{tabular}{ccddd}
One&Two&
\multicolumn{1}{c}{\textrm{Three}}&
\multicolumn{1}{c}{\textrm{Four}}&
\multicolumn{1}{c}{\textrm{Five}}\\
\hline
one&two&\mbox{three}&\mbox{four}&\mbox{five}\\
He&2& 2.77234 & 45672. & 0.69 \\
C\footnote{Some tables require footnotes.}
  &C\footnote{Some tables need more than one footnote.}
  & 12537.64 & 37.66345 & 86.37 \\
\end{tabular}
\end{ruledtabular}
\end{table}

The content of a table is typically a \texttt{tabular} environment,
giving rows of type in aligned columns.
Column entries separated by \verb+&+'s, and
each row ends with \textbackslash\textbackslash.
The required argument for the \texttt{tabular} environment
specifies how data are aligned in the columns.
For instance, entries may be centered, left-justified, right-justified, aligned on a decimal
point.
Extra column-spacing may be be specified as well,
although REV\TeX~4 sets this spacing so that the columns fill the width of the
table. Horizontal rules are typeset using the \verb+\hline+
command. The doubled (or Scotch) rules that appear at the top and
bottom of a table can be achieved enclosing the \texttt{tabular}
environment within a \texttt{ruledtabular} environment. Rows whose
columns span multiple columns can be typeset using the
\verb+\multicolumn{#1}{#2}{#3}+ command (for example, see the first
row of Table~\ref{tab:table3}).%

Tables~\ref{tab:table1}, \ref{tab:table3}, \ref{tab:table4}, and \ref{tab:table2}%
\begin{table}[b]
\caption{\label{tab:table2}
A table with numerous columns that still fits into a single column.
Here, several entries share the same footnote.
Inspect the \LaTeX\ input for this table to see exactly how it is done.}
\begin{ruledtabular}
\begin{tabular}{cccccccc}
 &$r_c$ (\AA)&$r_0$ (\AA)&$\kappa r_0$&
 &$r_c$ (\AA) &$r_0$ (\AA)&$\kappa r_0$\\
\hline
Cu& 0.800 & 14.10 & 2.550 &Sn\footnotemark[1]
& 0.680 & 1.870 & 3.700 \\
Ag& 0.990 & 15.90 & 2.710 &Pb\footnotemark[2]
& 0.450 & 1.930 & 3.760 \\
Au& 1.150 & 15.90 & 2.710 &Ca\footnotemark[3]
& 0.750 & 2.170 & 3.560 \\
Mg& 0.490 & 17.60 & 3.200 &Sr\footnotemark[4]
& 0.900 & 2.370 & 3.720 \\
Zn& 0.300 & 15.20 & 2.970 &Li\footnotemark[2]
& 0.380 & 1.730 & 2.830 \\
Cd& 0.530 & 17.10 & 3.160 &Na\footnotemark[5]
& 0.760 & 2.110 & 3.120 \\
Hg& 0.550 & 17.80 & 3.220 &K\footnotemark[5]
&  1.120 & 2.620 & 3.480 \\
Al& 0.230 & 15.80 & 3.240 &Rb\footnotemark[3]
& 1.330 & 2.800 & 3.590 \\
Ga& 0.310 & 16.70 & 3.330 &Cs\footnotemark[4]
& 1.420 & 3.030 & 3.740 \\
In& 0.460 & 18.40 & 3.500 &Ba\footnotemark[5]
& 0.960 & 2.460 & 3.780 \\
Tl& 0.480 & 18.90 & 3.550 & & & & \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{Here's the first, from Ref.~\onlinecite{feyn54}.}
\footnotetext[2]{Here's the second.}
\footnotetext[3]{Here's the third.}
\footnotetext[4]{Here's the fourth.}
\footnotetext[5]{And etc.}
\end{table}
show various effects.
A table that fits in a single column employs the \texttt{table}
environment.
Table~\ref{tab:table3} is a wide table, set with the \texttt{table*} environment.
Long tables may need to break across pages.
The most straightforward way to accomplish this is to specify
the \verb+[H]+ float placement on the \texttt{table} or
\texttt{table*} environment.
However, the \LaTeXe\ package \texttt{longtable} allows headers and footers to be specified for each page of the table.
A simple example of the use of \texttt{longtable} can be found
in the file \texttt{summary.tex} that is included with the REV\TeX~4
distribution.

There are two methods for setting footnotes within a table (these
footnotes will be displayed directly below the table rather than at
the bottom of the page or in the bibliography). The easiest
and preferred method is just to use the \verb+\footnote{#1}+
command. This will automatically enumerate the footnotes with
lowercase roman letters. However, it is sometimes necessary to have
multiple entries in the table share the same footnote. In this case,
there is no choice but to manually create the footnotes using
\verb+\footnotemark[#1]+ and \verb+\footnotetext[#1]{#2}+.
\texttt{\#1} is a numeric value. Each time the same value for
\texttt{\#1} is used, the same mark is produced in the table. The
\verb+\footnotetext[#1]{#2}+ commands are placed after the \texttt{tabular}
environment. Examine the \LaTeX\ source and output for
Tables~\ref{tab:table1} and \ref{tab:table2}
for examples.

Video~\ref{vid:PRSTPER.4.010101}
illustrates several features new with REV\TeX4.2,
starting with the \texttt{video} environment, which is in the same category with
\texttt{figure} and \texttt{table}.%
\begin{video}
\href{http://prst-per.aps.org/multimedia/PRSTPER/v4/i1/e010101/e010101_vid1a.mpg}{\includegraphics{vid_1a}}%
 \quad
\href{http://prst-per.aps.org/multimedia/PRSTPER/v4/i1/e010101/e010101_vid1b.mpg}{\includegraphics{vid_1b}}
 \setfloatlink{http://link.aps.org/multimedia/PRSTPER/v4/i1/e010101}%
 \caption{\label{vid:PRSTPER.4.010101}%
  Students explain their initial idea about Newton's third law to a teaching assistant.
  Clip (a): same force.
  Clip (b): move backwards.
 }%
\end{video}
The \verb+\setfloatlink+ command causes the title of the video to be a hyperlink to the
indicated URL; it may be used with any environment that takes the \verb+\caption+
command.
The \verb+\href+ command has the same significance as it does in the context of
the \texttt{hyperref} package: the second argument is a piece of text to be
typeset in your document; the first is its hyperlink, a URL.

\textit{Physical Review} style requires that the initial citation of
figures or tables be in numerical order in text, so don't cite
Fig.~\ref{fig:wide} until Fig.~\ref{fig:epsart} has been cited.

\begin{acknowledgments}
We wish to acknowledge the support of the author community in using
REV\TeX{}, offering suggestions and encouragement, testing new versions,
\dots.
\end{acknowledgments}

\section*{Acknowledgments}

The authors are thankful to L\'eonie Canet and Andrea Zoia for insightful discussions and comments. \\

\section{ A Field theoretic formulation of neutron physics}
\label{sec:FieldFormulation}
{ \color{red}
Branching processes are among the simplest models capable of characterizing the complex phenomenology of the evolution of the neutron population in fissile materials (REFS). These processes are also known to display a second order phase transition between an active and an inactive state (REFS). This has a direct and practical transcription in the realm of neutron physics, for the phase transition translate into the formation of self-sustained fission chain reactions or  criticality, at the core root of nuclear reactors operations. The exponential increase or decrease of the neutron population corresponding to the conditions of the apparitions of self-sustained chains is controlled by a single parameter, the reactivity $\rho$ of the medium.

Although branching processes are an invaluable tool to study the criticality of a nuclear reactor, this basic modelling is of course insufficient to properly encompass the complexity and variety of phenomena occurring in a nuclear reactor, especially when the latter is operating at a significant power level. A proficient tool for the analysis of nuclear reactor stability is to superimpose on the branching process a stochastic noise coming from the different noise sources that can alter the initial model. In other words, one suppose that the overall neutron population can be decoupled into
\begin{itemize}
    \item a smoothly varying neutron field $N(t)$. It correspond to the measured field observed in a reactor (through the operating power level for instance) and is only sensitive to the global reactivity $\rho$ of the reactor;

    \item a random source of noise $\eta(t)$, source of perturbations and driven by unknown or at least unspecified phenomena.
\end{itemize}
}
The random noise source term can be added to the mean field equation derived from the branching process. The latter then is transformed into a Stochastic Differential Equation (SDE),
\begin{equation}
\label{LangevinEquation}
    \cfrac{d N(t)}{dt} = \rho N(t) + \eta(t).
\end{equation}

Any observable quantity $\mathcal{O}[N]$ must be averaged over all possibles configurations for the noise term. This is formally achieved by performing a functional integration over the noise probability functional $P[\eta]$:
\begin{equation}
\label{observable}
    \langle \mathcal{O}[N] \rangle \, \propto \, \int \mathcal{D}\eta \, \mathcal{O}[N] \, P[\eta].
\end{equation}

Building on the seminal work of Janssen \cite{janssen1976} and De Dominicis \cite{dominicis1976}, one can take advantage of the structure of the problem to develop a field-theoretic formulation of the model, simplifying greatly the calculation of observables.

\subsection{The response functional formalism}
\label{subsec:ResponseFunctional}

The development of the field-theoretic version of the present model starts by noting that the $N$ appearing in equation \ref{observable} is constrained to be a solution of equation \ref{LangevinEquation}. The constraint (valid for all $t$),
\begin{equation}
    \label{constraint_N}
    \mathcal{C}[N] \, = \, \left\{\frac{dN(t)}{dt} \, - \, \rho \, N(t) \,  -  \, n(t)\right\} = 0,
\end{equation}
can be enforced with the help of a functional version of the resolution of the identity \cite{cardy_tauber},
\begin{equation}
    1 = \int \mathcal{D}N \, \prod_t  \delta \left( \mathcal{C}[N] \right).
\end{equation}
Performing a (functional) Fourier transformation, one can write,
\begin{eqnarray}
\label{constrained_field}
    1 & = & \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \, e^{- \int dt \ensuremath{\widetilde{N}}(t) \, \mathcal{C}[N]} \nonumber\\
    & = & \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \, e^{- \int dt \ensuremath{\widetilde{N}}(t) \, \left( \frac{dN(t)}{dt} \, - \, \rho \, N(t) \,  -  \, \eta(t)\right)},
\end{eqnarray}
where the purely imaginary auxiliary field $\ensuremath{\widetilde{N}}(t)$, known as the Martin-Siggia-Rose response field \cite{martin_siggia_rose} has been introduced. Inserting this Fourier transformed identity back into equation \ref{observable}, one obtains,
\begin{eqnarray}
\label{observable2}
    \langle \mathcal{O}[N] \rangle & \propto & \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \left\{ e^{- \int \ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt}  -  \rho   \right)N(t) \, dt} \right. \nonumber \\
    & \times &\mathcal{O}[N]
   \left.  \times  \int \mathcal{D}\eta \, e^{\int \ensuremath{\widetilde{N}}(t) \eta(t) dt} P[\eta] \right\}.
\end{eqnarray}
To go further, one needs to integrate over $\eta$ and thus specify the probability functional to which it is associated. In the following, it will be supposed that at any given moment in time, the noise possesses a Gaussian probability distribution
\begin{equation}
\label{noise_proba}
    P[\eta] \propto e^{ -\frac{1}{4}\int \eta(t) \Gamma^{-1} \eta(t) \, dt }
\end{equation}
with zero mean $\langle \eta(t) \rangle = 0$ and a variance given by,
\begin{equation}
    \langle n(t) n(t') \rangle = 2 \Gamma(N) \delta(t - t').
\end{equation}
Note that an explicit dependence of the correlation function $\Gamma(N)$ to the global power level (i.e. $N(t)$) is to be expected on physical ground : the expected noise in a reactor indeed is dependent on the operating power level.

Imposing a Gaussian profile has the merit of making the integral over the noise in equation \ref{observable2} tractable. One has,
\begin{equation}
\int \mathcal{D}n \, e^{\int \ensuremath{\widetilde{N}}(t) n(t) dt} P[n] \propto e^{ \Gamma \ensuremath{\widetilde{N}}^2}.
\end{equation}
Inserting this results into equation \ref{observable2}, one finally arrives at an expression that solely depends on the variance of the noise probability distribution:
\begin{equation}
\label{observable_final}
    \langle \mathcal{O}[N] \rangle \, = \mathcal{N}  \,  \int \mathcal{D}[i\ensuremath{\widetilde{N}}] \int \mathcal{D}N \, \mathcal{O}[N] \,
    e^{- \mathcal{S}[N,\ensuremath{\widetilde{N}}]},
\end{equation}
where both a a normalisation factor $\mathcal{N}$ and the so-called \textit{response functional} have been introduced,
\begin{equation}
\label{response_functional}
     \mathcal{S}[N,\ensuremath{\widetilde{N}}] = \int \left\{\ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt} - \rho \right)N(t) - \ensuremath{\widetilde{N}}(t) \Gamma \ensuremath{\widetilde{N}}(t)\, \right\} dt.
\end{equation}
One could even go further and add another constraint on the model of equation \ref{constraint_N} to impose an initial condition for the field $N(t)$. Imposing that $N(t_0) = N_0$, one can show that the response functional becomes,
\begin{eqnarray}
\label{response_functional_initial_cond}
     \mathcal{S}[N,\ensuremath{\widetilde{N}}]  & \, & = \int \left\{  \ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt} - \rho \right)N(t)          \right.\nonumber\\
      && \left. -  \ensuremath{\widetilde{N}}(t) \, \delta(t-t_0) \, N_0   -   \ensuremath{\widetilde{N}}(t) \Gamma \ensuremath{\widetilde{N}}(t)\, \right\} dt.
\end{eqnarray}

\section{Addition of counter-reactions - Directed percolation }
\label{sec:directed_percolation}

Up until now, the response functional formalism was used to retrieve most of the main results well known in nuclear reactor physics. Starting from the mean field equations derived from an underlying microscopic branching process, the proposed approach was capable of retrieving with relative ease the main results one would usually derive from solving the associated mean field equations.

But because it is built on a macroscopic level of description, the principal advantage of the present field-theoretic approach has not been exploited yet: the formalism can easily incorporate effects for which a microscopic description is not available. In reactor physics, this is for instance the case of the feedback mechanisms applying on the neutron population in a reactor operating a high power (Doppler effect, Xenon poisoning, etc.). These phenomena are a central aspect of the stability of nuclear reactors but their detailed modelling is generally difficult to implement in a microscopic neutron balance.

Because the basic object of the field-theoretic approach proposed in this paper is already the coarse grain field $N(t)$, the introduction of, at least, a phenomenological model is much easier. Although now a direct link between the mean field equations and an underlying stochastic process is lost, one still is capable of extracting information about the moments associated to the neutron population distribution.

One of the most paradigmatic feedback effect affecting the neutron population, and a key element in the study of reactors stability, is the Doppler broadening effect associated to the capture resonances of $^{238}$U. This effect originates from an increase in temperature observed in fissile materials when the neutron population increases. This induces a kinematic broadening of all of the resonant reactions in the medium, in particular the resonant capture of $^{238}$U (in low enriched reactors). As a consequence, when the temperature increases, more and more neutrons ends up being captured by the medium and the overall population is pulled back towards the reference operating power. The same mechanism applies in reverse when the temperature decreases.

 This phenomenon is quite difficult to include in the microscopic branching process followed by the neutron population but can be quite straightforwardly implemented under the form of a restoring force of the type $\sigma N(N - N_{\text{ref.}})$ in the diffusion equation \footnote{under the assumption that the characteristic time associated to the Doppler broadening process is much faster than the mean generation time of the neutron population.} \cite{bell_glasstone}. The resulting equation can thus be written,
\begin{equation}
    \label{stoch_diffusion_doppler}
    \frac{\partial }{ \partial t} N(\vec{x},t) - D ( \nabla^2 + r)  N(\vec{x},t)  +  \sigma N^2(\vec{x},t)  = 0,
\end{equation}
where the shifted parameter $ r = (\rho + \sigma N_{\text{ref.}})/D$ has been introduced, for convenience.

The response functional and the resulting theory to which this equation is associated has been extensively studied. It corresponds to the functional of the \textit{directed percolation} universality class \cite{janssen2004}, \cite{tauber_2005}:
\begin{equation}
    \mathcal{S} = \int d^d \vec{x} \int dt \ensuremath{\widetilde{N}} \left( \partial_t - D ( \vec{\nabla}^2 + r)  \right) N + \sigma \ensuremath{\widetilde{N}} N^2 - \lambda \ensuremath{\widetilde{N}}^2 N  .
\end{equation}

The inclusion of the term $\sigma N^2$ has a strong impact on the structure of the theory, for now the two following processes are allowed,
\begin{equation*}
    2 \lambda \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of a](o1);
    \diagram* {
    (f1) --[anti fermion] (a) --[fermion] (f2),
    (a) --[anti fermion] (o1)
    };
    \end{feynman}
  \end{tikzpicture} \hspace{2 cm} 2 \sigma
      \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [above right=of a](f1);
      \vertex [below right=of a](f2);
      \vertex [left=of a](o1);
    \diagram* {
    (f1) --[fermion] (a) --[anti fermion] (f2),
    (a) --[fermion] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
\end{equation*}
as well as their combination, i.e. the loop diagram,
\begin{equation}
    4 \sigma \lambda \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (z);
      \vertex [right=of z](a);
      \vertex [right=of a] (b);
      \vertex [above left=of a](f1);
      \vertex [below left=of a](f2);
      \vertex [right=of b](o1);
    \diagram* {
     (z) --[anti fermion] (a)
     -- [fermion,half left] (b)
     -- [fermion,half left] (a),
    (b) --[anti fermion] (o1)
    };
    \end{feynman}
  \end{tikzpicture} \, ,
\end{equation}
whose contribution leads to badly divergent behaviours of two different types \cite{tauber_2005}:
\begin{itemize}
    \item a divergence is observed in dimension $d > 2$, because the momentum  associated to the loop is not constrained by any conservation rule. It's integration thus runs over theoretically infinite momenta and therefore include small distances, short time unbounded fluctuations of the neutron population. This is a paradigmatic example of an Ultra-Violet (UV) divergence.

    \item another and more serious divergence occurs in dimension $d <4$. An Infrared (IR) divergence is observed in the long time, long distance limit. The highest dimension $d_c = 4$ for which the divergence is observed is called the \textit{critical dimension} of the problem.
\end{itemize}

Both types of divergences can be eliminated through the well-known renormalization procedure \cite{Delamotte_2004}, which, in essence, consists in the application of successive spatio-temporal scaling transformation on the problem, so that the small distance and UV divergent fluctuations be truly integrated over. The true macro-scale description of the problem is obtained at a fixed point of the iterated scaling process which generally has the induced benefit of taming the residual IR divergence of the problem. The whole renormalization procedure leaves the system with two fundamental characteristics:
\begin{itemize}
    \item The asymptotic behaviour of the theory does not depend on the fine grain details and, in particular, on the precise origin of the phenomena at play. The resulting theory is \textit{universal}.

    \item Because the theory resulting from the renormalization procedure is derived at a fixed point of a spatio-temporal scaling transformation, the observables extracted from it will typically present power law behaviours. An imprint of the system's fluctuations will generally subsist, through characteristic \textit{critical exponents} whose values differ from the usual predictions of the mean field theories. In other words, mean field theory breaks down because of the collective phenomena observed at the critical point. Only the renormalized field theoretic description of the system can correctly describe the evolution of observables in the vicinity of the critical point.
\end{itemize}

{\color{red} Some consequences are discussed below bla bla   }

\subsection{ 0D model }

{\color{red}
\begin{itemize}
    \item Mean number of neutrons : can we simply use the $\beta$ exponent ? $N \propto (-r)^\beta$ with
       \begin{equation}
           \beta = 1 - \frac{\epsilon}{6} (\epsilon= d - d_c)
       \end{equation}
       No experimental consequence can be envisaged. Neither r nor sigma can be used to verify the proportinality between power and r/sigma. A calculation with an external source maybe ? Can the RG procedure be done with an external source ?

     \item Variance : ??

     \item Comment eventually on percolation probability (= extinction probability) in the 0D model...

     \item Any other observable ? Relaxation time might have something to do with the $\alpha$ exponent

     \item  A model with precursors might be a nice refinement : might not be multi-species DP though because the precursors are not affected by the Doppler (?).
\end{itemize}
}
\subsection{3D observables}

{\color{red}

\begin{itemize}
    \item Mainly spatial correlations ? How to "use" the correlation length $\xi$ properly (what's the shape, $\exp( r/\xi)$ ?)? Can we conclude on neutron clustering ?

    \item Change in the typical cosine profile observed in reactors ? Is it possible to have an idea about that ?

    \item Boundary DP  ? \\
    Janssen H K, Schaub B, and Schmittmann B 1988 Z. Phys. B 72 111\\
    Howard M, Frojdh P, and Lauritsen K B 2000 Phys. Rev. E 61 167
\end{itemize}
}

\subsection{The DP conjecture}

{\color{red}
The critical behaviour of the problem might be different with a more detailled modelling of the neutron population : energy, exponential jumps etc, finite size effect and quenched disorder...

The DP conjecture \cite{Hinrichsen2000} state that the DP transition quite generally reveals itself in every systems for which
\begin{itemize}
    \item The model displays a continuous phase transition from a fluctuating active phase (neutrons alive) into a unique absorbing state (no neutrons).
    \item The transition is characterized by a positive one-component order parameter (survival probability).
    \item The dynamic rules involve only short-range processes.
    \item The system has no special attributes such as additional symmetries or quenched randomness
\end{itemize}

So that it might indeed shows up in neutronics ...
}

\section{Explicit form of the propagators in the two species model}
\label{app_b}
Setting $\Omega = 4 \beta \lambda_D + \lambda_D^2 + 2 \lambda_D r + r^2$, one has,
\begin{eqnarray}
    \Sigma_{NN}(t,t') &=& \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \left\{  \sqrt{\Omega} \cosh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right) \right. \nonumber\\
    & & \left.  + (r + \lambda_D) \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right)  \right\},
\end{eqnarray}

\begin{eqnarray}
    \Sigma_{MM}(t,t') &=& \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \left\{  \sqrt{\Omega} \cosh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right) \right. \nonumber\\
    & & \left.  - (r + \lambda_D) \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right)  \right\},
\end{eqnarray}

\begin{equation}
        \Sigma_{NM}(t,t') = 2 \beta \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right),
\end{equation}

\begin{equation}
    \Sigma_{MN}(t,t') = 2 \lambda_D \frac{e^{ \frac{1}{2} (r - \lambda_D)(t - t')  }}{\sqrt{\Omega}} \sinh\left(  \frac{1}{2} \sqrt{\Omega} (t - t') \right) .
\end{equation}

\section{Two species model - The importance of being delayed}

Another key element in the modelling  of the temporal behaviour of the neutron population in a nuclear reactor is the role played by the precursor population of delayed neutrons $M(t)$. Precursors are unstable fission products resulting from the fission of heavy nuclides. Their further decay, with rate $\lambda_D$, into a stable isotope result in the emission of a fraction $\beta$ of the total neutron population.

Including the coupling between the neutron and precursor populations, the point-like model studied up to now becomes (with $r = \rho - \beta)$,

\begin{equation}
\begin{cases}
        \cfrac{d N(t)}{dt} & = r N(t) + \lambda_D M(t) + \eta(t), \\
         \\
    \cfrac{d M(t)}{dt} & = \beta N(t) - \lambda_D M(t) + \zeta(t) ,
\end{cases}
\end{equation}
where external sources have been neglected, for simplicity. For the same reason, and although this does not correspond to what would have been expected on physical ground, the stochastic terms $\eta$ and $\zeta$ will be considered independent in the following.

The derivation of the response functional associated to this coupled system of equations follows exactly the recipe described in section \ref{sec:FieldFormulation}. Four independent fields $\ensuremath{\widetilde{N}} \ensuremath{\widetilde{M}} M$ and $N$ are now introduced into the theory and the response functional can be written $\mathcal{S} = \mathcal{S}_{\text{neut.}} + \mathcal{S}_{\text{prec}}$, with,
\begin{eqnarray}
     \mathcal{S}_{\text{neut.}} & = & \int \left\{\ensuremath{\widetilde{N}}(t)  \left( \frac{d}{dt} - r \right)N(t) \right. \nonumber\\
     & & \left. -  \lambda_D \ensuremath{\widetilde{N}}(t) M(t) - \lambda_1 \ensuremath{\widetilde{N}}^2(t) N(t) \right\} dt
\end{eqnarray}
and
\begin{eqnarray}
     \mathcal{S}_{\text{prec.}} & = & \int \left\{ \,\ensuremath{\widetilde{M}}(t)  \left( \frac{d}{dt} + \lambda_D \right)M(t) \right. \nonumber\\
     & & \left. - \beta \ensuremath{\widetilde{M}}(t) N(t) - \lambda_2 \ensuremath{\widetilde{M}}^2(t) M(t) \right\} dt.
\end{eqnarray}
The coupling between these two functional renders the perturbative calculations of observables much more involved. One has for instance to take into account the decay of precursors into neutrons or the creation of precursors following a fission event. This translate into a more general form for the basic objects of the theory: the propagators that must be taken into account now reads,
\begin{widetext}
\begin{eqnarray}
    \Sigma_{NN}(t,t')  & =&  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
    (a) --  (b) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture} =
     \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a](o1);
    \diagram* {
    (a) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.5 cm of c](o1);
    \diagram* {
    (a) -- (b) ,
    (b)--[scalar] (c),
    (c) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \left( \beta \, \lambda_D  \right)^2
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.7 cm of c,dot](d) {};
      \vertex [right=0.7 cm of d,dot](e) {};
      \vertex [right=0.5 cm of e](o1);
    \diagram* {
    (a) -- (b) ,
    (b)--[scalar] (c),
    (c) -- (d) ,
    (d)--[scalar] (e),
    (e) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + ... \, , \, \nonumber\\ & & \nonumber\\
  & = & e^{r(t-t')} + e^{r(t-t')} \beta \lambda_D \int_{t'}^t dt_1 e^{-(r + \lambda_D) t_1} \int_{t'}^{t_1} dt_2 e^{(r+\lambda_D) t_2} + ...\, , \, \nonumber\\ & & \nonumber\\
  & = &  e^{r(t-t')} \left\{1 + \int_{0}^{t-t'}\sqrt{\frac{\beta\lambda_D(t-t'-\tau)}{\tau}}e^{-(r + \lambda_D)\tau} I_1\left(2\sqrt{\beta\lambda_D \tau (t-t'-\tau)}\right) d\tau \right\} \\\nonumber \\ \nonumber \\
  & = & e^{ r (t-t')} e^{-(r+\lambda_D)(t - t')/2} \left\{ \cosh\left( \frac{ (r+\lambda_D)(t-t')}{2} \sqrt{ 1 + \frac{4 \beta \lambda_D}{(r+\lambda_D)^2}   }\right) \right\}\\
    \Sigma_{MM}(t,t')  & =&  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
    (a) -- [scalar] (b) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture} =
     \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a](o1);
    \diagram* {
    (a) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.5 cm of c](o1);
    \diagram* {
    (a) --[scalar] (b) ,
    (b)-- (c),
    (c) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \left( \beta \, \lambda_D  \right)^2
  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,dot](b) {};
      \vertex [right=0.7 cm of b,dot](c) {};
      \vertex [right=0.7 cm of c,dot](d) {};
      \vertex [right=0.7 cm of d,dot](e) {};
      \vertex [right=0.5 cm of e](o1);
    \diagram* {
    (a) --[scalar] (b) ,
    (b)-- (c),
    (c) --[scalar] (d) ,
    (d)-- (e),
    (e) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + ... \, , \, \nonumber\\ & & \nonumber\\
  & = & e^{-\lambda_d(t-t')} \left\{1 + \int_{0}^{t-t'}\sqrt{\frac{\beta\lambda_d(t-t'-\tau)}{\tau}}e^{(r + \lambda_d)\tau} I_1\left(2\sqrt{\beta\lambda_d \tau (t-t'-\tau)}\right) d\tau \right\} \\ \nonumber \\ \nonumber \\
    \Sigma_{NM}(t,t')  & =&  \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
    (a) -- [scalar] (b) --(o1)
    };
    \end{feynman}
  \end{tikzpicture} = \beta \left\{ \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b](o1);
    \diagram* {
    (a) --[scalar] (b),
    (b) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \, \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b,dot](c) {};
      \vertex [right=of c,dot](d) {};
      \vertex [right=of d](o1);
    \diagram* {
    (a) --[scalar] (b) ,
    (b)-- (c),
    (c)-- [scalar] (d),
    (d) -- (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + ...\right\} \nonumber\\ & & \nonumber\\
  & = & \beta e^{-\lambda_d(t-t')} \int_{0}^{t-t'} e^{(r + \lambda_d)\tau} I_0\left(2\sqrt{\beta\lambda_d \tau (t-t'-\tau)}\right) d\tau
  \end{eqnarray}
\end{widetext}
 \begin{widetext}
 \begin{eqnarray}
   \Sigma_{MN}(t,t') & = & \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[small, inline=(a)]
      \vertex (a);
      \vertex [right=0.5 cm of a,blob](b) {};
      \vertex [right=0.8 cm of b](o1);
    \diagram* {
    (a) --  (b) -- [scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture} = \lambda_D \left\{ \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b](o1);
    \diagram* {
    (a) -- (b),
    (b) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
  + \beta \, \lambda_D \, \,
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right=of a,dot](b) {};
      \vertex [right=of b,dot](c) {};
      \vertex [right=of c,dot](d) {};
      \vertex [right=of d](o1);
    \diagram* {
    (a) -- (b) ,
    (b)-- [scalar](c),
    (c)-- [] (d),
    (d) --[scalar] (o1)
    };
    \end{feynman}
  \end{tikzpicture}
 + ...\right\} \nonumber\\ & & \nonumber\\
  & = &  \lambda_d e^{r(t-t')} \int_{0}^{t-t'} e^{-(r + \lambda_d)\tau} I_0\left(2\sqrt{\beta\lambda_d \tau (t-t'-\tau)}\right) d\tau,
\end{eqnarray}
\end{widetext}
where $I_0$ and $I_1$ are the modified Bessel function of the first kind. These propagators are the building block of the theory and they must be taken into account for the calculation of observables. Imposing in the following the initial conditions $( N(0), M(0)) =(N_0, 0 )$, the mean number of neutrons and precursors can directly be written,
\begin{equation}
    \langle N(t) \rangle = N_0 \times \Sigma_{NN}(t,0)
\end{equation}
\begin{equation}
    \langle M(t) \rangle = N_0 \times \Sigma_{NM}(t,0)
\end{equation}

Imposing the initial condition $( N(0), M(0)) =(N_0, 0 )$, the variance of the neutron population may symbolically be written,
\begin{eqnarray}
    & & \langle N^2(t) \rangle_c   =  \nonumber\\ & & 2 \lambda_1 N_0
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- (f1)  -- (a) --  (f2) -- (ff2),
    (a) -- (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture} + 2 \lambda_2 N_0
      \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5 cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) --  (f1)  -- [scalar] (a) -- [scalar] (f2) -- (ff2),
    (a) -- [scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
\end{eqnarray}
whereas for the precursors, one has,
\begin{eqnarray}
    & & \langle N^2(t) \rangle_c   =  \nonumber\\ & & 2 \lambda_1 N_0
    \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- [scalar] (f1)  -- [] (a) -- []  (f2) -- [scalar] (ff2),
    (a) -- [] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture} + 2 \lambda_2 N_0
      \begin{tikzpicture}[baseline=(a)]
    \begin{feynman}[inline=(a)]
      \vertex (a);
      \vertex [right= 0.5 cm of a,blob](b) {};
      \vertex [above left=0.5 cm of a,blob](f1) {};
      \vertex [below left=0.5 cm of a,blob](f2) {};
      \vertex [above left=0.7 cm of f1](ff1);
      \vertex [below left=0.7 cm of f2](ff2);
      \vertex [right=0.7 cm of b](o1);
    \diagram* {
     (ff1) -- [scalar]  (f1)  -- [scalar] (a) -- [scalar] (f2) --[scalar]  (ff2),
    (a) -- [scalar] (b) --  (o1)
    };
    \end{feynman}
  \end{tikzpicture}
\end{eqnarray}