\section{Variation of the interface width scaling amplitude}
\label{secamp}


Once the squared interface width scales as $w_2 \simeq B t^{2\beta}$ in 1D KPZ interfaces, the amplitude $B$ can be obtained from the long time limit (when corrections to this scaling relation shall become negligible) of the ratios $w_2/t^{2\beta}$, which are shown in the insertion of Fig. 4a in the main manuscript. The values of $B$ estimated in this way for the SS model are depicted in Fig. S\ref{figS1}. Similar behavior is found for the other models. Thereby, in general, $B$ increases with both $\gamma$ (for fixed $\omega$) and $\omega$ (with $\gamma$ fixed). 

The convergence of the cumulant ratios to the GUE values in Fig. 4b of the main manuscript does not let room for doubt that the HDs are given by GUE, when $\gamma>1/z$. Henceforth, with the KPZ ansatz:
\begin{equation}
h=v_{\infty}t+s_{\lambda}(\Gamma t)^{\beta} \chi + \dots, 
\label{eqansatz}
\end{equation}
bearing in mind, which means that $w_2 \simeq (\Gamma t)^{2\beta} \langle \chi^2 \rangle_c$, we may conclude that $\Gamma=\left( B/\langle \chi^2 \rangle_c \right)^{1/2\beta}$ is also an increasing function of $\gamma$ and $\omega$, assuming that $\langle \chi^2 \rangle_c = 0.8132$ is a constant. The corresponding values of $\Gamma$ are also displayed (in right vertical axis) in Fig. S\ref{figS1}. It is noteworthy that for substrates expanding linearly in time ($\gamma=1$) the parameter $\Gamma$ (let us call it $\Gamma^*$) has a negligible dependence on $\omega$ and is the same as for flat (fixed-size - $\omega=0$) systems \cite{Ismael14}. Thus, for general $\gamma$, one may write
\begin{equation}
\Gamma(\gamma,\omega)=\Gamma^*+(\gamma-1) g(\gamma,\omega), 
\end{equation}
where $g(\gamma,\omega)$ is a positive function, so that $\Gamma < \Gamma^*$ ($\Gamma > \Gamma^*$) if $\gamma<1$ ($\gamma>1$), see Fig. S\ref{figS1}.

\begin{figure}[!h]
\includegraphics[width=7.cm]{FigS1.eps}
\caption{Variation of interface width amplitude $B$ (left) and the corresponding $\Gamma$ (right vertical axis) with the exponent $\gamma$ (main plot) and with the parameter $\omega$ (insertion) for the SS model. The dashed horizontal line indicates the exact value of $\Gamma$ (which is $\Gamma=1/2$) for this model in flat (fixed-size) substrates.}
\label{figS1}
\end{figure}





\section{Rescaling of the spatial covariance}
\label{secampcov}

Dynamic scaling predicts that $C_S/(\Gamma t)^{2\beta} \simeq \Phi(x)$, with $x \equiv A r^{2\alpha}/2 (\Gamma t)^{2\beta}$, where $A$ is the scaling amplitude of the height-height correlation function - defined as $C_H(r,t) \equiv \left\langle [h(x+r,t) - h(x,t)]^2 \right\rangle$ - with the distance $r$ [$C_H(r) \simeq A r^{2\alpha}$]. The scaling function $\Phi(x)$ is expected to be universal, being given by the Airy$_2$ (Airy$_1$) covariance in circular (flat) KPZ interfaces evolving in the plane. According to the KPZ theory, for 1D systems, $\Gamma = |\lambda| A^{2}/2$ \cite{Krug1992}, where $\lambda$ is the coefficient of the nonlinear term in the KPZ equation, usually referred to as the velocity excess. Thence, the ($\gamma,\omega$)-dependence in $\Gamma$ discussed above could imply in a similar dependence in $A$ and/or in $\lambda$. Unfortunately, it is hard to estimate accurately the value of $A$ from the scaling of $C_H$, which usually suffer from strong finite-size and -time corrections. We have indeed tried to do this, but the results (not shown) obtained do not allow us to conclude whether $A$ is a constant or ($\gamma,\omega$)-dependent. Furthermore, it is not clear to us how to determine the coefficient $\lambda$ in our expanding systems, because the tilting method \cite{barabasi} would not work here. Namely, by tilting the substrates, as they expand due to column duplications, their slope would decrease in time.

However, it does not matter what is the rescaling used, if it is appropriately applied to all data, it shall allow us to compare the covariances for different systems and to determine if they agree with Airy$_1$, Airy$_2$ or whatever. Thereby, in our work we compare the covariances by plotting $C_S/w_2$ against $r/a$, where $a$ is chosen so that at $r/a=1$ all curves have $C_S/w_2=0.2$. So, since $C_S(r=0)=w_2$, our rescaling obliges all curves to start at $C_S/w_2=1$ and to pass at $C_S/w_2=0.2$ at $r/a=1$. We remark that, if instead of the arbitrary $0.2$ point we had used $0.1$ or any other positive value, our conclusions on the covariances' behavior would still be the same.




\section{Derivation of the correction in the KPZ ansatz due to column duplication}
\label{seccorrec}

