\section{Introduction}

The basic operations of memory---imprinting, reading, and erasure---can occur rapidly.  Examples are the flipping of a magnetic domain, writing a nested sequence of values in return-point memory \cite{Sethna1993}, or aging and rejuvenation in glasses \cite{Jonason1998, Zou2010}.
The situation is drastically different for pulse-duration memories in traveling charge-density waves \cite{Fleming1986, Brown1986}, where a memory can be gradually encoded into or erased from the steady-state response of the system. If a series of equal-width current pulses is applied to the system, the nonlinear voltage response eventually becomes phase-locked to the \textit{ends} of the pulses. Thus, the system has a memory of the pulse duration.

Moreover, charge-density wave carriers can simultaneously store multiple memories. If a sequence of pulses is repeatedly applied, the system will remember each of the pulse widths. Counterintuitively, all but two will be forgotten in the steady state. However, if sufficient noise is present, all of the memories can be retained \cite{Coppersmith1997, Povinelli1999}. The defining features of multiple transient memories are thus: (i) initially, the system can learn multiple inputs, (ii) under continual application of these inputs, the system will forget almost all of them, and (iii) if sufficient noise is present, the system will remember the inputs indefinitely.

Recently, Keim \& Nagel \cite{Keim2011} showed how multiple transient memories could be generalized in a generic driven disordered system with noise. Turning to classical driven systems with disorder, they noted other examples of memories stored in a steady state. For example, granular materials that are driven by tapping \cite{Knight1995} or shearing \cite{Toiya2004} assume a steady-state density that is determined by the amplitude of the driving. With prior knowledge of this compaction behavior, simply measuring the pack height after many driving cycles identifies the driving parameter uniquely---i.e.,\ the system has a single memory in its steady state.

To demonstrate that particulate systems could possess multiple transient memories, Keim \& Nagel \cite{Keim2011} employed a simplified model that was developed by Cort\'e \textit{et al.}~\cite{Corte2008} as a kinematic description of experiments of non-Brownian suspensions under cyclic, low Reynolds-number shear \cite{Pine2005, Corte2008}. In such systems, particles were found to follow irreversible trajectories before self-organizing into a reversible steady state, where the particles retrace their trajectories exactly during every subsequent shear cycle.

In this system, a memory may be encoded by shearing a system cyclically between strains $\gamma = 0$ and $\gamma_1$. As the system self-organizes, more and more particle trajectories will become reversible for strains smaller than or equal to the strain amplitude $\gamma_1$. A memory is encoded as the sharp increase in particle irreversibility for shears infinitesimally larger than $\gamma_1$~\cite{Keim2011}. A configuration that is a reversible steady state for amplitude $\gamma_1$ is also reversible for any smaller shear amplitude, $\gamma<\gamma_1$. Hence, if multiple memories are stored in the system by shearing to multiple different amplitudes $\{\gamma_1, \gamma_2, ... \gamma_n\}$, $\gamma_n$ being the largest, then once the system reaches the final reversible steady state (i.e.,\ complete reversibility for shears up to $\gamma_n$), the smaller memories will have been forgotten. There then remains only the onset of irreversibility at $\gamma_n$.  We note here that the system has effectively two memories in the final steady-state configuration (in analogy with the charge-density-wave system).  One of these is, of course, at $\gamma_n$; the other is the maximum strain in the opposite direction, which in this work is fixed at $\gamma = 0$.  These values mark the two endpoints of the driving during the training cycles.

Crucially, if noise is present in the system, the suspension can never reach its final fixed point where all the particle trajectories would be fully reversible up to a strain amplitude $\gamma_n$. Thus the system explores a space of metastable states, and retains signatures of the inputs indefinitely. Because it is effectively maintained in a permanent transient, if the pattern of inputs is changed, the system will gradually acquire new memories and shed the old ones.  The ability to retain a memory of the small inputs is never lost.  This is a concrete example of the emergence of plasticity in memory formation.

While multiple transient memories are likely a generic phenomenon, a rigorous understanding of the conditions under which they are feasible has yet to have been expressed. Here, we address some of the issues underlying this question, by further exploring the sheared suspension system and focussing on three main themes: robustness, the structure of memories, and overdriving. We begin by describing the basic simulation algorithm in section II. In section III, we show that multiple transient memories are robust and are present under several different training algorithms. We then show that forgetting is gradual---a crucial property that distinguishes transient memories from other classes of memory. We demonstrate this over a range of parameters and algorithm variants by reading out a small memory after a single large shear has been applied. In section IV, we examine the spatial structure of particle configurations that possess memories in the Cort\'e algorithm, and we show how this and other effects lead to a broadening of the memory signature.

Finally, in section V, we examine sheared suspensions driven above the critical amplitude for irreversibility, $\gamma_c$, the largest strain amplitude at which the system can self-organize to a reversible steady state~\cite{Corte2008}. We show that even in the case of overdriving, memories can be retained in the system just as for smaller amplitude inputs.  Thus memories are not necessarily relegated to small amplitudes.  Morevover, driving above $\gamma_c$ can be harnessed as a source of noise that allows memories to be retained indefinitely.

The picture that emerges from this work is that multiple transient memories are a robust feature manifest in a range of simple models of suspensions under cyclic, low Reynolds-number shear, that the details of memory retention and erasure can be understood from the spatial structure of the particles, and that overdriving can provide another avenue for controlling memory formation.

\section{Description of Simulations}
\label{sec:algorithm}

We use a variety of models to study the rearrangements of particles in a viscous, non-Brownian suspensions under cyclic low Reynolds-number shear. One of these models was designed by Cort\'e \textit{et al.}~\cite{Corte2008}. A slight variant was used by Keim \& Nagel \cite{Keim2011} to show multiple transient memories in sheared suspensions. Although the models are very simplified, they are justified on the basis that they reproduce many behaviors seen in previous experiments \cite{Corte2008}.

The algorithm of Cort\'e \textit{et al.}~\cite{Corte2008} consists of three steps, illustrated in Fig.\ \ref{algorithm}: (1) $N$ discs of diameter $d$ are randomly placed in a square box of area $A_{box}$ with area fraction $\phi=N\pi (d/2)^2/A_\text{box}$. We choose $\phi=0.2$ and use periodic boundary conditions at the box edges. (2) The box is sheared affinely in the $x$ direction to an amplitude $\gamma_1$ so that a particle is translated a distance $\Delta x = \gamma_1 y$, where $y$ is the vertical coordinate of the particle center.  The particles are then returned to their original positions. (3) Any discs that overlap each other at any point during the shear (i.e.,\ if their centers come within a distance $d$) are given small, random displacements (or ``kicks''), once they have been returned to their unsheared positions. In each variant of the algorithm that we study, a different choice for the displacements is prescribed.

\begin{figure}
\includegraphics[width=2.5in]{algorithm.pdf}
\caption{
Simulation algorithm of Cort\'e \textit{et al.}~\cite{Corte2008}. $N$ discs of diameter $d$ are randomly placed in a square box with area fraction $\phi$ (left). During each cycle, the box is sheared affinely by translating the particle centers a distance $\Delta x= \gamma_1 y$ (right). They are then returned to their previous positions (left). Any particle pairs that overlap at any point during the cycle are given small displacements (left, dashed circles), which determine the new positions for the next shear cycle.  Figure adapted from \cite{Corte2008}.
\label{algorithm}
}
\end{figure}

This very simple model reproduces several important features of the sheared-suspension experiment. First, for small shear amplitude (where the definition of ``small'' depends on the area fraction), the system will eventually fall into a reversible steady state (i.e.,\ no particle pairs come within a distance $d$ of each other during the shear to $\gamma_1$). Second, for large enough shears, the particles never find a reversible configuration (on the timescale of the simulation), and instead each particle undergoes a random walk indefinitely. Third, there is a critical $\gamma_c$ (which depends on the area fraction, $\phi$) that separates these two regimes.  These features were also seen in the experiment \cite{Corte2008}.

