# Growing Brains: Co-emergence of Anatomical and Functional Modularity in
Recurrent Neural Networks

Ziming Liu
MIT & IAIFI
<EMAIL_ADDRESS>
Mikail Khona∗
MIT
<EMAIL_ADDRESS>
Ila R. Fiete
MIT
<EMAIL_ADDRESS>
Max Tegmark
MIT & IAIFI
<EMAIL_ADDRESS>
Equal contribution

###### Abstract

Recurrent neural networks (RNNs) trained on compositional tasks can exhibit
functional modularity [1, 2], in which neurons can be clustered by activity
similarity and participation in shared computational subtasks. Unlike brains,
these RNNs do not exhibit _anatomical modularity_ , in which functional
clustering is correlated with strong recurrent coupling and spatial
localization of functional clusters. Contrasting with functional modularity,
which can be ephemerally dependent on the input [2], anatomically modular
networks form a robust substrate for solving the same subtasks in the future.
To examine whether it is possible to grow brain-like anatomical modularity, we
apply a recent machine learning method, brain-inspired modular training
(BIMT), to a network being trained to solve a set of compositional cognitive
tasks. We find that functional and anatomical clustering emerge together, such
that functionally similar neurons also become spatially localized and
interconnected. Moreover, compared to standard $L_{1}$ or no regularization
settings, the model exhibits superior performance by optimally balancing task
performance and network sparsity. In addition to achieving brain-like
organization in RNNs, our findings also suggest that BIMT holds promise for
applications in neuromorphic computing and enhancing the interpretability of
neural network architectures.

##  1 Introduction

A powerful way for networks to generalize is through modularity: If seen and
unseen tasks in the world consist of combinations of subtasks, then a new task
can be quickly solved by decomposing it into the set of previously seen
subtasks, and tackling those based on prior learning. Recent work shows that
RNNs trained on a set of tasks drawn by combining subtasks from a common
dictionary begin to exhibit functional modularity, with similar activity
profiles across neurons responding to the same subtask. However, the formed
clusters were not anatomical. Anatomical clustering with localization of
function is a central feature of brains [3]: for example, visual processing
for object recognition is localized to the ventral visual pathway while the
initiation of voluntary movements is confined to a few motor and premotor
cortical regions. Anatomical modularization can facilitate continual learning:
if it takes the form of spatial localization, then new inputs can be easily
routed to the module; if it takes the form of recurrent connectivity, it
provides lasting substructures for solving specific computations on future
tasks. By contrast, functional clustering alone can be ephemeral, with
groupings that might be defined primarily by correlations in the inputs [2].
When the input correlations change, functional modules could disappear.

Recently, the method of brain-inspired modular training (BIMT) was proposed as
a way to make artificial neural networks modular and more interpretable [4].
The key idea of BIMT is to encourage local neural connections via two
optimization terms: distance-dependent weight regularization and discrete
neuron swapping. Here we ask whether BIMT can answer a fundamental question
about neuroscience: Can spatial constraints and wiring costs [5, 6, 7]
together lead to the emergence of anatomical modules that are also
functionally distinct?

We study how BIMT can lead to the emergence of spatial modules in a multitask
learning setting relevant to cognitive systems neuroscience with two sets of
combinatorially constructed tasks, 20-Cog-tasks [1] and Mod-Cog-tasks [2]. We
train recurrent neural networks (RNNs) on these tasks with BIMT in the
supervised setup. We observe brain-like spatial organization emerging in the
hidden layer of the RNN: neurons that are functionally similar are also
localized in space (Figure 1c). Such locality and sparsity are gained with no
sacrifice in performance, or even accompanied by an improvement in
performance. We introduce our methods in Section 2, and present results in
Section 3. Due to limited space, the main paper focus on 20-Cog-tasks [1] and
leave the results for Mod-Cog-tasks [2] in Appendix B.

##  2 Method

Figure 1: Training an RNN with BIMT on cognitive tasks. a: Visualization of
the network. Each line represents a weight; blue/red means positive/negative
weights; thickness corresponds to magnitudes. b: The hidden neurons are
clustered into functional modules. c: These functional modules (distinguished
by colors) are also clustered in space, visually resembling a brain. By
contrast, $L_{1}$ regularization leads to no anatomical modules. The network
has (d) good performance, (e) high sparsity and good locality (f). g and h:
Trade-off between performance (error) and sparsity. g uses the number of
active neurons as the sparsity measure, while h uses the wiring length as the
sparsity measure. In both cases, the Pareto frontier of BIMT is better than
that of $L_{1}$ regularization and no regularization.

Brain-inspired modular training (BIMT) Biological neural networks (e.g.,
brains) differ from artificial ones in that biological ones restrain neuronal
connections to be local in space, leading to anatomical modularity. Motivated
by this observation, [4] proposed brain-inspired modular training (BIMT) to
facilitate modularity and interpretability of artificial neural networks. The
idea is to embed neurons into a geometric space and minimize the total
connection cost by adding the connection cost as a penalty to the loss
function and swapping neurons if necessary. On a number of math and machine
learning datasets, they show that BIMT is able to put functionally relevant
neurons close to each other in space, just like brains. It is thus natural to
ask: Can BIMT give something back to neuroscience? In this work, we apply BIMT
to recurrent neural networks (RNN) for cognitive tasks, and show that the
neurons in the hidden layer are organized into modules which are both
anatomically and functionally distinct, just like brains.

RNN We take a simple recurrent neural network (RNN) in the context of systems
neuroscience, which is defined by

$\displaystyle\mathbf{h}_{t+1}=\phi(\mathbf{W}\mathbf{h}_{t}+\mathbf{W}_{\rm
in}\mathbf{u}_{t}+\mathbf{b}^{h}),$ (1)
$\displaystyle\mathbf{o}_{t+1}=\mathbf{W}_{\rm
out}\mathbf{h}_{t+1}+\mathbf{b}^{o},$

where $\mathbf{u}_{t}\in\mathbb{R}^{n_{u}}$,
$\mathbf{h}_{t}\in\mathbb{R}^{n_{h}\times n_{h}}$,
$\mathbf{o}_{t}\in\mathbb{R}^{n_{o}}$, $\mathbf{W}\in\mathbb{R}^{n_{h}\times
n_{h}}\times\mathbb{R}^{n_{h}\times n_{h}}$. We place $n_{h}\times n_{h}$
hidden neurons uniformly on a 2D grid $[0,1]^{2}$ (see Figure 1a), so the $ij$
neuron (the neuron in the $i^{\rm th}$ row and the $j^{\rm th}$ column) is
located at $(i/n_{h},j/n_{h})$. The $L_{1}$ distance between the $ij$ neuron
and the $mn$ neuron is thus $\mathbf{D}_{ij,mn}\equiv(|i-m|+|j-n|)/n_{h}$
111One can choose other distances, e.g., $L_{2}$ distance. We choose $L_{1}$
distance because it is more consistent with $L_{1}$ regularization.. We define
the RNN’s connection cost as

$\ell_{cc}=\underbrace{\|\mathbf{W}\|_{1}+\|\mathbf{W}_{\rm
in}\|_{1}+\|\mathbf{W}_{\rm
out}\|_{1}+|\mathbf{b}^{h}|_{1}+|\mathbf{b}^{o}|_{1}}_{\text{vanilla $L_{1}$
regularization}}+\underbrace{A\|\mathbf{D}\odot\mathbf{W}\|_{1}}_{\text{distance-
aware regularization}},$ (2)

where $\|\mathbf{M}\|_{1}\equiv\sum_{ij}|M_{ij}|$ and
$|\mathbf{v}|_{1}\equiv\sum_{i}|v_{i}|$ are matrix and vector $L_{1}$-norm,
respectively. $A$ is a hyper-parameter controling the strength of locality
constraint.

Cognitive tasks The 20-Cog-tasks are a set of simple cognitive tasks inspired
by experiments with rodents and non-human primates performed by systems
neuroscientists [1]. These tasks are designed to fall into families where each
family is defined by a set of computations drawn from a common pool of
computational primitives. Thus, the tasks have shared subtasks and an optimal
solution is to form clusters of neurons specialized to these subtasks and
share them across tasks, illustrated in Figure 2a. In our experiments, we set
both input rings to have 16 dimensions, so input dimension $n_{u}=16\times
2+20+1=53$ (where 20 for the one-hot length-20 task vector, and 1 for
fixation). Hidden neurons are aranged as $20\times 20$ grid ($n_{h}=20$), and
output dimension $n_{o}=16+1=17$. The prediction loss $\ell_{\rm pred}$ is the
cross-entropy between the ground truth and the predicted reaction.

BIMT loss simply combines the prediction loss and the connection cost, i.e.,
the total loss function is

$\ell=\ell_{\rm pred}+\lambda\ell_{cc},$ (3)

where $\lambda\geq 0$ is the strength of penalizing connection costs. When
$\lambda=0$, it boils down to train a fully-connected RNN without sparsity
constraint; when $\lambda>0$ but $A=0$, it boils down to train with vanilla
$L_{1}$ regularization. Besides adding connection costs as regularization,
BIMT allows swapping neurons to further reduce $\ell_{cc}$ by avoiding local
minima, e.g., if a neural network is initialized to be performing well but has
non-local connections, without swapping, the network would not change much and
maintain those non-local connections.

##  3 Results

We focus first on results in the 20-Cog-tasks ; qualitatively similar results
on the Mod-Cog-tasks are included in Appendix B.

###  3.1 BIMT learns a 2D brain that solves all 20-Cog-tasks