As demonstrated in Ref. \cite{Ismael14}, the KPZ non-linearity allied to the column duplication process yields a correction in the KPZ ansatz (Eq. \ref{eqansatz}), which is logarithmic when the substrate expands linearly in time. Here, we will generalize this by determining the main corrections for any $\gamma$. Let us start noticing that in systems with periodic boundary conditions the random column duplications cannot create a global curvature in the interface, so that $\langle \nabla^2 h \rangle=0$, but they can affect the nonlinear term $\langle (\nabla h)^2 \rangle$.

\begin{figure}[!t]
\includegraphics[width=7.cm]{FigS2.eps}
\caption{(Color online) Temporal variation of the correction $Q \equiv \partial_t \langle h \rangle - v_{\infty} + \beta s_{\lambda}\Gamma^{\beta} t^{\beta-1} \langle \chi \rangle$ for several models and parameters. The lines have the indicated slopes.}
\label{figS2}
\end{figure}

On a $d_s$-dimensional substrate, the average of the local gradient at time $t$ is given by
\begin{equation}
G_t=\frac{1}{L^{d_s}}\sum_{i=1}^{L^{d_s}} (\nabla h_i)^2.
\label{eqGt}
\end{equation}
If the number of duplications occurring in a time unity is (on average) $l=\gamma \omega t^{\gamma-1}$, then, at time $t+1$, one has
\begin{equation}
G_{t+1}=\frac{1}{(L+l)^{d_s}}\left[ \sum_{i=1}^{L^{d_s}}(\nabla h_i)^2+\sum_{i=L^{d_s}+1}^{(L+l)^{d_s}}(\nabla h_i)^2 \right],
\label{eqGt+1}
\end{equation}
where the first (second) summation in the rhs runs over the non-duplicated (duplicated) columns in that time interval. Disregarding the effects of particle deposition and considering the statistical equivalence of sites, one finds
\begin{equation}
\sum_{i=L^{d_s}+1}^{(L+l)^{d_s}}(\nabla h_i)^2 \approx\frac{d_s-1}{d_s} \left[ (L+l)^{d_s}-L^{d_s}\right]G_t,
\label{eqGradm}
\end{equation}
where the ratio $(d_s-1)/d_s$ is due to the local smoothing caused by the duplications, since immediately after the duplication of a column (say $j$) in a given direction $x_k$ one has $\left( \partial h_j/\partial x_k\right) =0$. From Eqs. \ref{eqGt} and \ref{eqGradm}, one may re-write \ref{eqGt+1} as
\begin{equation}
G_{t+1}\approx\left( 1-\frac{1}{d_s}\left[ 1-\left(\frac{L}{L+l} \right)^{d_s} \right]\right)G_t.
\end{equation}
Now, considering long times, such that $L\approx \omega t^{\gamma}$ and so $\frac{L+l}{L} \approx 1+\frac{\gamma}{t}$, and disregarding terms $\mathcal{O}(t^{-2})$, one obtains
$G_{t+1}-G_t\approx-\frac{\gamma}{t}G_t \quad \textnormal{or}\quad \frac{d G}{d t} \approx-\frac{\gamma}{t}G$,
which leads to $G_t\sim t^{-\gamma}$. Finally, since $\frac{\partial h}{\partial t} \sim G_t \sim t^{-\gamma}$ in the KPZ equation, a correction term $\zeta t^{1-\gamma}$ is expected in the KPZ ansatz (Eq. \ref{eqansatz}), where $\zeta$ is in principle a stochastic variable. Note that for $\gamma=1$ it turns out to be the logarithmic found in \cite{Ismael14}.

Therefore, the average height of the interfaces are expected to evolve as
\begin{equation}
\langle h \rangle =v_{\infty}t+s_{\lambda}(\Gamma t)^{\beta} \langle \chi \rangle + \langle \eta \rangle + \langle \zeta \rangle t^{1-\gamma} + \ldots, 
\label{eqansatz}
\end{equation}
where we have included also the well-known additional constant correction $\eta$ \cite{Takeuchi2010,Sasamoto2010,Alves11,Ismael14}. From the derivative of this equation in time, we find that $Q \equiv \partial_t \langle h \rangle - v_{\infty} + \beta s_{\lambda}\Gamma^{\beta} t^{\beta-1} \langle \chi \rangle \simeq (1-\gamma) \langle \zeta \rangle t^{-\gamma}$. This is indeed confirmed in Fig. S\ref{figS2}, where one sees that $Q$ decays in time as $Q \sim t^{-\gamma}$.



\section{Spatial covariance in random depositions out of the plane}
\label{secRD}

In order to analyze the spatial correlations generated by column duplications, we have performed simulations of a random deposition process \cite{barabasi} on substrates expanding as $\left\langle L(t) \right\rangle = L_0 + \omega t^{\gamma}$, following the same method used for the other models. In the RD model the aggregation rule is always simply $h_i \rightarrow h_i + 1$, where $i$ is a randomly chosen site, so that no correlation is generated by the deposition process. Indeed, $C_S(r,t)=0$ for any $r>0$ for this model in the flat ($\omega=0$) case. For expanding interfaces, however, we find non-trivial correlations. Although the rescaled covariance curves have negligible dependence on $\omega$ and on time, they change with the exponent $\gamma$ (see Fig. S\ref{figS3}). The curve for $\gamma=1.2$ is compared with those for the other models in Fig. 4d of the main manuscript.

\begin{figure}[!h]
\includegraphics[width=7.cm]{FigS3.eps}
\caption{(Color online) Rescaled spatial covariance for the random deposition model for the exponents $\gamma$ indicated.}
\label{figS3}
\end{figure}