Recently, Keim \& Nagel \cite{Keim2011} showed how such a system that is driven below its critical amplitude ($\gamma<\gamma_c$) can store multiple memories in its transient state, or in a steady state when random noise is added. In the current study, we will consider the different choices for the kicks applied to the overlapping particles, as well as different forms of driving.  We also consider the effect of overdriving, that is the effect of shears of amplitude $\gamma>\gamma_c$. In the present work $\phi=0.2$, and empirically $\gamma_c = 4.0$. Without loss of generality, we set the particle diameter $d=1$.

\section{Robustness of Multiple Transient Memories}

The capacity of a system with many metastable states and many degrees of freedom to ``learn'' multiple memories under cyclic driving appears to require dynamics with very few essential features~\cite{Keim2011}. Here we test and study the generality of this learning behavior under different algorithms. First, we evaluate several variants of the original simulation algorithm. We then consider one aspect of the kinematics that is necessary for multiple transient memories: the persistence of a memory over many cycles.

\subsection{Dependence on Kinematics}
\label{theme1}

\begin{figure}
\includegraphics[width=2.9in]{alt-kinematics.pdf}
\caption{
    (color online) Multiple memories observed using variants of the kinematics. Memory signal, $f_\mathrm{mov}''$, is plotted versus trial strain, $\gamma$. Red triangles above and below curves indicate values used in the training pattern. \textbf{(a)} Variant~\ref{variant:corte} in text, with an isotropic random kick for each collision, after $\sim 10^4$ cycles (average of 12 runs)~\cite{Corte2008}. \textbf{(b)} Variant~\ref{variant:tko}: each particle kicked at most once per cycle, regardless of number of collisions~\cite{Keim2011}, after $\sim 10^4$ cycles (12 runs); \textbf{(c)} Variant~\ref{variant:momcons}: equal and opposite kicks, after $\sim 10^4$ cycles (12 runs, actual height of clipped peak, $13$, trough, $-8$); \textbf{(d)} Variant~\ref{variant:repulsion}: equal repulsive kicks after $\sim 10^3$ cycles (12 runs, actual height of clipped peak, $62$, trough, $-62$); \textbf{(e)} Variant~\ref{variant:dilation}: Isotropic repulsion (simulating particle swelling) after 200 cycles (5 runs, height of trough: $-234$). A single cycle is plotted for each variant, taken from a period in that system's evolution when both memories are observable after each cycle, independent of which $\gamma_i$ was just applied.
\label{fig:variants}
}
\end{figure}

To test whether memories are present after some number of cycles, we evaluate $f_\text{mov}''(\gamma) \equiv d^2 f_\text{mov} / d\gamma^2$, where $f_\text{mov}(\gamma)$ is the fraction of particles that collide in the course of a trial deformation by $\gamma$. The quantity $f_\text{mov}(\gamma)$ is cumulative, including all particles colliding at shears up to $\gamma$, whereas $f_\text{mov}'(\gamma)$ measures the set of particles colliding at $\gamma$ specifically. The ``training'' from repeated shearing therefore depletes $f_\text{mov}'$ just below the training value and enhances it above. This change is observed as a peak in $f_\text{mov}''$ for each $\gamma_i$ when multiple transient memories are present, as in Fig.~\ref{fig:variants}(a). To find $f_\text{mov}''$, we sample $f_\text{mov}$ at intervals of 0.04 in $\gamma$ (except in Variant~\ref{variant:dilation} below). We then take the discrete first and second derivatives by subtracting the values at neighboring points.

Our simulations suggest that the multiple-memory behavior is insensitive to qualitative changes in the kinematics. We have assessed the robustness of the multiple-memory behavior under 5 variants of the kinematics, all with $\phi = 0.2$, $N=10^4$. We average over multiple runs, each with a different random initial condition. Representative results from each variant are plotted in Fig.~\ref{fig:variants}.

\begin{enumerate}
    \item Original kinematics of Cort\'e \textit{et al.}~\cite{Corte2008} (see Sec.~\ref{sec:algorithm}): \label{variant:corte} If, during a shear deformation cycle, two particles overlap, then each particle is given a kick with random direction, and magnitude drawn from a uniform distribution on [0, $\epsilon$] where here $\epsilon = 0.005$.
Particles are sheared cyclically with $\gamma_1 = 2$ and $\gamma_2 = 3$, and the smaller shear is applied 5 times for each application of the large shear.  Thus the pattern of shears is $3, 2, 2, 2, 2, 2$, and in Fig.~\ref{fig:variants}(a) this 6-cycle pattern is repeated for a total of $\sim 10^4$ cycles.

    \item ``Tag-kick-once''~\cite{Keim2011}: \label{variant:tko} We use the same algorithm as in  Variant~\ref{variant:corte} except that each colliding particle is kicked only once, independent of the number of other particles that it overlapped. The kick is drawn from the same distribution as in Variant~\ref{variant:corte}.

    \item Momentum conservation: \label{variant:momcons} Kicks to interacting particles are still drawn from the same isotropic random distribution as Variant~\ref{variant:corte}, and all parameters are the same. However, each particle in a colliding pair receives an equal and opposite kick. As in Variant~\ref{variant:corte}, kicks from multiple collisions are applied additively. Thus if any subset of the system does not interact with particles outside that subset, its center of mass will not move.

    \item Repulsion: \label{variant:repulsion} Particles in an interacting pair receive kicks away from each other, as defined by their positions at $\gamma = 0$. Kick magnitudes are equal, and are taken from the same distribution as in Variant~\ref{variant:corte}. Notably, we find that the system evolves at an approximately normal pace (including forgetting of the memory of the smaller shear value), yet it usually does not reach a reversible state within $2 \times 10^6$ cycles (for these parameters).

    \item Dilation and repulsion: \label{variant:dilation} Instead of a shearing motion, kinematics simulate uniform cyclic swelling of the particles to a diameter $\gamma d$.  As in the variants described above, we used a training pattern with two values of $\gamma$, with the smaller value repeated 5 times for each application of the larger value: $\gamma_i =$ 1.5, 1.1, 1.1, 1.1, 1.1, 1.1.  We use $\epsilon = 0.01$. Kicks to overlapping particles are repulsive and have equal randomly-selected magnitude, as in Variant~\ref{variant:repulsion}. For this algorithm we find that $\gamma_c \approx 1.8$, and we sample $f_\text{mov}$ at intervals of 0.02.

\end{enumerate}

As seen in Fig.~\ref{fig:variants}, each of these variants supports multiple transient memories. Our results thus show that the phenomenon is not sensitive to how particles interact, or even to the geometry of deformation.

\subsection{Gradual Forgetting}
\label{sec:gradual}

One important feature of transient memories is that the process of forgetting is gradual. This is observed in simulations of non-Brownian sheared suspensions \cite{Keim2011} and simulations of traveling charge-density waves \cite{Coppersmith1997, Povinelli1999}. For sheared suspensions, gradual forgetting is the property that a memory of $\gamma_i$ will persist even after one or more applications of a larger shear, $\gamma_j>\gamma_i$. This feature helps distinguish transient memories from other classes of memories. For example, in return-point memory in ferromagnets, smaller memories are wiped out the instant a larger field is applied \cite{Barker1983, Sethna1993}. The same holds for aging and rejuvenation in glasses; a memory of a temperature $T_i$ is erased as soon as the glass is heated above $T_i$ \cite{Jonason1998, Zou2010}.

In ref.~\cite{Keim2011}, the kick size was small ($0.005 \leq \epsilon \leq 0.1$). Here we show that forgetting is still gradual (albeit more rapid) when larger kicks are used ($0.01\leq \epsilon \leq 1$). Using the Cort\'e algorithm (Variant~\ref{variant:corte}), we train an initially randomized system ($N=10^5$, $\phi=0.2$) with a single strain amplitude, $\gamma_1=2$, to a reversible steady state. Then, we apply a single shear of strain amplitude $\gamma_2=3$, and attempt to read out the memory at $\gamma_1$. We average multiple systems together to improve the signal-to-noise ratio (between 3 and 400 systems), and we sample $f_\text{mov}$ at intervals of $\Delta \gamma = 0.1$.