We show results for $A=1.0$ and $\lambda=10^{-5}$ in Figure 1. a shows the
connectivity graph of the RNN. All the weights are plotted as lines whose
thicknesses are proportional to their magnitudes 222Some weights are visually
vanishing, because their magnitudes are too small, although we indeed plot
them., and blue/red means positive/negative weights. BIMT learns to prune away
peripheral neurons and concentrate important neurons only in the middle.
Following [1], we cluster neurons into functional modules based on their
normalized task variance (shown in b). These functional modules, colored by
different colors, are shown to be also anatomically modular, i.e., spatially
local (shown in c), visually resembling a brain. By contrast, $L_{1}$
regularization does not induce any anatomical module. d shows that the network
performs reasonably well. The network is sparse; e shows that it only contains
around 100 important neurons (measured by sum of task variances). The
connections in the hidden layer are mostly local, as we hoped; f shows that
weights decay fast as distance increases, which is similar to fixed local-
masked RNN [2]. However intriguingly, there is a second peak of (relatively)
strong connections around distances being 0.6, which is probably attributed to
inter-module connections.

###  3.2 Sparsity vs Accuracy Tradeoff

There is a Pareto frontier showing the trade-off between sparsity and
accuracy, shown in Figure 1g and h. We also compute the trade-off for networks
with vanilla $L_{1}$ regularization. We use two sparsity measures: the number
of active neurons (a neuron is active if its sum of tasks variances is larger
than $10^{-3}$), and the wiring length (sum of lengths of all active
connections; a connection is active if its weight magnitude is larger than
$10^{-2}$). Under both measures, BIMT is superior than $L_{1}$ regularization
in terms of having better Pareto frontier.

###  3.3 Anatomical modularity

Figure 2: (a) The myriad ways anatomical and functional modularity can present
itself in trained RNNs. (b) We test anatomical modularity for neural networks
with BIMT (left), L1 regularization (middle) or no regularization (right). We
propose two metrics, fraction of isolated neurons (top) and average
(functional) cluster size (bottom) to measure anatomical modularity. For both
metrics, smaller is better. We compare the trained network with networks whose
useful hidden neurons are randomly shuffled. No regularization and L1
regularization are in the distributions of randomly shuffled networks, while
BIMT is significantly out of distribution (smaller) than random ones,
indicating anatomical modularity.

Anatomical modularity means that neurons with similar functions are placed
close to each other in space. Because neurons of fully-connected layers have
permutation symmetries, there is no incentive for them to develop anatomical
modularity. By contrast, since BIMT penalizes connection costs, BIMT networks
potentially have anatomical modularity. In Figure 1c, each neuron’s functional
cluster is marked and there are clear spatial clusters in which all neurons
belong to the same functional cluster. Quantitatively, we propose two metrics:
(1) the fraction of isolated neurons. A neuron is isolated if none of its
(eight) neighbors belongs to the same functional cluster. (2) the average size
of functional clusters. For both metrics, the smaller the better. For
baselines, we randomly shuffle important hidden neurons 333A neuron is
important if the sum of its task variances is above $10^{-3}$.. Since
different random shuffling may yield different results, we try 10000 different
random seeds and plot histograms for the metrics. We compute the two metrics
for networks trained with BIMT, $L_{1}$ regularization, or no regularization
in Figure 2. Only BIMT networks are seen to be significantly out-of-
distribution from baselines, implying anatomical modularity. In future work,
we hope to explore improvements of BIMT to further increase the functional
modularity of the “brain" seen in Figure 1c.

Our main findings remain qualitatively similar when (i) the topology of hidden
layer is changed or (ii) the tasks are significantly harder, so that they
require more-involved recurrent connectivity within the network. Specifically,
in Appendix A, we include the results for 1D hidden layer. In Appendix B, we
find that results on a more-complex set of 84 cognitive tasks, the Mod-Cog-
tasks [2], show qualitatively similar anatomical clustering results,
demonstrating the robustness and generality of our core findings.

## Acknowledgement

ZL and MT are supported by IAIFI through NSF grant PHY-2019786, the
Foundational Questions Institute and the Rothberg Family Fund for Cognitive
Science. IRF is supported by the Simons Foundation through the Simons
Collaboration on the Global Brain, the ONR, the Howard Hughes Medical
Institute through the Faculty Scholars Program and the K. Lisa Yang ICoN
Center. MK acknowledges funding from the Department of Physics, MIT.

## References

  * Yang et al. [2019] Guangyu Robert Yang, Madhura R Joglekar, H Francis Song, William T Newsome, and Xiao-Jing Wang.  Task representations in neural networks trained to perform many cognitive tasks.  _Nature neuroscience_ , 22(2):297–306, 2019.
  * Khona et al. [2023] Mikail Khona, Sarthak Chandra, Joy J Ma, and Ila R Fiete.  Winning the lottery with neural connectivity constraints: Faster learning across cognitive tasks with spatially constrained sparse rnns.  _Neural Computation_ , pages 1–20, 2023.
  * Ferrier [1874] David Ferrier.  The localization of function in the brain.  _Proceedings of the Royal Society of London_ , 22(148-155):228–232, 1874.
  * Liu et al. [2023] Ziming Liu, Eric Gan, and Max Tegmark.  Seeing is believing: Brain-inspired modular training for mechanistic interpretability.  _arXiv preprint arXiv:2305.08746_ , 2023.
  * Chen et al. [2006] Beth L Chen, David H Hall, and Dmitri B Chklovskii.  Wiring optimization can relate neuronal structure and function.  _Proceedings of the National Academy of Sciences_ , 103(12):4723–4728, 2006.
  * Chklovskii et al. [2002] Dmitri B Chklovskii, Thomas Schikorski, and Charles F Stevens.  Wiring optimization in cortical circuits.  _Neuron_ , 34(3):341–347, 2002.
  * Chklovskii and Koulakov [2004] Dmitri B Chklovskii and Alexei A Koulakov.  Maps in the brain: what can we learn from them?  _Annu. Rev. Neurosci._ , 27:369–392, 2004.

Appendix

##  Appendix A 1D results for 20-Cog-tasks

In the main paper, we lay hidden neurons on a 2D grid, which is mainly because
the neocortex is effectively a two-dimensional sheet. Algorithm-wise, BIMT
allows any geometric space. For example, we can arrange hidden neurons into a
one-dimensional grid. Results are included in Figure 3,4 and 5.

Figure 3: Training an RNN with BIMT on cognitive tasks. a: Visualization of
the network. Each line represents a weight; blue/red means positive/negative
weights; thickness corresponds to magnitudes. The hidden-to-hidden connections
are zoomed-in in b, and displayed as a weight matrix in c, which is strongly
diagonal. d and e: We decompose important neurons into clusters based on task
variance [1]. The BIMT network has reasonable performance (f), high sparsity
(g) and high locality (h).

(a)

(b)

(c)

Figure 4: Sparsity-accuracy tradeoff and spatial structures for neural
networks. (a) BIMT has a similar Pareto frontier to the baseline L1
regularization, while fully-connected network is neither sparse nor optimally
performant. (b) Networks trained with $L_{1}$ penalty do not have spatial
structures, while (c) BIMT networks learn to prune away peripheral hidden
neurons, and have anatomical modularity (see Figure 5). Figure 5: (a) The
myriad ways anatomical and functional modularity can present itself in trained
RNNs. (b) We test anatomical modularity for neural networks with vanilla
training (left), L1 penalty (middle) or BIMT (right). We propose two metrics,
ratio of isolated neurons (top) and variance of indices within clusters
(bottom) to measure anatomical modularity. For both metrics, smaller is
better. We compare the trained network with networks whose hidden neurons are
randomly shuffled. Fully-connected networks and L1 networks are in the
distributions of randomly shuffled networks, while BIMT is significantly out
of distribution (smaller) than random ones, indicating anatomical modularity.

##  Appendix B Results for 84 Mod-Cog-tasks

The 20-Cog-tasks are a relatively simple set that can be solved with only
autapses in the hidden layers [2], raising the question of whether, when
subtasks require extensive recurrent connectivity in the hidden networks, BIMT
will lead to emergence of anatomical modularity. Khona et al. [2] extended
Yang 20 tasks [1] to a more diverse set of tasks (84 in total) that require
recurrent connectivity among multiple neurons within each module. We find that
the conclusions in the main paper also hold for this augmented suite of
cognitive tasks, demonstrating the robustness of our claims. Results are
included in Figure 6.

Figure 6: Results for 84 Mod-Cog tasks. a: the visualization of the
connectivity graph of the 1D RNN; b: neurons are grouped into functional
modules based on task variances; c: the network has high performance for most
tasks; d: the trade-off between performance and sparsity. The pareto frontier
of BIMT is similar to that of $L_{1}$, which means that anatomical modularity
can be achieved without sacrifice in performance or sparsity. e and f: RNNs
trained with BIMT clearly demonstrate anatomical modularity since the two
metrics (fraction of isolated neurons and variance of indices within clusters)
are significantly out-of-distribution of random baselines.


# Universal properties of the Kardar-Parisi-Zhang equation with quenched
columnar disorders

Astik Haldar<EMAIL_ADDRESS>Condensed Matter Physics Division, Saha
Institute of Nuclear Physics, HBNI, Calcutta 700064, India

###### Abstract

Inspired by the recent results on totally asymmetric simple exclusion
processes on a periodic lattice with short-ranged quenched hopping rates [A.
Haldar, A. Basu, Phys Rev Research 2, 043073 (2020)], we study the universal
scaling properties of the Kardar-Parisi-Zhang (KPZ) equation with short-ranged
quenched columnar disorder in general $d$-dimensions. We show that there are
generic propagating modes in the system that have their origin in the quenched
disorder and make the system anisotropic. We argue that the presence of the
propagating modes actually make the effects of the quenched disorder
irrelevant, making the universal long wavelength scaling property belong to
the well-known KPZ universality class. On the other hand, when these waves
vanish in a special limit of the model, new universality class emerges with
dimension $d=4$ as the lower critical dimension, above which the system is
speculated to admit a disorder-induced roughening transition to a
perturbatively inaccessible rough phase.

##  I Introduction