We show the results in Fig.\ \ref{gradualforgetting}. We plot $f_\text{mov}$ and its curvature, $f_\text{mov}''$, versus strain for each value of $\epsilon$ tested. A peak is clearly visible for small kicks. As $\epsilon$ is increased, the peak becomes broader and shorter. Nonetheless, there is still an identifiable peak in $f_\text{mov}''$ for kicks as large as $\epsilon=0.7$, as shown in the inset to Fig.\ \ref{gradualforgetting}(b). (Examining the $f_\text{mov}$ data at $\epsilon=0.7$ by eye, one might not identify the memory at $\gamma_1=2$. It is remarkable that the memory can be easily identified by taking two derivatives.)
The intuitively appealing picture that emerges is that the smaller the kick size, $\epsilon$, the slower the erasing of the memory at the lower shear value.

\begin{figure}[bt]
\centering
\begin{center}
\includegraphics[width=2.5in]{gradualforgetting.pdf}
\end{center}
\caption{(color online) Gradual forgetting of transient memories. \textbf{(a)} $f_\text{mov}$ versus $\gamma$ for systems that were driven once at a shear amplitude $\gamma_2=3$ after having been completely trained at $\gamma_1=2$. The kick size, $\epsilon$, is given in the legend, and the arrows show the direction of increasing $\epsilon$.
\textbf{(b)} The curvature of the data, $f_\text{mov}''$, shows that the memory at $\gamma_1=2$ is still present, even for large $\epsilon$. As $\epsilon$ is increased, the peak becomes broader and shorter. \textit{Inset}: The data for $\epsilon=0.7$ show that the memory can be read out even at this large a kick size. All data were taken from averaging over multiple initial configurations.}
\label{gradualforgetting}
\end{figure}

One way in which the memory might be surviving is due to the distribution of random kicks used. In the Cort\'e algorithm, each colliding particle is given a kick with a random size chosen uniformly between $0$ and $\epsilon$. Thus, even if $\epsilon$ is large, some particles are given small kicks, and these particles alone might store the smaller memory. To investigate this, we ran simulations where a kicks of fixed size $\epsilon$ are applied (in random directions) to colliding particles. In this case, we found that the memory at $\gamma_1=2$ was still present up to $\epsilon=0.3$.

We tested three other variants of the Cort\'e algorithm, to further test the generality of gradual forgetting. In each case the system is first trained to a reversible state at $\gamma_1=2$ and then a single shear of $\gamma_2=3$ is applied. The algorithms below were used in both the initial training to the reversible state, and in the final, disrupting shear.

\begin{enumerate}
    \setcounter{enumi}{5}
    \item Annealing at zero shear: After every cycle, overlapping particles are given random kicks until the system reaches a state with no overlaps (in the unsheared configuration). Gradual forgetting was observed at $\epsilon=0.2$ (no other values were tested).

    \item Pure shear: Particle centers are sheared along trajectories given by: $\Delta x=x(\cosh(\frac{1}{2}\gamma)-1)+y\sinh(\frac{1}{2}\gamma)$, $\Delta y=y(\cosh(\frac{1}{2}\gamma)-1)+x\sinh(\frac{1}{2}\gamma)$. This is in contrast to simple shear, where $\Delta x=\gamma y$ and $\Delta y=0$, as in the Cort\'e algorithm. (The relationship between pure shear and simple shear is more easily understood in terms of the instantaneous flow field, which is defined in terms of the shear rate, $\dot{\gamma} \equiv d\gamma(t)/dt$. Pure shear has velocity components $v_x=\frac{1}{2}\dot{\gamma} y$, $v_y=\frac{1}{2}\dot{\gamma} x$, whereas simple shear has $v_x=\dot{\gamma} y$, $v_y=0$.) Here gradual forgetting was observed at $\epsilon=0.3$ (larger $\epsilon$ were not tested).

    \item Dilation: The particles are swelled to radius $\sqrt{1+\gamma/\pi}$ (this radius is chosen so that the area of the interaction region matches that for shearing). Gradual forgetting was observed for $\epsilon=0.1$.

\end{enumerate}

Thus, the property of gradual forgetting is not relegated to small kicks and does not depend on the specific properties of the Cort\'e algorithm. It is a general feature of a range of kinematics of non-Brownian particles.

\section{Structure of a memory}
\label{theme2}

\subsection{Pair correlation function}

The two-dimensional (2D) pair correlation function, $g(x,y)$, affords us one method of characterizing the structure of a system with one or more memories. The pair correlation function is proportional to the probablity that two particles will be separated by the displacement $(x,y)$ and is normalized so that $g(x,y)=1$ for a uniform distribution of particles.  A complete memory of an amplitude $\gamma_1$ imposes a constraint on $g(x,y)$. Because there can be no overlapping particles in the unsheared configuration, there can be no particle centers within a radial distance 1 from the origin. Furthermore, a complete memory entails that no particle centers fall in this region when the system is sheared continuously up to $\gamma_1$. This constraint creates a propeller-shaped region where $g(x,y)=0$ for a complete memory, as shown in Fig.\ \ref{structure}(a). The data are averaged over 200 systems with $N=10^4$, $\phi=0.20$, and $\epsilon=0.1$. This excluded region is the set of points satisfying:
\begin{equation}
-\sqrt{1-y^2} \mp \gamma y \leq \pm x \leq \sqrt{1-y^2},
\label{excludedregion}
\end{equation}

\noindent  for $0 \leq \pm y \leq 1$.

\begin{figure}[bt]
\centering
\begin{center}
\includegraphics[width=3.4in]{structure.pdf}
\end{center}
\caption{(color online) Pair correlation function, $g(x,y)$, of trained systems, where $x$ and $y$ are the flow and gradient axes, respectively. \textbf{(a)} Complete memory of $\gamma_1=3$. All particles lie outside the excluded region created by the simple shear deformation. $g(x,y)$ has peaks at $(\pm 1,0)$ where there are cusps in the shape of the excluded region . \textbf{(b)} Two memories at $\gamma_1=2$ and $\gamma_2=3$, trained with 102 cycles following the pattern: $3,2,2,2,2,2$. Particles accumulate on the boundaries of the two regions given by eqn.\ \ref{excludedregion} for $\gamma_1=2$ and for $\gamma_2=3$. This accumulation marks the formation of two memories.  Data in (a, b) were averaged over 200 systems with $N=10^4$, $\phi=0.20$, and $\epsilon=0.1$. \textbf{(c)} Isolated pairs of particles trained with $\gamma_1=3$ (averaged over $1.22\times 10^7$ pairs, with $\epsilon=0.1$). The pair structure for the full simulation (a) and for isolated pairs (c) are qualitatively similar. Differences are caused by three-body effects.}
\label{structure}
\end{figure}

Two excluded regions for different values of $\gamma_i$ will nest: if $\gamma_1<\gamma_2$, then the excluded region for $\gamma_1$ is a subset of the excluded region for $\gamma_2$. This is equivalent to an observation made in the introduction and in ref.~\cite{Keim2011}: in the algorithms explored here, complete reversibility at one value of $\gamma$ dictates complete reversibility for any smaller strain amplitude. When the system has two or more memories, $g(x,y)$ is depleted of particles to a varying degree in each of the corresponding excluded regions. We show this in Fig.\ \ref{structure}(b) for a system driven at two amplitudes: $\gamma_1=2$ and $\gamma_2=3$. The data are averaged over 200 systems with $N=10^4$, $\phi=0.20$, and $\epsilon=0.1$.

In both the partial and complete memories, the pair correlation function has peaks at $(\pm 1,0)$. The peaks can be understood by the following argument. If two particles collide during a shear cycle, they will contribute to the pair correlation function at two points \textit{inside} the excluded region. These particles undergo a random walk until they exit the excluded region. Thus, isolated pairs of particles undergo a 2D diffusion process with an absorbing boundary. The curvature of the boundary of the excluded region of $g(x,y)$ influences the local density of particles exiting the region at that point---large positive curvatures mean few particles will escape there, while large negative curvatures lead to high accumulation at the boundary. The boundary of the excluded region of $g(x,y)$ has concave cusps at $(\pm 1,0)$ (i.e.,\ the curvature is negative infinity), which leads to a divergent density of particles at these points in the limit $\epsilon \to 0$.