Quenched disorders are known to significantly affect the universal scaling
properties of condensed matter systems and statistical models in equilibrium.
For example, the universal critical properties of classical $N$-component spin
models near its second order phase transitions can be affected by the quenched
disorder, leading to a new universality class Lubensky (1975); Mukherjee and
Basu (2020). On the other hand, quenched random fields can destroy the ordered
phase altogether Grinstein and Luther (1976). Scaling properties of interface
growth models are also found to be affected by quenched disorders Jeong et al.
(1996); Csahók et al. (1993); Barabási et al. (1996); Amaral et al. (1994).
Understanding of how quenched disorders affect universal properties in
nonequilibrium systems is far less understood than their equilibrium
counterparts, and remains a challenging problem. Lack of any general framework
to study the statistical properties of nonequilibrium systems has often
prompted physicists to construct and study simple conceptual models that are
easily analytically tractable, such that questions of basic principles can be
answered within well-controlled approximations.

Recently, the universal scaling properties of the density fluctuations in a
totally asymmetric simple exclusion process (TASEP) in a ring geometry with
quenched disordered hopping rates has been studied by mapping the problem into
a one-dimensional ($1d$) Kardar-Parisi-Zhang equation with columnar disorder
Haldar and Basu (2020). It has been shown that away from the half-filling, the
quenched disorder is irrelevant, and the universal long wavelength properties
of the system is identical to that of the $1d$ KPZ equation without any
disorder. On the other hand, close to half-filling, quenched disorder is
relevant, a new universality class different from the $1d$ KPZ universality
emerges, leading to the scaling exponents taking values different from their
$1d$ KPZ counterpart. These results qualitatively agree with the simulations
performed Tripathy and Barma (1997, 1998); de Queiroz and Stinchcombe (2008)
for driven lattice gases in presence of quenched disorder .

In this article we generalise the $1d$ KPZ equation with short-ranged quenched
columnar disorder Haldar and Basu (2020) to $d$-dimensions, and study its
universal scaling properties in the long wavelength limit. Our main results in
the article are as follows. (i) The generalised $d$-dimensional equation does
not remain invariant under Galilean transformation in presence of the
disorder, and there are generic underdamped propagating waves in the system.
The latter makes the system necessarily ansiotropic. (ii) Surprisingly, the
presence of the waves makes the disorder irrelevant (in a renormalisation
group or RG sense). As a result, the Galilean invariance and isotropy are
restored in the long wavelength limit; the universal scaling property
naturally belongs to the $d$-dimensional KPZ universality class. (iii) For
specific choices of the model parameters, the waves vanish and the system is
isotropic. In this case the disorder is relevant (in a RG sense), and the
Galilean invariance remains broken in the long wavelength limit. Consequently,
the universal scaling exponents belong to a new universality class different
from the pure KPZ equation. We argue that $d=4$ is the lower critical
dimension of the model in this case, such that above $d=4$, the model can
undergo a roughening transition to a (perturbatively inaccessible) rough
phase. We use a one-loop dynamic RG calculational framework for our work.

The rest of this article is organised in the following manner. We set up the
hydrodynamic equations for height of growing surfaces with quenched columnar
disorder in Sec. II. The KPZ universality has been discussed in Sec. III. We
discuss the scaling behaviour in linear theory in Sec. IV.1, and then
calculate the nonlinear effects on the scaling exponents by using RG framework
for this quenched disordered model in rest of Sec. IV. We summarise our
results in Sec. V. Some technical details are provided in Appendix for the
interested readers.

##  II KPZ equation with columnar quenched disorder

We begin by recalling the $1d$ hydrodynamic equation, derived in Ref. Haldar
and Basu (2020) for density fluctuations in a periodic TASEP with quenched
disordered hopping rates in the thermodynamic limit. The equation of motion
(EoM) for density fluctuation $\phi(x,t)$ reads Haldar and Basu (2020)

$\displaystyle\frac{\partial\phi(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle\nu\frac{\partial^{2}\phi(x,t)}{\partial
x^{2}}+\lambda_{1}\frac{\partial\phi(x,t)}{\partial
x}+\lambda_{m}\frac{\partial\delta m(x)}{\partial x}$ (1)
$\displaystyle-\frac{\lambda}{2}\frac{\partial\phi^{2}(x,t)}{\partial
x}+\lambda_{2}\frac{\partial}{\partial x}[\delta m(x)\phi(x,t)]$
$\displaystyle-\lambda_{3}\frac{\partial}{\partial x}[\delta
m(x)\phi^{2}(x,t)]+\frac{\partial f(x,t)}{\partial x}.$

Here, $\nu$ is a diffusion coefficient, $\lambda$, $\lambda_{m}$,
$\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ are different model parameters. In
particular, $\lambda_{1}$, $\lambda_{2}$ are both proportional to
$(2\rho_{0}-1)$ where $\rho_{0}$ is the mean particle density. Hence,
$\lambda_{1}$, $\lambda_{2}$ vanish in the half-filled limit ($\rho_{0}=1/2$)
whereas those have non zero value for away from the half-filled limit Haldar
and Basu (2020). Furthermore, $m(x)$ is the random quenched disorder that has
its origin in the quenched disorder hopping rates in the underlying lattice-
gas model. As explained in Ref. Haldar and Basu (2020), density fluctuation
$\phi(x,t)$ is driven by two conserved stochastic variables: quenched noise
$\partial_{x}\delta m(x)$ and annealed noise $\partial_{x}f(x,t)$. We write
$m(x)=m_{0}+\delta m(x)$, where $m_{0}=\overline{m(x)}$ is the mean value of
$m(x)$, and $\delta m(x)$ is local fluctuation of $m(x)$ about $m_{0}$. We
assume $\delta m(x)<m_{0}$ to ensure $m(x)$ remains positive everywhere; we
choose $\delta m(x)$ to be short-ranged Gaussian distributed with a variance

$\overline{\delta m(x)\delta m(0)}=2\tilde{D}\delta(x);$ (2)

$\tilde{D}$ is the strength of the quenched noise, and is definite positive.

Stochastic function $f(x,t)$ is an annealed noise which is assumed to be zero-
mean and Gaussian-distributed with a variance

$\langle f(x,t)f(0,0)\rangle=2D\delta(x)\delta(t).$ (3)

Here, $D>0$ is the noise strength, the analogue of temperature in
nonequilibrium systems. Here $\langle...\rangle$ implies the averages over
time-dependent noise distribution and an overline implies the averages over
quenched disorder distribution.

The variances of $f(x,t)$ and $\delta m(x)$ in the Fourier space are

$\displaystyle\langle
f(k_{1},\omega_{1})f(k_{2},\omega_{2})\rangle=2D\delta(k_{1}+k_{2})\delta(\omega_{1}+\omega_{2}),$
(4a) $\displaystyle\overline{\delta m(k_{1},\omega_{1})\delta
m(k_{2},\omega_{2})}=2\tilde{D}\delta(k_{1}+k_{2})\delta(\omega_{1}+\omega_{2})\delta(\omega_{1}).$
(4b)

Let us introduce a height variable via the transformation
$\phi(x,t)=\partial_{x}h$. This transforms (1) to

$\displaystyle\frac{\partial h(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle\nu\partial_{xx}h(x,t)+\lambda_{1}\partial_{x}h(x,t)+\lambda_{m}\delta
m(x)$ (5)
$\displaystyle-\frac{\lambda}{2}[\partial_{x}h(x,t)]^{2}+\lambda_{2}[\delta
m(x)\partial_{x}h(x,t)]$ $\displaystyle-\lambda_{3}[\delta
m(x)(\partial_{x}h(x,t))^{2}]+f(x,t).$

Here $h(x,t)$, a single-valued function, is the height of the $1d$ surface
with respect to an arbitrary baseline. Equation (5) gives the dynamics of a
growing $1d$ surface with quenched columnar disorder. When
$\lambda_{1}=0=\lambda_{2}$, Eq. (5) is invariant under inversion of space,
i.e., under $x\rightarrow-x$, whereas for $\lambda_{1}\neq 0,\,\lambda_{2}\neq
0$, it is not due to the presence of the propagating waves. These correspond,
respectively, to near the half-filled and away from the half-filled limits for
density fluctuation in Eq. (1).

We now generalise the problem to a $d$-dimensional growing surface in presence
of columnar quenched disorder. This can be done in several ways. We generalise
in a way to ensure that the $d$-dimensional generalised equation has a
definite symmetry, and reduces to the $1d$ case in a simple and unambiguous
way. Noticing that there are two terms (with coefficients $\lambda_{1}$,
$\lambda_{2}$) containing single spatial derivative on $h$ in (5) which imply
the absence of the space inversion symmetry in EoM, we generalise the EoM of
$h({\bf x},t)$ to $d$-dimensions by assuming breakdown of the inversion
symmetry along one particular dimension, while the system remains inversion-
symmetric along all other $(d-1)$ directions. We use $\parallel$ as subscript
to denote the special direction in which the inversion symmetry is absent. The
generalised EoM for $h({\bf x},t)$ ($\bf x$ now refers as position on
$d$-dimensional surface) has the generic form

$\displaystyle\frac{\partial h}{\partial t}$ $\displaystyle=$
$\displaystyle\nu\nabla^{2}h+\lambda_{1}\partial_{\parallel}h+\lambda_{m}\delta
m-\frac{\lambda}{2}[\bm{\nabla}h]^{2}$ (6) $\displaystyle+\lambda_{2}[\delta
m\partial_{\parallel}h]-\lambda_{3}[\delta m(\bm{\nabla}h)^{2}]+f.$

It is clear that Eq. (6) has a lower symmetry than the $d$-dimensional KPZ
equation. It is rotationally invariant in the $(d-1)$-dimensional subspace,
that excludes the “$\parallel$”-direction, unlike the $d$-dimensional KPZ
equation that is rotationally invariant in the full $d$-dimensional space. It
is instructive to write down the “Burgers” version of Eq. (6) for a “quenched
disorder Burgers velocity” ${\bf v}$ by using the transformation ${\bf
v}={\bm{\nabla}}h$:

$\displaystyle\frac{\partial{\bf v}}{\partial t}$ $\displaystyle=$
$\displaystyle\nu\nabla^{2}{\bf
v}+\lambda_{1}\bm{\nabla}v_{\parallel}+\lambda_{m}\bm{\nabla}\delta
m-\frac{\lambda}{2}\bm{\nabla}[v]^{2}$ (7)
$\displaystyle+\lambda_{2}\bm{\nabla}[\delta
mv_{\parallel}]-\lambda_{3}\bm{\nabla}[\delta m(v)^{2}]+\bm{\nabla}f.$

The variances of the noises in $d$-dimensional space are just the direct
generalisations, respectively, of (2) and (3):

$\displaystyle\overline{\delta m({\bf x})\delta
m(0)}=2\tilde{D}\delta^{d}({\bf x}).$ (8a) $\displaystyle\langle f({\bf
x},t)f(0,0)\rangle=2D\delta^{d}({\bf x})\delta(t).$ (8b)

Due to the presence of quenched disorders, Eq. (6) does not remain invariant
under the Galilean transformation ${\bf x}\mapsto{\bf x}+{\bf
v}_{0}t,\,t\mapsto t,\,h({\bf x},t\mapsto h({\bf x}+{\bf v}_{0}t)+{\bf
v}_{0}\cdot{\bf x})$. In Eq. (6) the term with linear derivative
$\lambda_{1}\partial_{\parallel}h$ implies the existence of underdamped
propagating modes. The surface would be anisotropic due to presence of the
propagating modes with picking up a specific direction, that without any loss
of generality is assumed to be along the $x$-direction, i.e., the “parallel
direction” $x_{\parallel}$ is identified with the $x$-direction. In the
absence of any quenched disorder $(\delta m({\bf x})=0)$, Eq. (10) reduces to
the well-known KPZ equation.

In what follows below, it is convenient to write $h({\bf x},t)=h_{1}({\bf
x})+h_{2}({\bf x},t)$ Haldar and Basu (2020). Time-independent function
$h_{1}({\bf x})$ then satisfies

$\displaystyle-\nu_{\psi}\nabla^{2}h_{1}-\lambda_{1\psi}\partial_{\parallel}h_{1}-\lambda_{m}\delta
m+\frac{\lambda_{\psi}}{2}(\bm{\nabla}h_{1})^{2}$
$\displaystyle-\lambda_{2\psi}[\delta
m(\partial_{\parallel}h_{1})]+\lambda_{3\psi}[\delta
m(\bm{\nabla}h_{1})^{2}]=0.$ (9)

Time-dependent function $h_{2}({\bf x},t)$ satisfies

$\displaystyle\frac{\partial h_{2}}{\partial t}$ $\displaystyle=$
$\displaystyle\nu_{\rho}\nabla^{2}h_{2}+\lambda_{1\rho}\partial_{\parallel}h_{2}-\frac{\lambda_{\rho}}{2}(\bm{\nabla}h_{2})^{2}+\lambda_{2\rho}[\delta
m(\partial_{\parallel}h_{2})]$ (10) $\displaystyle-$
$\displaystyle\lambda_{\rho\psi}(\bm{\nabla}h_{1})\cdot(\bm{\nabla}h_{2})-\lambda_{3\rho}[\delta
m(\bm{\nabla}h_{2})^{2}]$ $\displaystyle-\lambda_{3\rho\psi}[\delta
m\bm{\nabla}h_{1}\cdot\bm{\nabla}h_{2}]+f.$

Here, we use different sets of model parameters associated $h_{1}$ and $h_{2}$
to allow for different scalings of $h_{1}$ and $h_{2}$ in long wavelength
limit. The parameters $\nu_{\psi},\nu_{\rho}$ are come from $\nu$;
$\lambda_{\psi}$, $\lambda_{\rho}$, $\lambda_{\rho\psi}$ are from $\lambda$;
$\lambda_{1\psi}$, $\lambda_{1\rho}$ are from $\lambda_{1}$;
$\lambda_{2\psi}$, $\lambda_{2\rho}$ are from $\lambda_{2}$ and
$\lambda_{3\psi}$, $\lambda_{3\rho}$, $\lambda_{3\rho\psi}$ are from
$\lambda_{3}$.

Equations (9) and (10) obviously generalise the corresponding $1d$ equations
in Ref. Haldar and Basu (2020). Field $h_{1}({\bf x})$ represents the frozen
height of a surface which is entirely “driven” by the quenched disorder
$\delta m({\bf x})$, whereas $h_{2}$ represents the time-dependent height of a
growing surface, driven by the time-dependent additive noise $f$. Clearly, the
dynamics of the height field $h_{2}({\bf x},t)$ in Eq. (10) is affected by the
frozen height field $h_{1}({\bf x},t)$ that satisfies Eq. (9), and hence by
the quenched disorder.

We define two cases depending on the choices of parameters. Case-I: Here,
$\lambda_{1}\neq 0$, $\lambda_{2}\neq 0$ hence
$\lambda_{1\psi},\,\lambda_{1\rho},\,\lambda_{2\psi},\,\lambda_{2\rho}$ are
nonzero, leading to the presence of propagating modes. Case-II: here,
$\lambda_{1}=0$, $\lambda_{2}=0$ mean
$\lambda_{1\psi}=\lambda_{1\rho}=\lambda_{2\psi}=\lambda_{2\rho}=0$,
propagating modes vanish and the full $d$-dimensional isotropic property of
surface is restored. These two cases are corresponding to the away from half
filled limit and near to the half-filled limit of the periodic TASEP Haldar
and Basu (2020) respectively.

##  III Kardar-Parisi-Zhang universality class

We first briefly review the KPZ universality before analysing our disordered
model. The KPZ equation is given by Kardar et al. (1986)

$\frac{\partial h}{\partial
t}=\nu\nabla^{2}h-\frac{\lambda}{2}({\bm{\nabla}}h)^{2}+f.$ (11)

Here, $h({\bf x},t)$ is height of $d$-dimensional surface at any instant time
$t$ and driven by Gaussian white noise $f$ with satisfying (8b). This model
equation admits Galilean invariance due to transformation of the inertial
frame. Correlation function of the height fluctuations follows the scaling in
the long wavelength limit:

$\langle[h({\bf x},t)-h(0,0)]^{2}\rangle=|{\bf x}|^{2\chi}\varTheta(|{\bf
x}|^{z}/t).$ (12)

Here, $\chi$ and $z$ are roughness and dynamic exponents for $h$; $\varTheta$
is a dimensionless scaling function of its argument. One could find the
scaling behaviour of $h({\bf x},t)$ in the long time limit by defining a
dimensionless coupling constant $g=\frac{\lambda^{2}D}{\nu^{3}}$. The RG flow
equation for $g$ in the one-loop perturbative theory satisfies

$\frac{dg}{dl}=g\left[2-d+g\frac{4d-6}{d}\right].$ (13)

Notice from (13) that $g$ diverges at $d=3/2$; furthermore, between $3/2<d<2$
$g=0$, is the only physically acceptable solution. This is believed to be an
artefact of the one-loop perturbation theory Frey and Täuber (1994). The
Galilean invariance of the KPZ equation leads to $\chi+z=2$, an exact relation
between the scaling exponents. At $1d$, $z=3/2$ and $\chi=1/2$ are exactly
known analytically due to the fluctuation-dissipation-theorem in the model.
Coupling $g$ grows under rescaling and is marginally relevant at $d=2$, the
lower critical dimension of the KPZ equation. The KPZ equation admits a
roughening transition (smooth-to-rough phase transition) for dimensions higher
than $2d$; the rough phase is inaccessible in the perturbative theory.

##  IV Universal properties and scaling behaviour

The autocorrelation functions of $h_{1}({\bf x})$ and $h_{2}({\bf x},t)$ are
characterised by the universal scaling exponents. The autocorrelation
functions are given by

$\displaystyle C_{1}({\bf x})\equiv\overline{[h_{1}({\bf
x})-h_{1}(0)]^{2}}=|{x}_{\parallel}|^{2\chi_{1}}f_{1}\left(\frac{|{\bf
x}|^{\mu_{1}}_{\perp}}{x_{\parallel}}\right),$ (14) $\displaystyle C_{2}({\bf
x},t)\equiv\overline{\langle[h_{2}({\bf x},t)-h_{2}(0,0)]^{2}\rangle}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=|x_{\parallel}|^{2\chi_{2}}f_{2}(|{x_{\parallel}}|^{z}/t,\frac{|{\bf
x}|^{\mu_{2}}_{\perp}}{x_{\parallel}}).$ (15)

Here, $f_{1}$ and $f_{2}$ are two dimensionless scaling functions of their
arguments; ${\bf x}=(x_{\parallel},{\bf x}_{\perp})$ with $x_{\parallel}$
being the direction of propagation of the traveling waves and ${\bf
x}_{\perp}$ being $d-1$ transverse directions to the propagating direction of
waves. The scaling functions represent the generic anisotropy in the presence
of the propagating waves. Scaling exponents $\chi_{1}$ and $\chi_{2}$ are
roughness exponents for $h_{1}$ and $h_{2}$ respectively and $z$ is dynamic
exponent of $h_{2}$; exponents $\mu_{1},\,\mu_{2}$ are the anisotropy
exponents with $\mu_{2}=1$ in the pure $d$-dimensional KPZ equation (since it
is fully isotropic).

###  IV.1 Linear theory

The equal-time autocorrelation functions of $h_{1}({\bf x})$ and $h_{2}({\bf
x},t)$ can be written down exactly in the linear theory using the linear terms
in (9) and (10). These in the Fourier space for Case-I are

$\displaystyle\langle|h_{1}({\bf
k})|^{2}\rangle=\frac{2\tilde{D}\lambda^{2}_{m}}{\lambda_{1\psi}^{2}k_{\parallel}^{2}+\nu_{\psi}^{2}k_{\perp}^{4}},$
(16a) $\displaystyle\langle|h_{2}({\bf
k},t)|^{2}\rangle=\frac{D}{\nu_{\rho}k^{2}},$ (16b)

where in (16a) we have ignored $k_{\parallel}^{4}$ in comparison with
$k_{\parallel}^{2}$ in the long wavelength limit. Naturally, the exponents
$\chi_{1}=(2-d)/2$ and $\chi_{2}=(2-d)/2$ in the linear theory are known
exactly. Furthermore, we identify that $\mu_{1}=2$ and $\mu_{2}=1$ in the
linear theory.

Similarly, the equal-time autocorrelation functions for Case-II are also known
exactly in the linear theory. These are given by

$\displaystyle\langle|h_{1}({\bf
k})|^{2}\rangle=\frac{2\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}k^{4}}.$ (17a)
$\displaystyle\langle|h_{2}({\bf k},t)|^{2}\rangle=\frac{D}{\nu_{\rho}k^{2}}$
(17b)

in the Fourier space. We thus find that $\chi_{1}=(4-d)/2$ and
$\chi_{2}=(2-d)/2$ in linear theory; in this case $\mu_{1}=1=\mu_{2}$. $z$
remains at 2 and $\chi_{2}$ has the same value in the linear theory in both
the cases.

###  IV.2 Nonlinear effects

The nonlinear terms in (9) and (10) may change the scaling exponents found in
the linear theory. The exponents can no longer be found exactly in the
presence of the nonlinear terms. Furthermore, naïve perturbation theory
produces diverging corrections to the model parameters. To handle these
divergences systematically, we here use perturbative dynamical RG approach to
find the effects of anharmonic terms in scaling exponents for long wavelength
limit.

The RG method is well documented in the literature Chaikin and Lubensky
(2000); Barabási and Stanley (1995); Forster et al. (1977); we briefly outline
the steps of the RG analysis. The momentum shell dynamic RG consists of few
steps, at first the fields are expressed as the sum of low and high wavevector
modes: $h_{1}({\bf q})=h_{1}^{<}({\bf q})+h_{1}^{>}({\bf q})$ and $h_{2}({\bf
q},\omega)=h_{2}^{<}({\bf q},\omega)+h_{2}^{>}({\bf q},\omega)$. Here $<$, $>$
present low and high wavevector modes respectively. Let the upper limit value
of momentum is $\Lambda$ related to $1/a$, the microscopic cutoff in the
system. We integrate out the high wavevector Fourier modes of fields
restricted $\Lambda e^{-\delta l}<|{\bf q}|<\Lambda$, these are represented
diagrammatically by the Feynman diagrams which are given in Appendix A. This
integration reduces the upper cut off of wavevector to $\Lambda e^{-\delta l}$
from $\Lambda$, and then we restore the upper limit of wavevector to $\Lambda$
by rescaling of momentum by ${\bf q^{\prime}}={\bf q}e^{\delta l}$, along with
frequency by $\omega^{\prime}=\omega e^{z\delta l}$, where $z$ is the dynamic
exponent. The long wavelength parts of the fields are also rescaled; the
necessary details are available in Appendix A. This procedure ultimately
yields the RG flow equations, which in turn gives the RG fixed points.

###  IV.3 Universality class of Case-I

In this case, the surface is anisotropic due to the existence of the
travelling waves with nonzero $\lambda_{1\psi}$, $\lambda_{1\rho}$. We write
the model equations keeping the most relevant terms in the long wavelength
limit below.

$\displaystyle-\nu_{\psi}\nabla^{2}h_{1}-\lambda_{1\psi}\partial_{\parallel}h_{1}-\lambda_{m}\delta
m+\frac{\lambda_{\psi}}{2}(\bm{\nabla}h_{1})^{2}$
$\displaystyle-\lambda_{2\psi}[\delta m(\partial_{\parallel}h_{1})]=0,$ (18)

and

$\displaystyle\frac{\partial h_{2}}{\partial t}$ $\displaystyle=$
$\displaystyle\nu_{\rho}\nabla^{2}h_{2}+\lambda_{1\rho}\partial_{\parallel}h_{2}-\frac{\lambda_{\rho}}{2}(\bm{\nabla}h_{2})^{2}+\lambda_{2\rho}[\delta
m(\partial_{\parallel}h_{2})]$ (19) $\displaystyle-$
$\displaystyle\lambda_{\rho\psi}(\bm{\nabla}h_{1})\cdot(\bm{\nabla}h_{2})+f.$

The perturbative corrections of these parameters which involve disorder
(proportional to $\tilde{D}$) are finite; infrared divergent corrections are
due to the pure KPZ nonlinear term ($\frac{\lambda_{\rho}}{2}(\nabla
h_{2})^{2}$) only, see Appendix B. Now to determine the scaling of the
dimensionless effective coupling constants near the Gaussian fixed point, that
depend upon the disorder variance $\tilde{D}$. These are
$a_{0}=\frac{\lambda_{\psi}^{2}\tilde{D}\lambda_{m}^{2}}{\lambda_{1\psi}^{4}}$,
$a_{1}=\frac{\tilde{D}\lambda_{2\psi}^{2}}{\lambda_{1\psi}^{2}}$,
$a_{2}=\frac{\lambda_{\rho\psi}^{2}\tilde{D}\lambda_{m}^{2}}{\lambda_{1\rho}^{2}\lambda_{1\psi}^{2}}$
, $a_{3}=\frac{\lambda_{2\rho}^{2}\tilde{D}}{\lambda_{1\rho}^{2}}$,
$a_{4}=\frac{\lambda_{\psi}\lambda_{2\psi}\tilde{D}\lambda_{m}}{\lambda^{3}_{1\psi}}$.
Under the rescaling of the momentum, frequency and the fields, the flow of
these near the Gaussian fixed point Haldar and Basu (2020)

$\displaystyle\frac{da_{0}}{dl}=-da_{0},\,\frac{da_{1}}{dl}=-da_{1},$
$\displaystyle\frac{da_{2}}{dl}=-da_{2},\,\frac{da_{3}}{dl}=-da_{3},\,\frac{da_{4}}{dl}=-da_{4}.$
(20)

Thus all these coupling constants are irrelevant near the Gaussian fixed point
at all dimensions. The Galilean invariance of the dynamics of $h_{2}({\bf
x},t)$ is restored in the long wavelength limit due to the irrelevance of the
quenched disorder. The flow of the dimensionless coupling constant
$\tilde{g}=\frac{\lambda_{\rho}^{2}D}{\nu_{\rho}^{3}}$, which appears in the
pure KPZ problem as well, near the Gaussian fixed point is

$\frac{d\tilde{g}}{dl}=\tilde{g}[2-d].$ (21)

This is unstable for $d<2$, but stable for $d>2$. At any rate, even for $d>2$,
$\tilde{g}$ remains more relevant near the Gaussian fixed point than the other
dimensionless coupling constants as defined above, all of which have their
origins in the disorder. We can thus conclude that disorder is irrelevant (in
a RG sense) in the long wavelength limit at all dimensions for Case-I. This
immediately implies that the long wavelength scaling is governed by the KPZ
nonlinearity with the universal scaling belonging to the KPZ universality
class. Thus, a roughening transition is expected at $d>2$, just as it is for
the pure KPZ equation, to a rough phase identical to that for the pure KPZ
equation, whose scaling properties are not accessible to the standard
perturbative calculations. Furthermore and related to what we have just
concluded, isotropy is restored in the long wavelength limit of the
fluctuations of $h_{2}$. This implies $\mu_{2}=1$ in the renormalised
equation, as in the pure KPZ equation. This generalises one of the conclusions
of Ref. Haldar and Basu (2020) valid for $1d$. We therefore conclude that
while the Gaussian fixed point is stable with respect to perturbations by the
quenched disorder, the pure KPZ nonlinearity (with coupling $\tilde{g}$)
remains relevant and determines the universal scaling properties. Thus the
long wavelength properties of the model is given by the KPZ universality
class. Therefore lower critical dimension is 2 and a roughening transition
(i.e., a smooth-to-rough transition) exists as in the pure KPZ equation. The
upper critical dimension of Case-I, identical to that for the pure KPZ
equation, is not known. It is speculated to be 4 by some authors Lässig and
Kinzelbach (1997); Bhattacharjee (1998); Moser et al. (1991); Marinari et al.
(2002); Wiese (1998), although it remains debatable till the date.

###  IV.4 Universality class of Case-II

We now discuss the universal scaling properties of Case-II. In this case, the
underdamped propagating modes vanish (since $\lambda_{1}\rho=0$) in the
dynamics of $h_{2}({\bf x},t)$ as discussed before in Sec. II. Thus the
surface is isotropic. We write down below the model equations keeping most
relevant nonlinear terms:

$\displaystyle-\nu_{\psi}\nabla^{2}h_{1}-\lambda_{m}\delta
m+\frac{\lambda_{\psi}}{2}(\bm{\nabla}h_{1})^{2}=0.$ (22)

$\displaystyle\frac{\partial h_{2}}{\partial
t}=\nu_{\rho}\nabla^{2}h_{2}-\lambda_{\rho\psi}[\bm{\nabla}h_{1}\cdot\bm{\nabla}h_{2}]+f.$
(23)

The naïve perturbative fluctuation corrections to the model parameters diverge
in the long wavelength limit, which can be systematically handled within the
RG framework. Here we perform one loop RG analysis; fluctuation corrections of
parameters are represented by the one loop Feynman diagrams whose details are
available in Appendix C. The flow equations of the model parameters are

$\displaystyle\frac{d\nu_{\psi}}{dl}=\nu_{\psi}\left[d-2+(\frac{2}{d}-1)g_{1}\right].$
(24a)
$\displaystyle\frac{d\lambda_{m}}{dl}=\lambda_{m}\left[-\chi_{1}-\frac{3d}{2}+\frac{g_{1}}{2}\right].$
(24b)
$\displaystyle\frac{d\lambda_{\psi}}{dl}=\lambda_{\psi}\left[-d-2+\chi_{1}-\frac{2g_{1}}{d}\right].$
(24c)
$\displaystyle\frac{d\nu_{\rho}}{dl}=\nu_{\rho}\left[z-2+(\frac{2}{d}-1)g_{2}\right].$
(24d) $\displaystyle\frac{dD}{dl}=D\left[-d+z-2\chi_{2}+2g_{2}\right].$ (24e)
$\displaystyle\frac{d\lambda_{\rho\psi}}{dl}=\lambda_{\rho\psi}\left[z-2+\chi_{1}-\frac{2g_{2}}{d}\right].$
(24f)

Here, the dimensionless coupling constants are
$g_{1}=\frac{\lambda^{2}_{\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{4}_{\psi}}\frac{\Omega_{d}}{(2\pi)^{4}}$,
$g_{2}=\frac{\lambda^{2}_{\rho\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}\nu^{2}_{\rho}}\frac{\Omega_{d}}{(2\pi)^{4}}$
involving $\tilde{D}$, where $\Omega_{d}$ is angular part of the
$d$-dimensional loop integrals. The flow equations of $g_{1}$ and $g_{2}$ are

$\displaystyle\frac{dg_{1}}{dl}=g_{1}[4-d+(5-12/d)g_{1}],$ (25)
$\displaystyle\frac{dg_{2}}{dl}=g_{2}[4-d+(2-8/d)g_{2}+(3-4/d)g_{1}].$ (26)

Flow equations (25) and (26) clearly reveal that $d=4$ is the critical
dimension of both $g_{1}$ and $g_{2}$. Let us now calculate the fixed point
(FP) values of $g_{1}$ and $g_{2}$ at different dimensions by setting
$dg_{1}/dl=0=dg_{2}/dl$. The Gaussian FP is $g_{1}^{*}=0=g_{2}^{*}$,
furthermore one can find easily nontrivial FP as
$g_{1}^{*}=\frac{4-d}{12/d-5}$ from (25) and using this in (26)
$g_{2}^{*}=\frac{4-d}{12/d-5}$ same as $g_{1}^{*}$, shown in Fig. 1. Notice
that the nontrivial fixed point diverges at $d=12/5$ and physically
unacceptable between the dimensions $12/5$ and $4$.

Figure 1: Nontrivial fixed point value of $g_{1}$ as a function of $d$
(surface dimension) in one loop calculation. This is also same for $g_{2}$;
see text.

We find in detail below for different dimensions:

i) $d\leqslant 12/5$: In this regime the Gaussian FP $g_{1}^{*}=0=g_{2}^{*}$
is unstable. Thus $g_{1}^{*}$ and $g_{2}^{*}$ should be non-zero in the
statistical steady state, implying relevance (in a RG sense) of the quenched
disorder. The nontrivial FP is linearly stable. The scaling exponents at the
stable fixed point are given below.

(a) $d=1$: The exponents are $\chi_{1}=9/7$, $\chi_{2}=5/7$, $z=11/7$. Thus
$\chi_{2}+z=16/7>2$, in contrast to the pure $1d$ KPZ equation. The results
are consistent with Haldar and Basu (2020).

(b) $d=2$: The exponents are $\chi_{1}=2$, $\chi_{2}=2$, $z=2$. Also here
$\chi_{2}+z>2$, unlike the pure KPZ equation.

(ii) Dimension $12/5\leqslant d<4$: For this intermediate range of dimensions,
physically acceptable fixed point is only $g_{1}^{*}=0$, $g_{2}^{*}=0$, which
however is an unstable fixed point. The perturbation theory apparently breaks
down for $12/5\leqslant d<4$, which we believe is an artifact of the one-loop
expansion employed here. This is similar to what we find for $3/2<d<2$ in the
one-loop RG calculations for the pure KPZ equation; see Eq. (13) Frey and
Täuber (1994). A schematic RG flow diagram is shown in Fig. 2 below.

Figure 2: Schematic RG flow diagram for $d<4$. Stable and unstable fixed
points are shown; flow lines indicate the direction of the RG flows towards
stable FP (see text).

(iii) Dimension $d=4$: It is the lower critical dimension of the presented
model here. The nonlinear couplings are marginally relevant, are shown in (25)
and (26). This implies that at $d=4$, the coupling constants flow to infinity
in a finite RG “time” or for a finite system size, that is controlled by the
precise values of the unrenormalised model parameters (and hence non-
universal). This is reminiscent of the flow of the coupling constant at $2d$
in the pure KPZ equation.

(iv) Dimension $d>4$: Let $d=4+\epsilon$, where $\epsilon>0$ is assumed to be
small. The Gaussian fixed points are now stable, corresponding to $z=2$,
$\chi_{1}=2-d/2$, and $\chi_{2}=1-d/2$, corresponding to asymptotically smooth
phases. In addition, similar to the pure KPZ equation, we expect unstable
fixed points with both $g_{1}^{*},\,g_{2}^{*}\sim{\cal O}(\epsilon)$, that
would indicate a disorder induced roughening transition in the model. In fact,
we do find $g_{1}^{*}=\epsilon/2$. This however leaves $g_{2}$ marginal.
However, given the fact that $g_{2}^{*}=0$ is a stable fixed point (along with
$g_{1}^{*}=0)$, and $g_{1}^{*}$ has an unstable fixed point that is ${\cal
O}(\epsilon)$, we speculate that $g_{2}$ too shows a roughening transition
such that for large enough bare or unrenormalised value of $g_{2}$, it flows
to infinity under renormalisation. Our one-loop RG appears to be inadequate to
capture this behaviour; higher order perturbation theory or numerical
solutions of the dynamical equation could be useful in this regard. This
further indicates that $d=4$ is the lower critical dimension for both $g_{1}$
and $g_{2}$.

Some studies Lässig and Kinzelbach (1997); Bhattacharjee (1998); Moser et al.
(1991); Marinari et al. (2002); Wiese (1998) suggest $d=4$ as the upper
critical dimension on or above which the scaling exponents from the linear
theory holds, for the pure KPZ equation. But here $d=4$ is found to be the
lower critical dimension, which is higher than the lower critical dimension of
the pure KPZ equation. Since the upper critical dimension of a model should be
higher than the lower critical dimension, we expect the upper critical
dimension for Case II should be greater than $d=4$.

In the pure KPZ equation, the stable fixed point that governs the rough phase
is inaccessible in a perturbation theory, although such a fixed point should
exist on physical ground. Similarly, here too a globally stable fixed point
that characterises the rough phase must exist. Using an “Occam’s razor” style
argument, we speculate the RG flow diagram Fig. 3 that has the physically
acceptable and simplest flow topology. Notice that at this putative strong
coupling fixed point, both $g_{1}$ and $g_{2}$ are expected to be non-zero,
for if $g_{2}=0$ at this fixed point, the rough phase should become identical
to that of the pure KPZ equation. Such a possibility is however ruled out,
since a strong coupling perturbatively inaccessible rough phase appears
already at $2d$ for the pure KPZ equation, whereas for Case-II here, it can
appear only at $4d$ or above. Thus the rough phase here is expected to be
different from its counterpart for the pure KPZ equation.

Figure 3: Speculated global RG flows in the $g_{1}-g_{2}$ plane in an Occam’s
razor style argument. Perturbative accessible (small filled black circles) and
inaccessible (small coloured squares) are shown. RG flow lines with directions
are marked. Notice not all the unstable and in-principle perturbatively
accessible fixed points are captured by the one-loop theory; see text.

Finally, we give a schematic, pictorial comparison between the phases and
phase transitions in Case-I and Case-II as a function of dimension $d$ in Fig.
4 below.

Figure 4: Phases and phase transitions in the model as a function of $d$.
(top) Case-I and (bottom) Case-II. Lower ($d_{lc_{I}}=2$) and
($d_{lc_{II}}=4$) dimensions for Case-I and Case-II, respectively, and (yet
unknown) upper critical ($d_{uc_{I}}$) and $d_{uc_{II}}$ dimensions,
respectively, for Case-I and Case-II are marked. Further, $d_{uc_{II}}$ is
expected to be higher than $d_{uc_{I}}$; see text.

We expect the new universality class that emerges in the absence of the
propagating waves holds not only when these waves strictly vanish, but also
when these have “small” amplitudes, i.e., a nonuniversal “small” window of
speeds of propagation around zero, for reasons identical to those enunciated
in Ref. Haldar and Basu (2020). Essentially for sufficiently high wavevectors,
the propagating waves are subdominant to the damping terms. If under mode
elimination from the high wavevector limits, the theory gets renormalised
before, the crossover scale $\xi_{c}$ between the propagating modes and
damping terms, the “renormalised” theory at the crossover scale would like
what we have obtained for Case-II. Since at length scales larger than
$\xi_{c}$, the propagating modes are important, and no further renormalisation
takes place. As a result, the scaling behaviour corresponding to our Case-II
above will be displayed by the system, with a non-KPZ scaling behaviour. If,
on the other hand, the theory does not get “renormalised” till the crossover
scale $\xi_{c}$ is reached, further renormalisation of the theory at scales
larger than $\xi_{c}$ takes place, which comes entirely from the KPZ
nonlinearity. This naturally leads to a scaling behaviour identical to the
pure KPZ equation; see Ref. Haldar and Basu (2020) for more detailed
discussions.

##  V Summary and outlook

To summarise, we have studied here the universal scaling properties of
$d$-dimensional KPZ equation with quenched columnar disorder. For this, we
have generalised the $1d$ continuum equation constructed and studied in Ref.
Haldar and Basu (2020). We have obtained a number of interesting results. For
instance, we show that the columnar disorder in general leads to the loss of
the Galilean invariance of the pure KPZ equation and generation of underdamped
propagating waves in the system, which in turn makes the system anisotropic.
Interestingly, these waves render the quenched disorder irrelevant (in a RG
sense); as a result, the universal scaling exponents belong to the
$d$-dimensional KPZ universality class. Since the pure KPZ equation is
isotropic, this implies that the long wavelength scaling properties of the
model are actually isotropic. We argue that the rough phase of the model that
is inaccessible in a perturbation theory, is statistically identical to that
for the pure KPZ equation in the long wavelength limit. Thus isotropy becomes
an effective symmetry in the long wavelength limit. For certain choices of the
model parameters, the propagating waves vanish. In that limit the model is
already isotropic. Furthermore, the quenched disorder in the absence of the
waves is now relevant (in a RG sense). We show that the model now has $d=4$ as
the lower critical dimension. We have calculated the roughness and the dynamic
scaling exponents within a one-loop approximation, which belong to a
universality class hitherto unknown. We argue that above $d=4$ the model
undergoes a roughening transition from a smooth to rough phase that is
analogous to the well-known roughening transition in the KPZ equation above
$d=2$. This rough phase, although not accessible in a perturbation theory, is
statistically different from its counterpart in the pure KPZ equation. We have
argued that the upper critical dimension in Case II should be higher than 4.

We have considered only Gaussian-distributed short-ranged quenched disorder.
Our calculational framework can be extended to Gaussian-distributed long-range
disorder in a straight forward manner. While precise values of the scaling
exponents should depend upon the scaling of the variance of the long-range
disorder, by using the logic outlined above we can generally argue that the
universal scaling properties will crucially depend upon whether or not there
are underdamped propagating waves. We expect our results here will provide
impetus to future theoretical work along this direction.

##  VI acknowledgements

The author would like to acknowledge Prof. Abhik Basu (SINP, Kolkata) for
helpful discussions, valuable suggestions and critical comments to prepare the
manuscript.

##  Appendix A Perturbation Theory: Generating Functional and Scaling

The Generating Functional Zinn-Justin (2002); Täuber (2014); DeDominicis and
Martin (1979); Bausch et al. (1976) is defined as

${\mathcal{Z}}=\int{\cal D}\hat{h}_{1}{\cal D}h_{1}{\cal D}\hat{h}_{2}{\cal
D}h_{2}{\cal D}\delta m\exp\\{-{\cal
S}[\hat{h}_{1},h_{1},\hat{h}_{2},h_{2},\delta m]\\}.$ (27)

$\hat{h}_{1},\,\hat{h}_{2}$ are the conjugate fields to $h_{1},\,h_{2}$
respectively. Here, ${\cal S}[\hat{h}_{1},h_{1},\hat{h}_{2},h_{2},\delta m]$
is the action functional. Two types of terms are present in the action ${\cal
S}$: linear terms and nonlinear terms. The perturbation theory is set up by
expanding the nonlinear terms present in the action.

The perturbative step of evaluating fluctuation corrections is followed by
rescaling momentum (equivalently, space) and frequency (equivalently, time) as
${\bf q}\rightarrow b{\bf q}$ and $\omega\rightarrow b^{z}\omega$
respectively. The long wavelength parts of the fields are rescaled as follows:

$\displaystyle\hat{h}_{1}({\bf q})=b^{-\chi_{1}-d}\hat{h}_{1}(b{\bf
q}),\;h_{1}({\bf q})=b^{d+\chi_{1}}h_{1}(b{\bf q}),$ $\displaystyle\delta
m({\bf q})=b^{d/2}\delta m(b{\bf q}),$ $\displaystyle\hat{h}_{2}({\bf
q},\omega)=b^{z-\chi_{2}}\hat{h}_{2}(b{\bf q},b^{z}\omega),$ $\displaystyle
h_{2}({\bf q},\omega)=b^{d+z+\chi_{2}}h_{2}(b{\bf q},b^{z}\omega).$ (28)

##  Appendix B Renormalisation Group analysis for Case-I

The action functional corresponding to equations (18) and (19) is

$\displaystyle{\cal S}[\hat{h}_{1},h_{1},\hat{h}_{2},h_{2},\delta m]=\int_{\bf
x}\hat{h}_{1}[-\nu_{\psi}\nabla^{2}h_{1}-\lambda_{1\psi}\partial_{\parallel}h_{1}$
$\displaystyle-\lambda_{m}\delta
m+\frac{\lambda_{\psi}}{2}(\bm{\nabla}h_{1})^{2}-\lambda_{2\psi}[\delta
m(\partial_{\parallel}h_{1})]+\frac{\delta m^{2}}{4\tilde{D}}$
$\displaystyle+\int_{{\bf
x},t}\hat{h}_{1}[-D\hat{h}_{2}+\partial_{t}h_{2}-\nu_{\rho}\nabla^{2}h_{2}-\lambda_{1\rho}\partial_{\parallel}h_{2}+\frac{\lambda_{\rho}}{2}(\bm{\nabla}h_{2})^{2}$
$\displaystyle-\lambda_{2\rho}[\delta
m(\partial_{\parallel}h_{2})]+\lambda_{\rho\psi}(\bm{\nabla}h_{1})\cdot(\bm{\nabla}h_{2})].$
(29)

The two point functions found from linear terms present in (29) and nonlinear
vertices in (29) are presented diagramatically in Fig. 5 and Fig. 6
respectively.

Figure 5: Diagramatic representations of two point functions. Arrows in
propagators ((a), (e)) constitute casuality information.

Figure 6: Diagramatic representations of the anharmonic vertices terms in
action 29.

The bare two point functions in the harmonic theory neglecting nonlinear terms
from (29) are

$\displaystyle\langle\hat{h}_{1}(-{\bf k})h_{1}({\bf
k})\rangle=\frac{1}{\nu_{\psi}k^{2}+ik_{\parallel}\lambda_{1\psi}}.$ (30a)
$\displaystyle\langle|h_{1}({\bf
k})|^{2}\rangle=\frac{2\tilde{D}\lambda_{m}^{2}}{k_{\parallel}^{2}\lambda_{1\psi}^{2}+\nu_{\psi}^{2}k^{4}}.$
(30b) $\displaystyle\langle h_{1}(-{\bf k})\delta m({\bf
k})\rangle=\frac{2\tilde{D}\lambda_{m}}{\nu_{\psi}k^{2}-ik_{\parallel}\lambda_{1\psi}}.$
(30c) $\displaystyle\langle\hat{h}_{2}(-{\bf k},-\omega)h_{2}({\bf
k},\omega)\rangle=\frac{1}{-i\omega+ik_{\parallel}\lambda_{1\rho}+\nu_{\rho}k^{2}}.$
(30d) $\displaystyle\langle|h_{2}({\bf
k},\omega)|^{2}\rangle=\frac{2D}{(\omega-
k_{\parallel}\lambda_{1\rho})^{2}+\nu_{\rho}^{2}k^{4}}.$ (30e)

The fluctuation corrections of the propagators and correlators of
$h_{1},\,h_{2}$ are represented by one-loop Feynman diagrams in Fig. 7 and
Fig. 8 respectively.

Figure 7: One loop Feynman diagrams for the corrections of the propagator and
correlator of $h_{1}$: (a),(b) for propagator and (c),(d),(e) for correlator
corrections.  Figure 8: One loop Feynman diagrams for the corrections of the
correlator and propagator of $h_{2}$: (a),(b) for correlator and (c),(d) for
propagator corrections.

The one loop contributions to the corrections of parameters in (29) are
finite. Those are of the forms

$\displaystyle{\rm Fig.}(\ref{h1-corrections}a)$
$\displaystyle=\frac{\lambda^{2}_{\psi}\tilde{D}\lambda^{2}_{m}}{i\lambda^{3}_{1\psi}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{3}_{\parallel}}[k^{2}q^{2}+({\bf k}\cdot{\bf
q})^{2}$ $\displaystyle-q^{2}{\bf k}\cdot{\bf
q}(2+3k_{\parallel}/q_{\parallel})].$ (31a) $\displaystyle{\rm
Fig.}(\ref{h1-corrections}b)$
$\displaystyle=\frac{2\lambda^{2}_{2\psi}\tilde{D}}{\lambda_{1\psi}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}.$ (31b) $\displaystyle{\rm Fig.}(\ref{h1-corrections}c)$
$\displaystyle=\frac{\lambda^{2}_{\psi}(\tilde{D}\lambda^{2}_{m})^{2}}{\lambda^{4}_{1\psi}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}.$ (31c) $\displaystyle{\rm Fig.}(\ref{h1-corrections}d)$
$\displaystyle=\frac{2\lambda^{2}_{2\psi}\tilde{D}^{2}\lambda^{2}_{m}}{\lambda^{2}_{1\psi}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}.$ (31d) $\displaystyle{\rm Fig.}(\ref{h1-corrections}e)$
$\displaystyle=\frac{2\lambda_{\psi}\lambda_{2\psi}\tilde{D}\lambda^{2}_{m}}{\lambda^{3}_{1\psi}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{q^{2}}{q^{2}_{\parallel}}.$ (31e)

$\displaystyle{\rm Fig.}(\ref{h2-corrections}a)$
$\displaystyle=\frac{2D\tilde{D}\lambda_{\rho\psi}^{2}\lambda_{m}^{2}}{\lambda^{2}_{1\rho}\lambda^{2}_{1\psi}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{q^{2}}{q^{2}_{\parallel}}.$ (32a) $\displaystyle{\rm
Fig.}(\ref{h2-corrections}b)$
$\displaystyle=-\frac{2D\lambda^{2}_{2\rho}}{\lambda^{2}_{1\rho}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{q^{2}_{\parallel}}{q^{2}}.$ (32b) $\displaystyle{\rm
Fig.}(\ref{h2-corrections}c)$
$\displaystyle=\frac{i\tilde{D}\lambda^{2}_{\rho\psi}\lambda^{2}_{m}}{\lambda_{1\rho}\lambda^{2}_{1\psi}}\int\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{qq^{2}-\parallel}[-2({\bf k}\cdot{\bf
q})^{2}-k^{2}q^{2}]$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\frac{q({\bf
k}\cdot{\bf q})}{q^{2}_{\parallel}}[2+k/q+k_{\parallel}/q_{\parallel}].$ (32c)

##  Appendix C Renormalisation Group analysis for Case-II

Action functional corresponding to equations (22) and (23) is

$\displaystyle{\cal S}=\int_{\bf
x}\hat{h}_{1}[-\nu_{\psi}\nabla^{2}h_{1}-\lambda_{m}\delta
m+\frac{\lambda_{\psi}}{2}(\bm{\nabla}h_{1})^{2}]+\frac{\delta
m^{2}}{4\tilde{D}}$ $\displaystyle+\int_{{\bf
x},t}\hat{h}_{2}[-D\hat{h_{2}}+(\partial_{t}h_{2}-\nu_{\rho}\nabla^{2}h_{2}+\lambda_{\rho\psi}\bm{\nabla}h_{1}\cdot\bm{\nabla}h_{2})].$
(33)

Two point functions of the harmonic theory neglecting nonlinear terms fom (33)
are

$\displaystyle\langle\hat{h}_{1}(-{\bf k})h_{1}({\bf
k})\rangle=\frac{1}{\nu_{\psi}k^{2}}.$ (34a) $\displaystyle\langle|h_{1}({\bf
k})|^{2}\rangle=\frac{2\tilde{D}\lambda_{m}^{2}}{\nu_{\psi}^{2}k^{4}}.$ (34b)
$\displaystyle\langle h_{1}(-{\bf k})\delta m({\bf
k})\rangle=\frac{2\tilde{D}\lambda_{m}}{\nu_{\psi}k^{2}}.$ (34c)
$\displaystyle\langle\hat{h}_{2}(-{\bf k},-\omega)h_{2}({\bf
k},\omega)\rangle=\frac{1}{-i\omega+\nu_{\rho}k^{2}}.$ (34d)
$\displaystyle\langle|h_{2}({\bf
k},\omega)|^{2}\rangle=\frac{2D}{\omega^{2}+\nu_{\rho}^{2}k^{4}}.$ (34e)

The two point functions for this case are also represented diagramatically in
Fig.5. For this case, the nonlinear vertices are presented by diagrams (6a),
(6b) and (6d) in Fig.6.

In this case, the one loop diagrams in Fig.(7a), (7c) contribute to the
corrections of the propagator and correlator of $h_{1}$ respectively. The
diagrams in Fig.(8a), (8c) contribute to the corrections of correlator and
propagator of $h_{2}$ respectively. Fig.9, Fig.10 contribute to vertices
corrections in (33).

The one loop contributions to corrections of parameters in (33) are

$\displaystyle{\rm
Fig.}(\ref{h1-corrections}a)=\frac{\lambda^{2}_{\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{3}_{\psi}}[\frac{2}{d}-1]\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}.$ (35a) $\displaystyle{\rm
Fig.}(\ref{h1-corrections}c)=\frac{\lambda^{2}_{\psi}\tilde{D}^{2}\lambda^{2}_{m}}{\nu^{4}_{\psi}}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}.$ (35b)

$\displaystyle{\rm
Fig.}(\ref{h2-corrections}a)=\frac{2D\lambda^{2}_{\rho\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}\nu^{2}_{\rho}}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}.$ (36a) $\displaystyle{\rm
Fig.}(\ref{h2-corrections}c)=\frac{\lambda^{2}_{\rho\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}\nu_{\rho}}[1-2/d]\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}.$ (36b)

Figure 9: One loop Feynman diagram for the vertex $\lambda_{\psi}$ correction
for case-II.

$\displaystyle{\rm Fig.}(\ref{h1-vertex-
correction}a)=\frac{2\lambda^{3}_{\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{4}_{\psi}d}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}.$ (37a) $\displaystyle{\rm
Fig.}(\ref{h1-vertex-
correction}b)=-\frac{\lambda^{3}_{\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{4}_{\psi}d}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}.$ (37b)

Figure 10: One loop Feynman diagram for the vertex $\lambda_{\rho\psi}$
correction for case-II.

${\rm Fig.}(\ref{h2-vertex-
correction})=\frac{2\lambda^{3}_{\rho\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}\nu^{2}_{\rho}d}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}.$ (38)

Finally one loop corrected parameters in (33) are

$\displaystyle\nu^{<}_{\psi}=\nu_{\psi}\left[1+(\frac{2}{d}-1)\frac{\lambda^{2}_{\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{4}_{\psi}}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}\right].$ (39a)
$\displaystyle\lambda_{m}^{<}=\lambda_{m}\left[1+\frac{\lambda^{2}_{\psi}\tilde{D}\lambda^{2}_{m}}{2\nu^{4}_{\psi}}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}\right].$ (39b)
$\displaystyle\lambda_{\psi}^{<}=\lambda_{\psi}\left[1-\frac{\lambda^{2}_{\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{4}_{\psi}}\frac{2}{d}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}\right].$ (39c)
$\displaystyle\nu^{<}_{\rho}=\nu_{\rho}\left[1+(\frac{2}{d}-1)\frac{\lambda^{2}_{\rho\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}\nu^{2}_{\rho}}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}\right].$ (39d) $\displaystyle
D^{<}=D\left[1+\frac{2\lambda^{2}_{\rho\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}\nu^{2}_{\rho}}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}\right].$ (39e)
$\displaystyle\lambda_{\rho\psi}^{<}=\lambda_{\rho\psi}\left[1-\frac{\lambda^{2}_{\rho\psi}\tilde{D}\lambda^{2}_{m}}{\nu^{2}_{\psi}\nu^{2}_{\rho}}\frac{2}{d}\int^{\varLambda}_{\varLambda/b}\frac{d^{d}{\bf
q}}{(2\pi)^{d}}\frac{1}{q^{4}}\right].$ (39f)