This picture ignores three-body effects, which could spread out these peaks. To test the relative strength of three-body effects, we mimic the extreme low area-fraction limit of the Cort\'e kinematics. Namely, we simulate an ensemble of isolated pairs of particles. We evolve the pairs by giving both particles a random kick if they collide (which is equivalent to leaving one particle at the origin and giving the other particle \emph{two} random kicks). The systems are each evolved to a reversible state. We show the results in Fig.\ \ref{structure}(c), where we plot $g(x,y)$ accumulated over $1.22\times 10^7$ systems of 2 particles with $\epsilon=0.1$.
The diffusion picture qualitatively captures the peaks in the pair correlation function at $(\pm 1,0)$. Comparing $g(x,y)$ from the diffusion simulation to the full simulation indicates that three-body effects subdue the concentration of particles along much of the boundary of the excluded region, except near $(\pm 1,0)$.

\subsection{Memory Step Function}

\begin{figure*}[bt]
\centering
\begin{center}
\includegraphics[width=7.0in]{memorystep.pdf}
\end{center}
\caption{(color online)  Step function in the memory. \textbf{(a,b,c)} $f_\text{mov}$ versus strain, $\gamma$, for systems with complete memory of $\gamma_1=2$, versus kick size, $\epsilon$, trained with different algorithm variants: \textbf{(a)} radial swelling of isolated pairs of particles to diameter $\gamma_1$, \textbf{(b)} shearing of isolated pairs of particles, and \textbf{(c)} Cort\'e algorithm. In (a), as $\epsilon$ approaches zero, the memory approaches a step function in $f_\text{mov}$, since all pairs are within $2\epsilon$ of colliding at the end of the simulation. In (b), the step function is smoothed out and the height of the jump shortened. In (c), three-body effects shorten and smooth the jump even more. In all cases, a smaller kick size, $\epsilon$, steepens $f_\text{mov}$ at $\gamma_1$. Data in (a,b) are from $5\times 10^4$ isolated pairs of particles. Data in (c) are single systems of $N=10^5$ particles (except for $\epsilon=0.001$, where $N=10^4$). The arrow in (b) shows the direction of increasing $\epsilon$, which is the same in (a) and (c). \textbf{(d)} The size of the jump in $f_\text{mov}$ at $\gamma_1=2$ for radial swelling of isolated particles ($\bullet$), for shearing of isolated particles ($\vartriangle$), and for the Cort\'e algorithm ($\blacktriangle$). \textbf{(e)} Slope of $f_\text{mov}$ on the right-hand side of $\gamma_1$, scaled by $\epsilon$. In both (d) and (e), the data appear to plateau at a finite value for all three algorithms in the limit of small $\epsilon$. The plateau value indicates the degree to which the behavior is approximated by a step function.}
\label{memorystep}
\end{figure*}

When a system stores a partial or complete memory of some amplitude, many pairs of particles are nearly touching at that shear value. That is, the training has pushed the particles just far enough so that they do not collide at the shear amplitude used in the training.  In examining the pair correlation function, $g(x,y)$, we see that this observation can be stated in the following way: there is a high concentration of particles within a distance $2\epsilon$ from the boundary of the excluded region in $g(x,y)$. Therefore, many particle pairs that behave reversibly for a shear to $\gamma_1$, will collide when the system is driven to $\gamma_1+2\epsilon$. In the case of a complete memory, this means that $f_\text{mov}(\gamma_1)=0$, whereas $f_\text{mov}(\gamma_1+2 \epsilon)$ is finite. Consequently, one might expect a step function in $f_\text{mov}$ at $\gamma_1$ for small $\epsilon$. Counter to this expectation, none of the data reported by Keim \& Nagel \cite{Keim2011} show a step function in $f_\text{mov}$ at the memory amplitudes.

To understand the apparent absence of a step function in $f_\text{mov}$, we consider a simpler system, where a step-function in $f_\text{mov}$ is clearly present. In Fig.\ \ref{memorystep}(a), we show $f_\text{mov}$ for simulations of $5\times 10^4$ isolated pairs of particles that initially overlap, trained by cyclic swelling (to diameter $\gamma$). We vary the kick size from $\epsilon=0.3$ down to $0.001$. The $f_\text{mov}$ data approach a step-function at $\gamma_1=2$ in the limit of small $\epsilon$.

The step function is widened for isolated pairs of particles under \emph{shear}: if a particle lands near the points $(\pm 1, 0)$, a very large shear may be necessary for a collision. For example, a particle that is displaced to the point $(1+\epsilon, 0)$ will not collide with its partner, even for arbitrarily large strain amplitudes. If displaced to $(1+\epsilon,\delta)$ (where $\delta>0$ is small), a shear of $\gamma \approx \epsilon/\delta$ is required for a collision---the smaller the displacement off the y-axis, $\delta$, the larger $\gamma$ must be for a collision. Although this effect is specific to a vanishingly small region of configuration space of particle pairs (i.e.,\ an $\epsilon$ neighborhood of the points $(\pm 1, 0)$), this configuration is where $g(x,y)$ is sharply peaked, as was shown in Fig.\ \ref{structure}. We show $f_\text{mov}$ for $5\times 10^4$ isolated pairs of particles that initially overlap, trained with cyclic shear of strain amplitude $\gamma_1=2$ in Fig.\ \ref{memorystep}(b). The size of the jump is less than unity because some pairs separate such that they will never collide, even for arbitrarily large shear amplitudes (by separating either vertically, or into the first or third quadrants of the domain of $g(x,y)$) .

Finally, three-body effects smooth the step function even more. In Fig.\ \ref{memorystep}(c), we show $f_\text{mov}$ for a single system with $N=10^5$ particles and area fraction $\phi=0.2$ driven at $\gamma_1=2$ to a reversible state, under the Cort\'e algorithm. (We used $N=10^4$ particles for $\epsilon=0.001$, due to the prohibitively long computing times of a larger system.) As $\epsilon$ decreases, $f_\text{mov}$ becomes steeper at $\gamma_1$, but the curve does not have a punctuated jump, except perhaps at the smallest kick size tested, $\epsilon=0.001$.

We have given a qualitative account of the disappearance of the memory step function in $f_\text{mov}$ for a complete memory. To be more quantitative, we articulate two measures of the degree to which a particular $f_\text{mov}$ curve is approximated by a step function. The first measure is the size of the jump in $f_\text{mov}$ for a test shear to $\gamma_1$ versus a test shear to $\gamma_1+2\epsilon$. For a complete memory, the size of the jump is simply $f_\text{mov}(\gamma_1+2\epsilon)$, since $f_\text{mov}(\gamma_1)=0$. We show this value versus $\epsilon$ in Fig.\ \ref{memorystep}(d). Whereas swelling of isolated pairs of particles gives a perfect step (i.e.,\ a jump of unity for all values of $\epsilon$), shearing of isolated pairs and including three-body interactions limit the jump to a smaller value for small $\epsilon$.

The second measure is the slope of $f_\text{mov}$ on the right-hand-side of $\gamma_1$, denoted by $f_\text{mov}'(\gamma_1)_+$.
For small enough $\epsilon$, the dominant length scale that dictates this slope will be $\epsilon$, due to a separation in length scales between the kick size, $\epsilon$, and the particle size. We plot the product $\epsilon f_\text{mov}'(\gamma_1)_+$ in Fig.\ \ref{memorystep}(e). The data plateau to a finite value for small $\epsilon$ under all three algorithms, indicating that the slope approaches this prefactor times $1/\epsilon$. The prefactor indicates how well $f_\text{mov}$ is approximated by a step function, which increases as we take out three-body effects and change to a swelling system.

Both quantitative measures show that there is a step-function contribution to  $f_\text{mov}$ in all cases.  However, this step function is broadened by at least two effects. Our results indicate
that three-body effects are the smaller contribution to the broadening and shortening of the step function, when compared with the effect of the cusp in the excluded region for the shearing algorithm. This suggests that an experimental system that employs particle swelling may have a sharper signature of memory formation (assuming all other effects stay constant).

\subsection{Faster training with $\gamma=3,2,2,2,2,2$ is an artifact of drift}

Finally, we note one surprising feature found when training two memories in the sheared suspension simulations. In ref.~\cite{Keim2011}, when a system of $N=10^4$ particles with area fraction $\phi=0.2$ was trained simultaneously at two strain amplitudes, $\gamma_1=2$ and $\gamma_2=3$, with the pattern of shears, $3,2,2,2,2,2$, the system organized to a fully reversible state in $\sim 30\,000$ cycles, compared with $\sim 55\,000$ cycles if the system was trained with just a single amplitude, $\gamma_1=3$. We have investigated this effect further, and we find that the effect is statistically significant, and grows with system size (we tested system sizes from $N=11$ to $N=10^5$). The effect is  more pronounced in Variant~\ref{variant:tko} (``Tag-kick-once"), versus Variant~\ref{variant:corte} (original kinematics of Cort\'e \textit{et al.}).

We find that the effect may be attributed to the ability of particle pairs to drift: since in those two algorithms two overlapping particles each receive a random kick in a random direction, their center of mass may diffuse. This mechanism speeds up training under a more ``gentle'' protocol of $\gamma_1=2$ interspersed with $\gamma_2=3$, which allows higher-density regions to relax while minimizing the effect of particle pairs wandering off and disturbing regions that had already been trained to a reversible state. When we incorporate momentum conservation into these algorithms (Variant~\ref{variant:momcons}, and another variant similar to ``Tag-kick-once'' and including momentum conservation) so that diffusion of interacting particles is suppressed, we find that the effect is almost completely eliminated.

\section{Driving above the threshold for irreversibility}
\label{theme3}

\begin{figure}[bt]
\includegraphics[width=3.4in]{over-memories.pdf}
\caption{
(color online) Multiple-memory behavior at strain values above the critical threshold for irreversibility. Here $\gamma_c = 4.0$; $\gamma_1 = 4.3$ and $\gamma_2 = 4.6$ (red triangles at top). Driving above $\gamma_c$ prevents the system from reaching a reversible steady state with only one memory, and so both training inputs are retained indefinitely. The smaller shear value decays negligibly after $\sim 2 \times 10^6$ cycles. The curve plotted here is from cycle $5.0 \times 10^6$ and is an average over 100 initial conditions.
\label{fig:over_memories}
}
\end{figure}

While we have treated different training values of $\gamma$ as more or less equivalent, in fact as $\gamma$ is increased, it takes longer and longer for the system to reach a reversible steady state. Cort\'e \textit{et al.}~\cite{Corte2008} demonstrated that the characteristic time for self-organization diverges at a value $\gamma_c$; this is a critical transition in the conserved directed percolation class~\cite{Menon2009}. For $\gamma_1 > \gamma_c$, reversible (non-interacting) arrangements of the particles may exist, but the system cannot reach them through its self-organization process~\cite{Corte2008}. In the simulations of memory effects reported above and in ref.~\cite{Keim2011}, all $\{\gamma_i\}$ were kept below $\gamma_c$, so that in the absence of external noise, the system would always evolve to a reversible steady state with only a single memory.

However, if we permit one or more $\gamma_i > \gamma_c$ in the training protocol, the system can never reach a reversible steady state; some particles are kicked in each cycle. This restriction bears a resemblance to that imposed by external noise~\cite{Povinelli1999}. Keim and Nagel~\cite{Keim2011} showed that when uncorrelated noise is added to the particles' positions on each cycle, with a scale $\epsilon_\text{noise} \ll \epsilon$, self-organization to reversibility is averted and multiple memories may persist indefinitely. Similarly, we expect that in the absence of noise, a system driven with some $\gamma_i > \gamma_c$ should sustain multiple memories indefinitely, if it supports the formation of multiple memories at all. Because our simulations must run for finite time, we consider memories to be phenomenologically indefinite when they do not decay over timescales much longer than the time required to enter an apparent steady state; in this regime, changes in memory strength are dominated by fluctuations and no decay is discernible.

\subsection{Forming Memories Above $\gamma_c$}

Figure~\ref{fig:over_memories} shows that multiple memories are indeed possible with $\gamma_i > \gamma_c$. We use the same parameters and kinematics as in Fig.~\ref{fig:variants}(b) (Variant~\ref{variant:tko}). The training strain amplitudes are $\gamma_1 = 4.3, \gamma_2 = 4.6$, with training pattern $\gamma_1, \gamma_2, \gamma_2, \gamma_2, \gamma_2, \gamma_2$, repeat\dots, and $\gamma_c = 4.0$. We observe that after $\sim 2 \times 10^6$ cycles, both memories appear to be stable and present in each cycle, in contrast to the ``forgetting'' behavior reported for $\{\gamma_i\} < \gamma_c$ in ref.~\cite{Keim2011}.  (To compute the curve, $f_\text{mov}$ was sampled at intervals of $\Delta \gamma = 0.0429$, so that sampling positions coincided with both training values.)

\subsection{Overdriving as Noise}

\begin{figure}[bt]
\includegraphics[width=3.4in]{over-noise.pdf}
\caption{
(color online) Multiple memories retained by occasional application of a large strain. Training pattern is $\gamma_i =$ 3, 2, 2, 2, 2, 2 (red triangles at top). Parameters and kinematics are the same as in Fig.~\ref{fig:variants}(b), but with $\epsilon = 0.015$. Curves are offset by $1$ for clarity.
    \textbf{(a)} Without application of any strain greater than $\gamma_c = 4.0$, the $\gamma_1 = 2$ memory is essentially lost after $7 \times 10^5$ cycles shown here (average of 48 runs). It is completely gone (by definition) when the simulation reaches a steady state that is reversible for $\gamma_1 = 3$ ($\sim1.2 \times 10^6$ cycles).
    \textbf{(b, c)} Simulation with $\gamma_\text{large} = 100$ applied every $\sim 100$ cycles. Large-amplitude shearing takes the role of noise, preventing organization to a reversible steady state and so indefinitely maintaining plasticity, the system's susceptibility to multiple memories. Curves show memories after $\sim 2 \times 10^6$ cycles, (b) before and (c) after application of the large-amplitude strain. The peaks in (c) at the two training amplitudes are still well defined showing that this protocol ``refreshes'' the system without destroying memories.
Each curve is an average over 96 runs.
\label{fig:over_noise}
}
\end{figure}

A second, more practical use of $\gamma > \gamma_c$ is to simply take the role of external noise in a training protocol that otherwise has only $\{\gamma_i\} < \gamma_c$. We have realized this principle by intermittently applying a $\gamma_\text{large} \gg \gamma_c$, which for suitably small $\epsilon$ and infrequent application, disrupts the particles' self-organized positions to prevent a reversible steady state, but not so much as to completely destroy already-formed memories. We wish to apply this large shear every $\Delta t_\text{large}$ cycles, to mimic noise applied every cycle with a scale $\sim 0.1 \epsilon$, which is known to sustain memories indefinitely~\cite{Keim2011}. We use the scaling of displacement for a random walk, $\Delta x \sim t^{1/2}$, to arrive at $\Delta t_\text{large} \sim 100$. Since in Variant~\ref{variant:tko} of the kinematics multiple collisions in a cycle do not add to displacements, our large strain amplitude need only satisfy $\gamma_\text{large} \gg \gamma_c$; here we choose $\gamma_\text{large} = 100$. For kinematics in which multiple collisions are additive, such as Variant~\ref{variant:corte}, $\gamma_\text{large}$ would have to be selected more carefully to avoid erasure of existing memories, using the considerations outlined in sec.~\ref{sec:gradual}.

The results of overdriving used as noise are shown in Fig.~\ref{fig:over_noise}, where multiple memories are sustained indefinitely, long after an equivalent simulation without large-amplitude shearing reaches a steady state with just one memory. This technique has the advantage that the sustaining ``noise'' is added in the same way that the driving signal is applied, rather than relying on an additional noise mechanism that may be difficult to introduce or control.

\section{Conclusion}

In summary, we have explored multiple transient memories in simulations of sheared suspensions, focussing on three main themes: robustness, structure, and overdriving. We have shown that multiple transient memories are a robust feature that are manifest in a variety of very simple models of suspensions under cyclic, low Reynolds-number shear.  The details of memory formation can be understood from the spatial correlations of the particles.

We find that overdriving can provide another means for controlling memory formation and retention.  Remarkably, memories can be stored, not only for training at $\gamma<\gamma_c$ but also at amplitudes $\gamma>\gamma_c$.  Moreover, we have found that a large shear that exceeds the critical value, $\gamma_c$, can be used as an alternative to noise to allow memories to persist indefinitely.

As demonstrated by Keim \& Nagel \cite{Keim2011} and recapitulated over a wider range of algorithm variants in this work, the requirements for transient memories are relatively minimal. The salient components appear to be a system that relaxes to a reversible state under cyclic driving during a transient period of irreversibility, and an ordering of reversible states, meaning that a system reversible under a strain amplitude $\gamma_1$ is reversible under smaller shears, $\gamma<\gamma_1$. Thus, we expect other sheared disordered systems with these properties to exhibit transient memories as well, such as granular materials \cite{Toiya2004, Mueggenburg2005}, colloids \cite{Ackerson1988, Haw1998}, foams \cite{Lundberg2008}, and filament networks \cite{schmoller10}, and under other driving protocols, such as tapping \cite{Knight1995}. These results, combined with the previous studies of charge-density waves, suggest that multiple transient memories may be a generic feature in a wide class of driven disordered systems.

\begin{acknowledgments}
We thank Thomas A.\ Caswell for assistance in optimizing the simulations. This work was supported by NSF Grant DMR-1105145 and by NSF MRSEC DMR-0820054 (J.D.P.\ and S.R.N.), and by NSF MRSEC DMR-1120901 (N.C.K.).  Use of computation facilities funded by the US Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, Award No.\ DE-FG02-03ER46088 is gratefully acknowledged.
\end{acknowledgments}


\section{Introduction}

The quantum Einstein gravity (QEG) can be a good candidate to describe the gravitational
interactions in the framework of quantum field theory
\cite{Reuter:1996cp,Reuter:2007rv,Reuter:2012id}. The action of the model is integrated
out for all the possible paths of the field variables, which role is now played by the
metrics. The Einstein-Hilbert action of QEG contains only the
cosmological and the Newton couplings \cite{Reuter:2001ag}. The functional renormalization group (RG)
method allows us to perform the path integration systematically for the degrees of freedom, and gives
us the scaling behavior of the couplings starting from the high energy, ultraviolet (UV)
region down to the low energy infrared (IR) regime \cite{Wetterich:1992yh,
Berges:2000ew,Polonyi:2001se}. The method provides an integro-differential
equation for the effective action, the Wetterich equation.
Besides QEG the RG technique is widely used recently in many areas gravitational issues starting
from black holes \cite{Falls:2010he,Koch:2013owa} to cosmological problems \cite{Hindmarsh:2011hx,Kaya:2013bga}.
The QEG is asymptotically safe \cite{Percacci:2007sz,Reuter:2012xf,Nagy:2012ef,Gies:2013pma}, which implies
that there is a UV non-Gaussian fixed point (NGFP) in the phase space with a finite number of
relevant couplings.

The Wetterich equation contains a regulator function to remove the divergences of the momentum
integrals. In the IR limit the effective action should not depend on
the regulator, because it is an artificial term that is put by hand into the action.
However during the solution of the Wetterich equation we use several approximations
and truncations which may introduce some regulator-dependence. In the RG method mostly the optimized
Litim's regulator is used \cite{Litim:2000ci,Litim:2001up}, especially due to its analytic form. The sensitivity of the
physical quantities on the regulator parameters should be minimal \cite{Canet:2002gs}, since the regulator itself
is an artificial element of the RG method and the results should be independent on it ideally.

There are several types of regulators which are widely used, e.g. the power-law, the exponential
and the Litim ones, which can be get by the limiting cases of the so-called compactly supported
smooth (css) regulator \cite{Nandori:2012tc,Nandori:2013nda}. This regulator provides us such a wide class of regulators
which contain all the important types of regulator functions and enables us to investigate
the sensitivity of the physical quantities on a broader range of set of regulators and their parameters.

We calculate the critical exponent $\nu$ of the correlation length and the anomalous dimension $\eta$
around the fixed points of QEG. We obtain that $\eta=2$ independently of the regulators at the UV NGFP, but the value of
the other exponent $\nu$ can be arbitrary there. The strong truncation of the action might cause such a nonuniversal behavior.
Using the css regulator, we also calculate the value of the exponent $\nu$ for the 3-dimensional $O(1)$ model.
We investigate how the exponents depend on the css regulator parameters and search for those
parameter regions where the exponents are practically parameter-independent or have at least a minimal
sensitivity of them. We also determine how the exponents scale around the IR fixed point.

The paper is organized as follows. In \sect{sec:model} the investigated models, the RG method,
and the regulators are introduced. In \sect{sec:uv} the UV critical behavior is discussed
for the QEG. The crossover criticality is discussed in \sect{sec:co} for the 3-dimensional $O(1)$
model, while we turn back to QEG  to treat its IR criticality in \sect{sec:ir}.
Finally, in \sect{sec:sum} the conclusions are drawn up.

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

The Wetterich equation for the effective average action $\Gamma_k$ is \cite{Wetterich:1992yh}
\beq\label{potev}
\dot \Gamma_k = \hf\mbox{Tr}\frac{\partial_t R_k}{\Gamma''_k+R_k}
\eeq
where the dot denotes the derivative w.r.s. to the `RG time' $t=\ln k$,
furthermore the prime is the differentiation with respect to the field
variable and the trace Tr denotes the integration over all momenta and the summation
over the internal indices. The function $R_k$ plays the role of the IR regulator.
The Einstein-Hilbert effective action is
\beq\label{EHtrunc}
\Gamma_k = \frac1{16\pi G_k}\int d^d x\sqrt{g}(-R+2\Lambda_k),
\eeq
with the dimensionful Newton constant $G_k$ and the cosmological constant $\Lambda_k$.
The RG equations are formulated by the dimensionless couplings, i.e.
$\lambda = \Lambda_k k^{-2}$ and $g = G_k k^{d-2}$. Since the determinant of the metric occurs only in Eq. \eqn{EHtrunc}
and nowhere does it any more in our paper, the usage of $g$ for the Newton coupling below shall not be confusing.
The action contains the first two terms in the Taylor expansion of the curvature $R$.
If one inserts \eqn{EHtrunc} into the Wetterich equation then one obtains the
evolution equations for the couplings, which are derived and given in \cite{Reuter:2001ag} and have the form
\bea
\dot\lambda &=& 2(2-\eta)\lambda+\hf(4\pi)^{1-d/2}g\nn
&&\times [2d(d+1)\Phi^1_{d/2}(-2\lambda)-8d\Phi^1_{d/2}(0)]\nn
&&-d(d+1)\eta\t \Phi^1_{d/2}(-2\lambda)],\nn
\dot g &=& (d-2+\eta)g,
\eea
with the anomalous dimension
\beq
\eta = \frac{gB_1(\lambda)}{1-gB_2(\lambda)}.
\eeq
The functions $B_1(\lambda)$ and $B_2(\lambda)$ are
\bea
B_1(\lambda) &=& \frac13(4\pi)^{1-d/2}[d(d+1)\Phi^1_{d/2-1}(-2\lambda)\nn
&&-6d(d-1)\Phi^2_{d/2}(-2\lambda)-4d\Phi^1_{d/2-1}(0)\nn
&&24\Phi^2_{d/2}(0)],\nn
B_2(\lambda) &=& -\frac16(4\pi)^{1-d/2}[d(d+1)\t \Phi^1_{d/2-1}(-2\lambda)\nn
&&-6d(d-1)\t\Phi^2_{d/2}(-2\lambda)],
\eea
with the threshold functions
\bea\label{thrs}
\Phi_n^p(\omega) &=& \frac1{\Gamma(n)}\int_0^\infty dy y^n\frac{r-yr'}{(y(1+r)+\omega)^p},\nn
\t\Phi_n^p(\omega) &=& \frac1{\Gamma(n)}\int_0^\infty dy y^n\frac{r}{(y(1+r)+\omega)^p},
\eea
where $y=p^2/k^2$ and $r=r(y)$ is the dimensionless regulator $r=R/p^2$.
The dimensionless css regulator has the form
\beq
\label{rcss}
r_{css} = \frac{s_1}{\exp[s_1 y^b/(1 -s_2 y^b)] -1} \theta(1-s_2 y^b),
\eeq
where $b\ge 1$ and $s_1$, $s_2$ are positive parameters.
Although this regulator cannot provide analytic form of the evolution equations,
it gives a broader range of regulators.
The limiting cases of the css regulator provide us the following commonly used regulator functions
\bea\label{limregs}
\lim_{s_1\to 0} r_{css} &=& \left(\frac{1}{y^b} -s_2\right) \theta(1-s_2 y^b),\nn
\lim_{s_1\to 0,s_2\to 0} r_{css} &=& \frac1{y^b},\nn
\lim_{s_2\to 0} r_{css} &=& \frac{s_1}{\exp[s_1 y^b]-1}.
\eea
where the first limit gives the Litim's optimized regulator when $s_2=1$, the second gives the
power-law regulator, and the third gives the exponential one, if $s_1=1$.
It provides us a possibility of simultaneous optimization
of some physical quantities among the Litim's, the power-law and the exponential
regulators which can be continuously deformed from one to the other by only
two parameters $s_1$ and $s_2$. We note that the $b=1$ case satisfies the normalization
conditions \cite{Reuter:2001ag}
\beq\label{reglim}
\lim_{y\to 0} yr = 1 ~\mbox{and}~\lim_{y\to\infty} yr = 0.
\eeq

\section{Ultraviolet criticality}\label{sec:uv}

QEG has two phases, as is shown in \fig{fig:qegphase}.
\begin{center}
\begin{figure}[ht]
\epsfig{file=fig1.ps,width=6cm,angle=-90}
\caption{\label{fig:qegphase}
The phase structure of QEG is shown for css regulator with the parameters $s_1=s_2=b=1$.
The grey points denote the critical points. The thick line represents the separatrix.
}
\end{figure}
\end{center}
The phase space contains 3 fixed points. The UV NGFP is said to be UV attractive,
however it repels the trajectories from the viewpoint of the RG flow. The trajectories start
there and flow towards the hyperbolic Gaussian fixed point (GFP). The trajectories with positive $\lambda$ in the
IR limit constitute the (symmetry-) broken phase of the model which contains an attractive IR fixed point
\cite{Donkin:2012ud,Christiansen:2012rx,Nagy:2012rn,Demmel:2012ub}. The
other trajectories belong to the symmetric phase. One can linearize the RG equations around
the fixed points. In the case of the UV NGFP the eigenvalues of the corresponding stability matrix
can be written as
\beq\label{uveig}
\theta_{1,2}^{UV} = \theta'\pm i\theta'',
\eeq
where the real part of the exponent can be related to the critical exponent $\nu$
of the correlation length, i.e., $\nu=1/\theta'$. We calculated
the position of the fixed point and the corresponding exponent in QEG with Einstein-Hilbert
action in $d=4$. In \fig{fig:sur} we plotted how the reciprocal of the exponent
$1/\nu$ depends on the parameters of the css regulator. The figure demonstrates the appearing
limiting regulators in \eqn{limregs}.
\begin{center}
\begin{figure}[ht]
\epsfig{file=fig2.ps,width=6cm,angle=-90}
\caption{\label{fig:sur}
The reciprocal of the exponent $1/\nu$ is calculated for various parameters of
the css regulator, with $b=3$. The value of $1/\nu$ at the corner $s_1=s_2=0$ corresponds to
{\bf P}ower-law regulator result, the corner $s_1=0,~s_2=1$ gives the {\bf L}itim's regulator
result and the corner $s_1=1,~s_2=0$ results in the {\bf E}xponential regulator value.
}
\end{figure}
\end{center}
We look for the extremum of $1/\nu$ as the function of $s_1$, $s_2$ and $b$. The increase of
the parameters $s_1$ and $s_2$ gives monotonically increasing exponents, their large
values give logarithmically growing exponents. Therefore, we have to look for a minimum of $1/\nu$.
The limit of the power law regulator $s_1=s_2=0$ gives the minimal value for $1/\nu$, when $b=3$.
The exponent also grows with increasing $b$, therefore one might expect a local extremum at $b=1$.
We calculated $1/\nu$ at the Litim's limit ($s_1\to 0$)
for different values of $s_2$ is shown in \fig{fig:qegnulin} for $b=1$.
\begin{center}
\begin{figure}[ht]
\epsfig{file=fig3.ps,width=6cm,angle=-90}
\caption{\label{fig:qegnulin}
The exponent $1/\nu$ is calculated for various values of parameters of the css regulator.
In the inset the different curves correspond to different values of
$s_2=1,5,10,20,50,100,200,500,1000$, increasing from below. There we chose $b=1$.
}
\end{figure}
\end{center}
Although the optimization of the exponent $\nu$ was performed earlier \cite{Litim:2003vp,Fischer:2006fz},
the figure shows that the css regulator can mimic such a great set of regulators
(including the Litim's, the power-law, and the exponential ones) that it can generate arbitrary values of $1/\nu$
and enables us to perform the optimization program among these fundamental regulators simultaneously.
The UV exponent is known to be $1/\nu=3$ \cite{Hamber:1991ba}, which can be easily reached and exceeded up to practically
$\infty$. It implies that one can tune the regulator parameters to get the physical value of the
exponent, e.g. according to \fig{fig:qegnulin} the choice $b=1$, $s_1=251$ and $s_2=1000$ does the trick.
We note that the strong truncation of the action in the curvature cannot give the real value, therefore
the setting up of regulator parameters in order to get $1/\nu=3$ would not leave us to
the optimal form of the regulator.
Furthermore the exponent shows high sensitivity on the parameters, which also implies that it cannot
be considered as the physical exponent, since we usually look for the parameter regime
where the physical quantities show minimal sensitivity on the regulator \cite{Litim:2001fd,Canet:2002gs}.
We cannot find any absolute extremum in the parameter space of $b$, $s_1$ and $s_2$.
However we have an inflection point at $b=1$, $s_1\to 0$ and $s_2\approx 1/2$, where we have
minimal sensitivity to the regulator parameters.
It corresponds to the Litim's regulator of the form
\beq\label{optreg}
r_{opt} = \left(\frac1{y} -\hf\right) \theta(1-y/2),\nn
\eeq
The value of the exponent there should be considered as the physical exponent, that is
$1/\nu\approx 1.472$, which is less than the half of the real value.
Other RG equations might provide better values for $\nu$ even in the framework
of the Einstein-Hilbert action \cite{Litim:2006dx,Litim:2003vp}, however in this model
if we consider this large set of regulator functions the principle of minimal sensitivity
gives this result. The power-law limit of the css regulator
does not exist for small $b$ in the case of $d=4$, because the RG equation gives
UV divergent integrals, however the finite small value of $s_2$ serves as a UV cutoff of the
loop integral. As $s_2$ gets smaller and smaller the value of $1/\nu$ decreases quickly and
tends to zero and negative values. When $1/\nu=0$ then the eigenvalue $\theta'$ in \eqn{uveig}
is purely imaginary, then the UV NGFP is not attractive any more and the trajectories evolve
towards a limit cycle as was obtained in \cite{Litim:2012vz,Bonanno:2012dg}. Further decrease of $s_2$ changes
the sign of $\theta'$ to negative making the UV NGFP a UV repulsive one.

We also calculated the position of the fixed point in the phase space. If $s_1\to\infty$
the values of $\lambda^*$ and $g^*$ become $s_2$ independent and scale according to
$\lambda^*\sim s_1^{-0.89}$ and $g^*\sim s_1^{0.89}$ implying that the product $\lambda^*g^*$
tends to a constant, $\lim_{s_1\to\infty}\lambda^*g^* \approx 0.133$. The regulator in \eqn{optreg} gives
$\lambda^*g^* = 0.136$, which value is practically the same. In the limit $s_1\to\infty$
the UV NGFP disappears, since $g^*\to\infty$ making the QEG nonrenormalizable.

\section{Crossover criticality}\label{sec:co}

There is a crossover (CO) fixed point in QEG, namely the hyperbolic GFP.
However its critical behavior is trivial since the eigenvalues corresponding to the
linearized RG flows equal the negative of the canonical dimension of the
couplings, and the inclusion of further couplings or taking into account
corrections from the gradient expansion beyond the local potential approximation (LPA)
do not change their values. In order to test the CO scaling criticality we treat the 3-dimensional
$O(1)$ model around the Wilson-Fisher (WF) fixed point and the corresponding critical
exponent $\nu$. We start with the effective action
\beq\label{o3d1effac}
\t V = \sum_{i=1}^N \frac{ \t g_i}{(2i)!}\phi^{2i},
\eeq
with the dimensionless couplings $\t g_i$ and restrict ourselves to the truncations $N=2$ and $N=4$.
In LPA the evolution equations for the first two couplings read as
\bea
\dot{\t g}_1 &=& -2 \t g_1 + \t g_2 \bar \Phi^2_{3/2}(\t g_1), \nn
\dot{\t g}_2 &=& -\t g_2 + 6 \t g_2^2 \bar \Phi^3_{3/2}(\t g_1),
\eea
 where the threshold function is introduced according to
\beq\label{thrs2}
\bar \Phi^p_n(\omega) = \frac{1}{(4\pi)^n\Gamma(n)}
\int_0^\infty dy y^{n+1}\frac{r'}{(y(1+r)+\omega)^p}.
\eeq
The evolution equations also contain regulator-dependence. As it is well-known the model plays
the role of the testing ground for any new improvement of renormalization. The exponent $\nu$
of the WF fixed point is well-known, and the optimization of the RG regulator
was first performed around the WF fixed point giving the Litim's regulator, furthermore the value of $\nu$
was calculated for every available type of regulators in the literature.
However by using the css regulator we can get an optimized exponent among
fundamental regulators which was not investigated so far. We considered the RG
evolution for 2 and 4 couplings. The truncation of the scalar potential
to only two terms in the Taylor expansion is very strong, but
we also have only two couplings in QEG with the same level of truncation.
The values of the exponents are plotted in \fig{fig:d3o1nu}.
\begin{center}
\begin{figure}[ht]
\epsfig{file=fig4.ps,width=6cm,angle=-90}
\caption{\label{fig:d3o1nu}
The exponent $\nu$ is calculated for 2 and 4 couplings and
shown by the curves at the bottom and at the top, respectively. The Litim's regulator gives the
extremal exponent in both cases at the point of the crossing of the grey lines.
}
\end{figure}
\end{center}
For 2 couplings we obtained that we have a maximum of the exponent $\nu$ at
$b=s_2=1$ and $s_1\to 0$ which corresponds to the Litim's regulator in \eqn{limregs}.
There the sensitivity of $\nu$ on the regulator parameters is minimal.
Fortunately even the truncation with $N=2$ can provide us a real physical exponent
although it is smaller than the proper value of $\nu$.
We repeated the calculations for $N= 4$ couplings. We got another extremum but now it is
a minimum and it corresponds to the Litim's regulator, too. Further increase of the
number of the couplings and the inclusion of the evolving anomalous dimension $\eta$
may give the Litim's regulator, too.

\section{Infrared criticality}\label{sec:ir}

In the broken phase of field theoretical models we usually find an IR fixed point.
This is the case in QEG, too \cite{Donkin:2012ud,Christiansen:2012rx,Nagy:2012rn}.
We calculated the exponent $\nu$ by the dynamically induced correlation
\cite{Nagy:2012rn,Nagy:2012np,Nagy:2012qz} around the IR fixed point.
We obtained that the exponent $\nu$ around the IR fixed point equals the one
calculated around the GFP and it is $\nu=1/2$. We also determined the scaling of
the anomalous dimension $\eta$, the results are plotted in \fig{fig:qegeta}.
\begin{center}
\begin{figure}[ht]
\epsfig{file=fig5.ps,width=6cm,angle=-90}
\caption{\label{fig:qegeta}
The anomalous dimension $\eta$ is plotted as the function of the scale. In the vicinity
of the fixed points the $\eta$ shows different scaling behavior. We chose $s_1=s_2=b=1$
for the css regulator.
}
\end{figure}
\end{center}
From the evolution of the Newton constant $g$ it is trivial that in the UV NGFP we have
$\eta=-2$. This corresponds to marginal scaling of $\eta$ in \fig{fig:qegeta}. Then
as the flow approaches the GFP it tends to zero as a power law function $k^2$, where
the exponent 2 is universal.
Going further it starts to diverge as $(k-k_c)^{-\alpha}$, where $\alpha$ is around 3/2 and
shows a slight parameter dependence. The value $k_c$ denotes
the critical scale where the evolution stops due to the appearing singularity.
From the threshold function in \eqn{thrs} in the RG equations the singularity condition is
\beq
y(1+r)-2\lambda = 0,
\eeq
e.g. the $\lambda^*=1/2$, when $r$ is the Litim's regulator. The exponents are summarized
in \tab{tab:exps}.
\begin{table}
\begin{center}
\begin{tabular}{|c||c|c|c|}
\hline
exponent & {\bf UV} & {\bf G} & {\bf IR} \\
\hline
\hline
$\nu$ & 1.472 & 1/2 & 1/2 \\
\hline
$\eta$ & -2 & $k^2$ & $k^{-7/4}$ \\
\hline
\end{tabular}\caption{\label{tab:exps}
The summary of the critical exponents $\nu$ and the anomalous dimensions $\eta$ at the  various
fixed points of QEG.
}
\end{center}
\end{table}
The IR scaling of $\eta$ with the exponent $7/4$ belongs to the Litim's regulator.

\section{Summary}\label{sec:sum}

By using the functional RG method the critical exponent $\nu$
of the correlation length and the anomalous dimension $\eta$ are calculated for the QEG
with the css regulator. We considered the Einstein-Hilbert
action with the Newton constant $g$ and the cosmological constant $\lambda$.
The value of $\nu$ is calculated around the UV NGFP, the GFP, and the IR fixed point.
We obtained that $\nu$ can be arbitrary around the UV NGFP. There can be several reasons
for this result. On the one hand although the regulator is a mass-like term, it possesses a complicated
momentum dependence which might introduce non-local interactions into the action.
Furthermore the deep IR physics cannot be altered by the regulator, but the UV limit can show strong
regulator dependence. On the other hand the Einstein-Hilbert action contains only
two couplings. The inclusion of further couplings in the Taylor expansion in the curvature
\cite{Lauscher:2002sq,Rechenberger:2012pm,Bonanno:2013dja}
may restrict the value of $\nu$ to a certain interval. The several possible extensions
\cite{Eichhorn:2013ug,Henz:2013oxa} can
also give some restriction to the value of $\nu$. We showed that the exponent $\nu=1.472$
can be got if we look for that value which has minimal sensitivity on the regulator
parameters at the Litim's regulator.

The CO scaling around the GFP gives $\nu=1/2$ independently on the approximations
that were used. We demonstrated, that the CO scaling in the 3-dimensional $O(1)$ model
can give a wide range of the value of $\nu$ however an extremum appears in the parameter
space, which corresponds to the Litim's regulator, too.
The IR fixed point gives $\nu=1/2$ which is inherited from the CO GFP.

We also showed that the anomalous dimension takes the value $\eta=-2$ at the UV NGFP, and
it scales in an irrelevant manner by a universal exponent when the RG flow approaches the GFP,
while it is relevant in the IR, with an exponent which shows moderate dependence on the regulator parameters.

\section*{Acknowledgments}

This research was realized in the frames of T\'AMOP 4.2.4.A/2-11-1-2012-0001 ``National Excellence
Program – Elaborating and operating an inland student and researcher personal support
system'' The project was subsidized by the European Union and co-financed by the
European Social Fund.