## References

  * Lubensky (1975) T. C. Lubensky, Physical Review B 11, 3573 (1975).
  * Mukherjee and Basu (2020) S. Mukherjee and A. Basu, Phys. Rev. Research 2, 033423 (2020).
  * Grinstein and Luther (1976) G. Grinstein and A. Luther, Physical Review B 13, 1329 (1976).
  * Jeong et al. (1996) H. Jeong, B. Kahng, and D. Kim, Physical review letters 77, 5094 (1996).
  * Csahók et al. (1993) Z. Csahók, K. Honda, E. Somfai, M. Vicsek, and T. Vicsek, Physica A: Statistical Mechanics and its Applications 200, 136 (1993).
  * Barabási et al. (1996) A.-L. Barabási, G. Grinstein, and M. Munoz, Physical review letters 76, 1481 (1996).
  * Amaral et al. (1994) L. A. N. Amaral, A.-L. Barabási, and H. E. Stanley, Physical review letters 73, 62 (1994).
  * Haldar and Basu (2020) A. Haldar and A. Basu, Phys. Rev. Research 2, 043073 (2020).
  * Tripathy and Barma (1997) G. Tripathy and M. Barma, Physical review letters 78, 3039 (1997).
  * Tripathy and Barma (1998) G. Tripathy and M. Barma, Physical Review E 58, 1911 (1998).
  * de Queiroz and Stinchcombe (2008) S. de Queiroz and R. Stinchcombe, Physical Review E 78, 031106 (2008).
  * Kardar et al. (1986) M. Kardar, G. Parisi, and Y.-C. Zhang, Physical Review Letters 56, 889 (1986).
  * Frey and Täuber (1994) E. Frey and U. C. Täuber, Phys. Rev. E 50, 1024 (1994).
  * Chaikin and Lubensky (2000) P. M. Chaikin and T. C. Lubensky, _Principles of condensed matter physics_ , vol. 1 (Cambridge university press Cambridge, 2000).
  * Barabási and Stanley (1995) A.-L. Barabási and H. E. Stanley, _Fractal concepts in surface growth_ (Cambridge university press, 1995).
  * Forster et al. (1977) D. Forster, D. R. Nelson, and M. J. Stephen, Physical Review A 16, 732 (1977).
  * Lässig and Kinzelbach (1997) M. Lässig and H. Kinzelbach, Physical review letters 78, 903 (1997).
  * Bhattacharjee (1998) J. Bhattacharjee, Journal of Physics A: Mathematical and General 31, L93 (1998).
  * Moser et al. (1991) K. Moser, J. Kertész, and D. E. Wolf, Physica A: Statistical Mechanics and its Applications 178, 215 (1991).
  * Marinari et al. (2002) E. Marinari, A. Pagnani, G. Parisi, and Z. Rácz, Physical review E 65, 026136 (2002).
  * Wiese (1998) K. J. Wiese, Journal of statistical physics 93, 143 (1998).
  * Zinn-Justin (2002) J. Zinn-Justin, _Quantum field theory and critical phenomena_ , vol. 113 (Clarendon Press, Oxford, 2002).
  * Täuber (2014) U. C. Täuber, _Critical dynamics: a field theory approach to equilibrium and non-equilibrium scaling behavior_ (Cambridge University Press, 2014).
  * DeDominicis and Martin (1979) C. DeDominicis and P. Martin, Physical Review A 19, 419 (1979).
  * Bausch et al. (1976) R. Bausch, H.-K. Janssen, and H. Wagner, Zeitschrift für Physik B Condensed Matter 24, 113 (1976).