\section*{Methods}

The samples we used were ultra-high mobility SiGe/Si/SiGe quantum wells. The peak electron mobility in these samples reaches 240~m$^2$/Vs. The approximately 15~nm wide silicon (001) quantum well was sandwiched between Si$_{0.8}$Ge$_{0.2}$ potential barriers. The samples were patterned in Hall-bar shapes with the distance between the potential probes of 150~$\mu$m and width of 50~$\mu$m using standard photo-lithography. The long side of the Hall bar corresponded to the direction of current parallel to the [110] or [$\bar1$10] crystallographic axis. The in-plane magnetic field was perpendicular to the current to exclude the influence of ridges on the quantum well surface on the resistance measured at high electron densities (for more details, see Refs.~\cite{melnikov2015ultra,melnikov2017unusual}). Measurements were carried out in an Oxford TLM-400 dilution refrigerator with a base temperature of $\approx 25$~mK. Data were taken by a standard four-terminal lock-in technique in a frequency range 0.5--11~Hz in the linear regime of response.

\subsection*{Data availability.} The data that support the findings of this study are available from the corresponding author upon reasonable request.


\section{Introduction}

In this paper we establish the  local Lipschitz regularity of weak solutions of a certain class of quasilinear, degenerate parabolic equations in the Heisenberg group $\mathbb H^n$, or more in general in contact subRiemannian manifolds. In particular we extend to the non-stationary setting the early work \cite{Manfredi-Mingione, Mingione-Zatorska-Zhong}, by introducing  a new, more elegant approach to the regularity problem.

In a cylinder $Q=\Omega\times (0,T)$, where $\Omega\subset\mathbb H^n$ is an open set and $T>0$, we consider the equation
\begin{equation}\label{maineq-zero}
\p_t u = \sum_{i=1}^{2n} X_i A_i(x, \nabla_0 u) \quad \quad  \text{ in  \quad }Q=\Omega\times (0,T),
\end{equation}
modeled on the  parabolic $p$-Laplacian
\begin{equation}\label{plap-para} \p_t u= \sum_{i=1}^{2n} X_i \bigg(|\nabla_0 u|^{p-2} X_i u\bigg),\end{equation}
where  $ 2\le p \le 4$ and $X_1,...,X_{2n}$ denote the horizontal left invariant frame in $\mathbb H^n$.
In a previous study \cite{CCG}, Garofalo and the first two listed authors have extended techniques originally introduced by the third listed author \cite{Zhong} to establish $C^\infty$ smoothness for weak solutions to \eqref{maineq} in the range $2\le p <\infty$ under some additional non-degeneracy hypothesis.
 In the present paper we show that in the restricted range of the nonlinearity $2\le p \le 4$, and without the additional non-degeneracy assumptions, one can obtain Lipschitz regularity of weak solutions.

We indicate with $x=(x_1,...,x_{2n},x_{2n+1})$ the variable point in $\mathbb H^n$. We will occasionally denote the variable $x_{2n+1}$ in the center of the group with the letter $z$. Consequently, we will indicate with $\p_i$ partial differentiation with respect to the variable $x_i$, $i=1,...,2n$, and use the notation $Z = \p_z$ for the partial derivative $\p_{x_{2n+1}}$. The notation $\nabla_0 u=\sum_{i=1}^{2n} X_i u X_i \cong (X_1 u,...,X_{2n} u)$ denotes the so-called {\it horizontal gradient} of the function $u$, where
$$X_i = \partial_i - \frac{x_{n+i}}2 \partial_{z},\ \ \ \ \ \
X_{n+i} = \partial_{n+i} + \frac{x_{i}}{2} \partial_{z},\ \ \ \  i=1,...,n.$$
As it is well-known, the $2n+1$ vector fields $X_1,...,X_{2n}, Z$ are connected by the following commutation relation: for every couple of index $i, j$, if $j=i+n$, then $[X_i,X_j] = Z$; all other commutators being trivial.

\noindent{\bf Structural assumptions.}  The relevant  assumptions on the vector-valued function $$(x,\xi)\to A(x,\xi)=(A_1(x,\xi),...,A_{2n}(x,\xi))$$ are that
 there exist $2\le p \le 4$,  and  $0<\lambda'\le \Lambda' <\infty$ such that for a.e. $x\in\Omega, \xi\in \mathbb R^{2n}$ and for all $\eta\in \mathbb R^{2n}$, one has
\begin{equation}\label{structure-zero}
\begin{cases}
 \lambda' |\xi|^{p-2} |\eta|^2 \le  \p_{\xi_j} A_i(x,\xi) \eta_i \eta_j \le \Lambda' |\xi|^{p-2} |\eta|^2,
 \\
|A_i(x,\xi)| +  |\p_{x_j} A_i(x,\xi)| \le  \Lambda' |\xi|^{p-1}.
\end{cases}
\end{equation}

Given an open set $\Omega\subset \mathbb H^n$, we indicate with $W^{1,p}(\Omega)$ the Sobolev space associated with the $p$-energy $\mathscr E_{\Omega,p}(u) = \frac 1p \int_{\Omega}|\nabla_0 u|^p$, i.e., the space of all functions $u\in L^p(\Omega)$ such that their distributional derivatives $X_i u$, $i=1,...,2n,$ are also in $L^p(\Omega)$. The corresponding norm
is $||u||^p_{W^{1,p}(\Omega)} = ||u||_{L^p(\Omega)}+||\nabla_0 u||_{L^p(\Omega)}.$ We
will add the subscript {\it loc} for the local versions of such spaces, and denote by $W^{1,p}_0(\Omega)$ the completion of $C^\infty_0(\Omega)$ with respect to such norm.  A function $u\in L^p((0,T), W_{{\rm loc}}^{1,p}(\Omega))$
is a weak solution of \eqref{maineq} in the cylinder $\Omega\times (0,T)$ if
\begin{equation}\label{weak}
\int_0^T \int_\Omega u \phi_t -  \sum_{i=1}^{2n} A_i(x,\nabla_0 u) X_i \phi  =0,
\end{equation}
for every $\phi\in C^{\infty}_0(Q)$.  Our main result is a Lipschitz regularity estimate for weak solutions, on parabolic cylinders $Q_{\mu,r}$ (see Definition \ref{cylinder}).

\begin{thrm}\label{main1}
Let $A_i$ satisfy the structure conditions \eqref{structure-zero} and
let $u\in L^p((0,T), W_{{\rm loc}}^{1,p}(\Omega))$ be a weak solution of \eqref{maineq-zero}  in $Q=\Omega\times (0,T)$.  If $2\le p \le 4$ then $|\nabla_0 u| \in L^\infty_{{\rm loc}}(Q)$ and $\partial_t u, Zu \in L^q_{{\rm loc}}(Q)$ for every $1\le q<\infty$. Moreover, one has
that  for any $Q_{\mu, 2r}\subset Q$,
\begin{equation}\label{desired}
\sup_{Q_{\mu, r}}  |\nabla_0u|
\le C  \max\Big( \Big(\frac{1}{\mu r^{N+2}}\int\int_{Q_{\mu,2r}} (\delta+|\nabla_0u|^2)^{\frac{p}{2}}\Big)^{\frac{1}{2}}, \mu^{\frac{p}{2(2-p)}}\Big),
\end{equation}
where $C=C(n,p,\lambda, \Lambda, r,\mu)>0$.  In the special case where there is no direct dependence on the space variable, i.e.  $A_i(x,\xi)=A_i(\xi)$, the parameters dependence  is more explicit, with
$C=C(n,p,\lambda, \Lambda) \mu^{\frac 1 2}>0$.
\end{thrm}

The boundedness of the gradient of solutions to equation (\ref{maineq-zero}) in the setting of Euclidean spaces is well-known \cite{DF, Ch, DB}. The above theorem shows that this is also the case in the setting of Heisenberg group for the range $2\le p\le 4$. We believe that it is true for all range $1<p<\infty$ if the solution is bounded as in the Euclidean setting.
The proof of the main theorem is based on the Caccioppoli type estimate in Proposition \ref{cor1}, which is new in the subelliptic context. The Lipschitz regularity follows  then through a Moser type iteration, which we present in detail in Section \ref{lipschitz estimate}. The Caccioppoli estimate in Proposition \ref{cor1} is derived using two main ingredients: The first of these  is an approximation scheme, that allows us to invoke the regularity results from \cite{CCG}, thus dealing with smooth approximants that can be differentiated directly, without recurring to fractional difference quotients. One of the  original contributions of the present paper is that we can avoid the extra assumption of Riemannian approximation which is needed in \cite{CCG} (hypotheses (1.6) and (1.7) in \cite{CCG}), and in fact we prove that our structure hypotheses \eqref{structure-zero} imply that such approximation always exists. The second ingredient is a Poincar\'e-type inequality for smooth functions, which was originally established in \cite{CLM}. This latter estimate, in Lemma \ref{SobXU}, is the only point in the paper where we are forced to impose the limited range $2\le p \le 4$. In fact, we believe that such constraint is not needed for the Caccioppoli inequality in Proposition \ref{cor1} and we plan to return to this point in future work.

To the best of our knowledge, the present paper and \cite{CCG} are the first instances in the literature of the study of higher regularity for weak solutions of the non stationary $p$-Laplacian PDE   in the sub-Riemannian setting. Both are based on  techniques introduced by Zhong in \cite{Zhong}. By contrast, the stationary case has been far more developed, mostly  in the Heisenberg group case. We mention here the contributions of  Domokos \cite{Dom}, Manfredi, Mingione \cite{Manfredi-Mingione}, Mingione, Zatorska-Goldstein and Zhong \cite{Mingione-Zatorska-Zhong}, Ricciotti \cite{Ricciotti}, \cite{Ricciotti1} and eventually those in  \cite{Zhong,MZ}, where   the horizontal $C^{1,\alpha}$ regularity in the full range $1<p<\infty$ is proved. Regularity in more general  contact sub-Riemannian manifolds, including the rototraslation group,  has been recently established by  two of the authors and coauthors \cite{CCLDO} and independently by Mukherjee
\cite{M} based on an extension of the techniques in \cite{Zhong}. More recently Domokos and Manfredi \cite{DomMan} have studied regularity in higher steps groups and in some special non-group structures.

The structure of the paper is as follows: In Section \ref{preliminaries} we review some preliminary definitions and results from \cite{CCG} and lay out the approximation scheme, thus reducing the problem to finding  estimates for smooth solutions $u_\delta$ of approximating regularized equations, which are stable as $\delta\to 0$. From that point on, we will simplify the notation by dropping the script $\delta$ and by focusing on the case $A_i(x,\xi)=A_i(x)$, thus highlighting how in this case we can obtain more explicit constants in the right hand side of our estimates. In Section \ref{old energy estimates} we recall some energy type estimates from \cite{CCG}.  In Section \ref{new energy estimates} we show that   derivatives of weak solutions along the center are in every $L^q_{loc}$ space, $q\ge 2$, uniformly in $\delta>0$ and establish the key Caccioppoli type inequality in Proposition \ref{cor1}. In this section we need to use the limitation $2\le p \le 4$ in the proof of a Sobolev type estimate. We conjecture that with the exception of  Lemma \ref{SobXU}, all other estimates continue to hold in the range $2\le p<\infty$. Using  Proposition \ref{cor1}, in Section \ref{lipschitz estimate} we prove that the solutions are locally Lipschitz continuous in the subRiemannian metric (i.e. the horizontal gradient is in $L^\infty_{loc}$) uniformly in $\delta$.
We note explicitly that the Moser iteration in Section
\ref{lipschitz estimate} involves a Sobolev type estimate, which is also stable as $\delta\to 0$, in view of the results in \cite{CC}. Section \ref{time-derivative} addresses the higher integrability of the time derivatives $\p_t u$ of weak solutions.

\noindent {\it Acknowledgements} The authors are grateful to Nicola Garofalo for  conversations around the topics of this paper. Indeed, the results presented here are a development of  some initial work that the authors did jointly with him.

\section{Approximating weak solutions via regularizations}\label{preliminaries}

As mentioned in the introduction, our strategy for the proof of the Lipschitz regularity is to locally
approximate the weak solutions of \eqref{maineq-zero} with smooth solutions  $u_\delta$ of less degenerate PDE such as \eqref{maineq} and prove  estimates on such approximate solutions that are uniform as $\delta\to 0$.  This approximation is built using both the regularity results in \cite{CCG},  recalled below in Theorem \ref{CCG}, and a Riemannian approximation scheme (see \cite{Montgomery,CC} and references therein). We start by recalling the main points of the latter.
First, we will use interchangeably the notation $Z$ and $X_{2n+1}$ for the generator of the center of the Lie algebra.  The left invariant sub-Riemannian metric $(\mathbb H^n,g_0)$  defined by $\langle X_i, X_j \rangle_0 =\delta_{ij}$, for $i,j=1,...,2n$, can be approximated in the Gromov-Hausdorff sense through a sequence of Riemannian metrics $g_\epsilon$, for $\epsilon\to 0^+$, defined by imposing that
$X_1,...,X_{2n},\epsilon Z$ is an orthonormal $g_\epsilon-$frame for all $\epsilon>0$. In the terminology of \cite{Montgomery}, the metrics $g_\epsilon$ {\it tame} the metric $g_0$.  We relabel the vectors in this frame as $X_1^\epsilon, ..., X_{2n+1}^\epsilon$. The corresponding gradient
\[\nabla_\epsilon u= \sum_{i=1}^{2n} X_i u X_i + \epsilon^2 Zu Z= \sum_{i=1}^{2n+1} X_i^\epsilon u X_i^\epsilon
\]
has the obvious property that
$\nabla_\epsilon u \to (\nabla_0 u, 0)$ as $\epsilon \to 0$.  We note explicitly that
\[| \nabla_\epsilon u|_{\epsilon}^2:=| \nabla_\epsilon u |_{g_\epsilon}^2=\sum_{i=1}^{2n} (X_i u)^2+ \epsilon^2 (Zu)^2\to |\nabla_0 u|_0,
\]
as $\epsilon\to 0$. For $\delta>0$, the $\delta-$regularized Riemannian $p-$Laplacian, i.e. the operator related to the Euler-Lagrange equations for the $p$-energy $\int |\nabla_\epsilon u|_\varepsilon^p dx$, is
\begin{equation}\label{epsilon-plap} L_p^\epsilon u := \sum_{i=1}^{2n+1} X_i^\epsilon (  [\delta+| \nabla_\epsilon u|_\epsilon^2]^{\frac{p-2} {2}}X_i^\epsilon u),\end{equation}
and provides a  natural (quasilinear) elliptic regularization of the subelliptic $p$-Laplacian.

Next, we recall the regularity theorem proved in \cite{CCG}.
\begin{thrm}[\cite{CCG}]\label{CCG}  For $\Omega\subset \mathbb H^n$, $2\le p < \infty$, and  $\delta> 0$, assume that the
the functions $A_{i,\delta}:\Omega\times \mathbb R^{2n}\to \mathbb R$, $i=1,...,2n$ satisfy the following structure conditions:

\noindent {\bf (i)} For some $\lambda, \Lambda>0$ depending only on $\lambda', \Lambda'$, one has
\begin{equation}\label{structure1}
\begin{cases}
\lambda (1+|\xi|^2)^{\frac{p-2}{2}} |\eta|^2 \le  \p_{\xi_j} A_{i, \delta} (x,\xi) \eta_i \eta_j \le \Lambda (1+|\xi|^2)^{\frac{p-2}{2}} |\eta|^2,
\\
|A_i (x,\xi)| +  |\p_{x_j} A_{i, \delta} (x,\xi)| \le  \Lambda (1+|\xi|^2)^{\frac{p-1}{2}}.
\end{cases}
\end{equation}
\noindent {\bf (ii)} We assume that one can approximate $A_{i,\delta}$ by  a $1$-parameter family of regularized approximants $A^\epsilon_\delta(x,\xi)=(A_{1,\delta}^\epsilon(x,\xi),...,A_{2n+1,\delta}^\epsilon(x,\xi))$  defined for a.e. $x\in \Omega$ and every $\xi \in \mathbb R^{2n+1}$, and such that
 for a.e. $x\in\Omega, $  for all $\xi =  \sum_{i=1}^{2n} \xi_i X_i^\epsilon + \xi_{2n+1} X_{2n+1}^\epsilon$ , and $\xi^\epsilon =  \sum_{i=1}^{2n} \xi_i X_i^\epsilon + \epsilon \xi_{2n+1} X_{2n+1}^\epsilon$   one has uniformly on compact subsets of $\Omega$,
\begin{equation}\label{structure2}
(A_{1,\delta}^\epsilon(x,\xi^\epsilon),...,A_{2n+1,,\delta}^\epsilon(x,\xi^\epsilon))\ \underset{\epsilon\to 0^+}{\longrightarrow}\    (A_{1,,\delta}(x, \xi_1,...,\xi_{2n}), ..., A_{2n,\delta} (x, \xi_1,...,\xi_{2n}),0),
\end{equation}
and furthermore
\begin{equation}\label{structure-epsilon}
\begin{cases}
\lambda (1+|\xi|^2)^{\frac{p-2}{2}} |\eta|^2 \le  \p_{\xi_j} A_{i,\delta}^\epsilon(x,\xi) \eta_i \eta_j \le \Lambda (1+|\xi|^2)^{\frac{p-2}{2}} |\eta|^2,
\\
|A_{i,\delta}^\epsilon(x,\xi)| +  |\p_{x_j} A_{i,\delta}^\epsilon(x,\xi)| \le  \Lambda (1+|\xi|^2)^{\frac{p-1}{2}},
\end{cases}
\end{equation}
 for all $\eta\in \mathbb R^{2n+1}$, and
 for some $0<\lambda\le \Lambda <\infty$ independent of $\epsilon$.
 Let $u_\delta\in L^p((0,T), W_{{\rm loc}}^{1,p}(\Omega))$  be a weak solution of  \begin{equation}\label{maineq}
\p_t u_\delta = \sum_{i=1}^{2n} X_i A_{i, \delta}(x, \nabla_0 u_\delta). \end{equation} in $Q=\Omega\times (0,T)$.
 If $\delta>0$ then $u_\delta$  is $C^\infty$ smooth in $Q$.
 \end{thrm}

One of our main contributions in the present paper is that one can avoid the assumptions \eqref{structure2}, and \eqref{structure-epsilon}, and build the Riemannian approximation using solely the structure condition \eqref{structure1}. Our result is stated in the following proposition

\begin{prop}\label{approximation theorem}
Let $A_i$ be as in \eqref{structure-zero}. For every $\delta>0$ there exists $A_{i,\delta}$ such that
\begin{equation}\label{approx-structure}
A_{\delta} (x,\xi)\ \underset{\delta\to 0^+}{\longrightarrow}\    A(x, \xi),
\end{equation}
satisfying the hypothesis \eqref{structure1}, \eqref{structure2}, and \eqref{structure-epsilon} with constants depending only on the original $\lambda',\Lambda'$.  Moreover, if a function $A_{i,\delta}$ satisfies \eqref{structure1}, then it also satisfies \eqref{structure2}, and \eqref{structure-epsilon} with constants depending only on the original $\lambda^\prime,\Lambda^\prime$.
\end{prop}

In view of Theorem \cite{CCG}, the latter yields immediately the following
\begin{cor}\label{approx}
Let $u$ be a weak solution of \eqref{maineq-zero} in $Q=\Omega\times (0,T)$, with the structure conditions \eqref{structure-zero}. For any sub-cylinder $Q_1=\Omega_1\times (t_1,t_2) \subset \subset  \Omega\times (0,T)$,  there exists a  sequence $\{u_{\delta}\}$  of smooth solutions of the
 regularized problem
\begin{equation}\label{maineq}
\p_t u_\delta = \sum_{i=1}^{2n} X_i A_{i, \delta}(x, \nabla_0 u_\delta) \quad \quad  \text{ in  \quad }Q_1, \text{ and } u_\delta=u \quad \text{ on }\partial_p  \ Q_1\end{equation}
converging to $u$, as $\delta\to 0^+$,  uniformly on compacts subsets of $Q_1$ and weakly in the $W^{1,p}$-norm.  Here we have denoted by $\partial_p Q_1=\Omega_1\times \{t=t_1\}\cup \partial \Omega_1 \times (t_1,t_2)$ the parabolic boundary of $Q_1$. The functions $A_\delta$ satisfy \eqref{approx-structure}, \eqref{structure1}, \eqref{structure2}, and \eqref{structure-epsilon} with constants depending only on the original $\lambda',\Lambda'$.
\end{cor}
\begin{proof}
Let $u$ be a weak solution of  \eqref{maineq-zero} in $Q$. In view of the results in \cite{ACCN} we know that the
solution is H\"older continuous in compact subsets.  For $\delta>0$,  let $A_{i,\delta}$ be as in the statement of Proposition \ref{approximation theorem}, and consider  the unique weak  solution
$u_{\delta}$ of \eqref{maineq}.  In view of the comparison principle, the uniform continuity  and Caccioppoli inequalities for $\{u_\delta\}$ proved in \cite{ACCN} and of \eqref{approx-structure}, \eqref{structure1} one can easily see that $u_\delta \to u$ uniformly on compact subsets of $Q_1$ and weakly
in the $W^{1,p}$-norm. In order to conclude the proof, we  only need to observe that $u_\delta$ are smooth (with the regularity possibly depending on $\delta>0$ of course) thanks to
Theorem \ref{approximation theorem} and Theorem \ref{CCG}.
In fact,  one can apply the results from \cite{CCG}, to derive regularity estimates that are uniform in the parameter $\epsilon$, thus  yielding that the family $u_{\delta, \varepsilon}$ has a subsequence converging to a solution of the problem
\eqref{maineq} as $\varepsilon \to 0$, which coincides with $u_\delta$ in view the comparison principle.
\end{proof}

We are left with the task of proving Proposition \ref{approximation theorem}. To better illustrate the argument of the proof  we present the special, simpler,  case of the $p$-Laplacian \eqref{plap-para}, i.e. for $\xi  = \sum_{i=1}^{2n} \xi_i X_i\in \mathbb R^{2n}$, and $x\in \mathbb H^n$,
$$A_i(x,\xi)=|\xi|^{p-2}\xi_i, \quad \quad i=1,...,2n.$$
 In this case, we consider for each $\delta>0$, the  following functions
 \begin{equation}\label{plap-para-delta}
 A_{i,\delta}(x,\xi)= (\delta+|\xi|^2)^{\frac{p-2}{2}}\xi_i, \quad \quad i=1,...,2n,
 \end{equation}
 and for each $\epsilon>0$, and $\xi  = \sum_{i=1}^{2n+1} \xi_i X_i^\epsilon\in \mathbb R^{2n+1}$,
 \begin{equation}\label{plap-para-delta-epsilon }
  A_{i,\delta}^\epsilon(x,\xi)= (\delta+\|\xi\|_{g_\epsilon(x)}^2)^{\frac{p-2}{2}}\xi_i,\quad \quad i=1,...,2n+1.
 \end{equation}
 While the quantity $||\cdot||_{g_\epsilon(x)}$ a priori depends on $x\in \mathbb H^n$, we remark that when $\xi$ is a left invariant vector field, since $g_\epsilon$ is left invariant as well, the dependence of $||\xi||_{g_\epsilon(x)}$ on the point $x$ vanishes.

\begin{proof}[Proof of  Proposition \ref{approximation theorem}]
Following the intuition from the example above, we construct the approximates through a two steps process.
For $0<\delta<1$, let us define
 \begin{equation}\label{approx1}
 A_{\delta} (x,\xi)= A (x,\xi) + \lambda\delta ^{\frac{p-2}{2}} \xi
 \end{equation}
It is clear that
\begin{equation}\label{approx-structure}
A_{\delta} (x,\xi)\ \underset{\delta\to 0^+}{\longrightarrow}\    A(x, \xi)
\end{equation}
and furthermore, for some $\lambda, \Lambda>0$ depending only on $\lambda', \Lambda'$, one has the estimate \eqref{structure1}.

For each  $\xi  = \sum_{i=1}^{2n+1} \xi_i X_i^\epsilon\in \mathbb R^{2n+1}$, and $\epsilon,\delta>0$ we set
\begin{equation}
\label{approx2} A_{i,\delta, \varepsilon} (x,\xi) =  \tilde A_{i}(x,\xi_H) + \lambda  (\delta + |\xi|_\epsilon^2)^{\frac{p-2}{2}} \xi_i,
\end{equation}
for $i=1,...,2N+1$.
 Here we have denoted $\xi_H=(\xi_1,...,\xi_{2n})$,  $\tilde A=(A,0)\in \mathbb R^{2n+1}$, and
$|\xi|_\epsilon^2= \sum_{i=1}^{2n+1} \xi_i^2 $.

Clearly
 for a.e. $x\in\Omega, $ and for all $\xi^\epsilon =  \sum_{i=1}^{2n} \xi_i X_i^\epsilon + \epsilon \xi_{2n+1} X_{2n+1}^\epsilon$  one has uniformly on compact subsets of $\Omega \times (0,T)$,
$$
(A_{1, \delta, \varepsilon} (x,\xi^\epsilon),...,A_{2n+1, \delta, \varepsilon} (x,\xi^\epsilon))\ \underset{\varepsilon\to 0^+}{\longrightarrow}\    (A_{1, \delta}(x, \xi_1,...,\xi_{2n}), ..., A_{2n, \delta} (x, \xi_1,...,\xi_{2n}),0),
$$
where $A_\delta$ is defined as in \eqref{approx1}.
In addition one can see that there exist constants $\lambda, \Lambda>0$ depending on $\lambda^\prime, \Lambda^\prime$ such that
$$
\begin{cases}
\lambda (1+|\xi|_\epsilon^2)^{\frac{p-2}{2}} |\eta|_\epsilon^2 \le  \p_{\xi_j} A_{i, \delta, \epsilon} (x,\xi) \eta_i \eta_j \le \Lambda (1+|\xi|_\epsilon^2)^{\frac{p-2}{2}} |\eta|_\epsilon^2,
\\
|A_i (x,\xi)| +  |\p_{x_j} A_{i, \delta,\epsilon} (x,\xi)| \le  \Lambda (1+|\xi|_\epsilon^2)^{\frac{p-1}{2}},
\end{cases}
$$
 for all $\eta =\sum_{i=1}^{2n+1} \eta_i X_i^\epsilon\in \mathbb R^{2n+1}$.
\end{proof}

\begin{rmrk}\label{notation}
For the rest of the paper we will always consider solutions $u_\delta$ of the Dirichlet problem \eqref{maineq} with $\delta>0$,  in a cylinder $ D\times (\tau_1,\tau_2)\subset \subset Q_1$, with $D\subset \subset \Omega_1$ and $[\tau_1,\tau_2]\subset (t_1,t_2)$. For the sake of notation we will drop the subscript $\delta$ from $u_\delta$ and $A_{i,\delta}$, and with a slight abuse of notation write  $Q=\Omega \times (0,T)$ instead of $D\times (\tau_1,\tau_2)$. To further simplify the formulation of the estimates, we will assume that $A_i(x,\xi)= A_i(\xi)$, as in this case we can obtain sharper constants, and so we highlight these more involved aspects of the proofs. The more general case is handled in a similar fashion, and does not lead to explicit constants on the right hand side of the estimates.  We note explicitly that all constants are independent of the parameter $\delta>0$.
\end{rmrk}

\section{Preliminary energy estimates}\label{old energy estimates}
We recall the basic Caccioppoli inequalities proved in \cite{CCG}. These inequalities apply to a smooth solution $u$ of the approximating equation \eqref{maineq} with $\delta>0$,  in a cylinder $ Q\subset \subset Q_1$. In what follows we will implicitly assume that all constants on the right hand side of the inequalities depend on $n, p$, on the structure constants, $\lambda, \Lambda$ {\it but not on $\delta$.}

\begin{lemma} \label{diff-PDE} Let $u$ be a solution of \eqref{maineq} in $Q$, with $\delta>0$. If we set $v _l = X _lu$, with $ l = 1, 2, . . . , 2n$, and $s_l = (-1)^{[l/n]}$ then the function $v _l$ is  a  solution of
\begin{equation}\label{eqderivX}
\p_t v _l
= \sum_{i,j=1}^{2n} X _i\Big( A_{i, \xi_j}  (\nabla_0u) X _l X _ju \Big) + s_l Z (A_{l+s_ln}  (\nabla_0u)) .\end{equation}

\end{lemma}

\begin{lemma}\label{eqZ} Let $u$ be a solution of \eqref{maineq} in  $Q$, with $\delta>0$. The function
$Zu$  is  then a  solution of the equation
$$\partial_t Zu=
\sum_{i,j=1}^{2n}
X _i(A _{i, \xi_j}(\nabla_0u)X _j Zu).
$$
\end{lemma}

First we recall a Caccioppoli estimates for derivatives of the solution along the center of the group.

\begin{lemma}[Lemma 3.4 \cite{CCG}]\label{stimaZu} Let $u$ be a solution of \eqref{maineq} in $Q$ with $\delta>0$.
For every $ \beta \geq 0$ and  non-negative $\eta\in C^{1}([0,T], C^\infty_0 (\Omega))$ vanishing on the parabolic boundary of $Q$, one has
\begin{equation*}
\begin{aligned}
\int_{t_1} ^{t_2}\int_\Omega
(\delta+|\nabla_0u|^2)^{\frac{p-2}{2}}|Zu|^{\beta} | \nabla_0Zu|^2\eta^{4+\beta}
\le & C
\int_{t_1} ^{t_2}\int_\Omega
(\delta+|\nabla_0u|^2)^{\frac{p-2}{2}}|Zu|^{\beta+2}|\nabla_0\eta|^2 \eta^{2+\beta}\\
&
+ C\int_{t_1}^{t_2}\int_\Omega |Zu|^{\beta+2} |\partial_t \eta |  \eta^{3+\beta},
\end{aligned}
\end{equation*}
where $C=C(\lambda, \Lambda)>0$.
\end{lemma}

Second, we recall a Caccioppoli estimate for the horizontal derivatives.

\begin{lemma}[Lemma 3.5 \cite{CCG}] \label{lemma3.4} Let $u$ be a weak solution of \eqref{maineq} in $Q$, with $\delta>0$. For every  $\beta \geq 0 $ and non-negative $\eta\in C^{1}([0,T], C^\infty_0 (\Omega))$ vanishing on the parabolic boundary of $Q$,
we have
\begin{equation*}
\begin{aligned}
&\frac{1}{\beta+2}\sup_{t_1<t<t_2}\int_\Omega (\delta+ |\nabla_0u|^2)^{\frac{\beta+2}{2}}\eta^2 +\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{(p-2+\beta)/2} |\nabla_0^2u|^2 \eta^2\\
\leq & \, C
\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{\frac{p+\beta}{2}} ( |\nabla_0\eta|^2 + |Z\eta| \eta)
+
\frac{C}{\beta+2}
\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{\frac{\beta+2}{2}}|\partial_t \eta| \eta \\
& +
C(\beta + 1)^4
\int_{t_1} ^{t_2} \int_\Omega
(
\delta + |\nabla_0u|^2)^{\frac{p-2+\beta}{
2}} |Zu|^2 \eta^2,
\end{aligned}
\end{equation*}
where $C = C(n, p, \lambda, \Lambda) > 0$, independent of $\delta$.
\end{lemma}

\section{Main Caccioppoli inequality}\label{new energy estimates}

The main result of this section is a Caccioppoli inequality, Proposition \ref{cor1}, for the horizontal derivatives of the weak solutions of \eqref{maineq}, with $\delta>0$. To do this,  we first need to prove an  estimate for  the derivative along the center $Zu$ in Lemma \ref{SobXU} and Lemma \ref{mainlemma}.  All estimates are  uniform in $\delta>0$, and the constants are stable as $\delta\to 0$.

We begin by recalling a Poincar\'e-like interpolation inequality from \cite{CLM}. In the proof, we will need the restriction $2\le p\le 4$ and this is the only use we make of this hypothesis in the paper.

\begin{lemma}\label{SobXU} Assume that $2\le p \leq 4$ and let $u\in C^2(Q)$.
There exists a constant $C>0$ depending only on $n, p$ such that for every $ \beta \geq 0$ and non-negative $\eta\in C^{1}([0,T], C^\infty_0 (\Omega))$ vanishing on the parabolic boundary of $Q$, we have
\begin{equation}\label{interpolation}
\begin{aligned}
\int_{t_1} ^{t_2}\int_\Omega &
|Zu|^{p+\beta} \eta^{p+\beta}
\le  C(p+\beta) ||\nabla_0\eta||_{L^\infty}\int\int_{spt(\eta)} (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}}\\
&  +
C(p+\beta)\int_{t_1} ^{t_2}\int_\Omega
(\delta + |\nabla_0u|^2)^{(p-2)/2}
|Zu|^{\beta} | \nabla_0Zu|^2\eta^{4+\beta}
\end{aligned}
\end{equation}

\end{lemma}
\begin{proof}
We denote
\begin{equation}\label{LR}
L=\int_{t_1} ^{t_2}\int_\Omega
|Zu|^{p+\beta} \eta^{p+\beta}, \quad R=\int\int_{spt(\eta)} (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}}.
\end{equation}
We estimate $L$ from above as follows. Fix $l=1,2,...,n$. Note that
\[
Zu=X_lX_{n+l}u-X_{n+l}X_lu.
\]
We can write
\[ |Zu|^{p+\beta}=|Zu|^{p-2+\beta}Zu (X_lX_{n+l}u-X_{n+l}X_lu).
\]
Then integration by parts gives us
\begin{equation}\label{ML1}
\begin{aligned}
L=&\int_{t_1} ^{t_2}\int_\Omega |Zu|^{p-2+\beta}Zu (X_lX_{n+l}u-X_{n+l}X_lu)\eta^{p+\beta}\\
=&-(p-1+\beta)\int_{t_1} ^{t_2}\int_\Omega |Zu|^{p-2+\beta}(X_lZuX_{n+l}u-X_{n+l}ZuX_lu)\eta^{p+\beta}\\
&-(p+\beta)\int_{t_1} ^{t_2}\int_\Omega|Zu|^{p-2+\beta}Zu(X_{n+l}uX_l\eta-X_luX_{n+l}\eta)\eta^{p-1+\beta}\\
\le & \, 2(p+\beta)\int_{t_1} ^{t_2}\int_\Omega |\nabla_0u||Zu|^{p-2+\beta}|\nabla_0Zu|\eta^{p+\beta}\\
&+2(p+\beta)\int_{t_1} ^{t_2}\int_\Omega |\nabla_0u||Zu|^{p-1+\beta}|\nabla_0\eta|\eta^{p-1+\beta}=I_1+I_2
\end{aligned}
\end{equation}
We will estimate the integrals $I_1,I_2$ on the right hand side of \eqref{ML1} by H\"older's inequality. First for $I_2$, we have
\begin{equation}\label{ML2}
I_2\le 2(p+\beta) ||\nabla_0 \eta||_{L^\infty} R^{\frac{1}{p+\beta}}L^{\frac{p-1+\beta}{p+\beta}},
\end{equation}
where $L$ and $R$ are as in \eqref{LR}.

Second, for $I_1$, we have
\begin{equation}\label{ML3}
I_1\le 2(p+\beta) M^{\frac{1}{2}}R^{\frac{4-p}{2(p+\beta)}} L^{\frac{2p-4+\beta}{2(p+\beta)}},
\end{equation}
where
\[
M=\int_{t_1} ^{t_2}\int_\Omega
|\nabla_0u|^{p-2}|Zu|^{\beta} | \nabla_0Zu|^2\eta^{4+\beta}.
\]

This yields
\begin{equation}\label{last-eq}
L\le C(p+\beta) ||\nabla_0 \eta||_{L^\infty} R^{\frac{1}{p+\beta}}L^{\frac{p-1+\beta}{p+\beta}}
+C (p+\beta) M^{\frac{1}{2}}R^{\frac{4-p}{2(p+\beta)}} L^{\frac{2p-4+\beta}{2(p+\beta)}},
\end{equation}
from which the conclusion follows immediately through Young's inequality.
\end{proof}

The previous Poincar\'e-like inequality can be applied to solutions of \eqref{maineq} and through invoking  Lemma \ref{stimaZu} lead us to the following key estimate.

\begin{lemma}\label{mainlemma}
Let $u$ be a solution of \eqref{maineq} in $Q$, with $\delta>0$ and $2\leq p \leq 4$.
Then for every $ \beta \geq 0$ and non-negative $\eta\in C^{1}([0,T], C^\infty_0 (\Omega))$ vanishing on the parabolic boundary of $Q$, we have
\begin{equation}\label{mainestimate}
\begin{aligned}
\Big(\int_{t_1} ^{t_2}\int_\Omega &
|Zu|^{p+\beta} \eta^{p+\beta} \Big)^{\frac{1}{p+\beta}}
\le  C(p+\beta) ||\nabla_0\eta||_{L^\infty}\Big(\int\int_{spt(\eta)} (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}}\Big)^{\frac{1}{p+\beta}}\\
&  +
C(p+\beta)||\eta\partial_t\eta||_{L^\infty}^{\frac{1}{2}}|spt(\eta)|^{\frac{p-2}{2(p+\beta)}}\Big(\int\int_{spt(\eta)}  (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}} \Big)^{\frac{4-p}{2(p+\beta)}}
\end{aligned}
\end{equation}
\end{lemma}
\begin{rmrk}\label{limit-remark}  Suspending temporarily the notation established in Remark \ref{notation}, we denote by $u_\delta$ the solutions of the approximating equation \eqref{maineq}.  In particular, Lemma \ref{mainlemma} establishes the local $L^q$  integrability of  $Z u_\delta$, the derivative along the center of the approximating solutions, with uniform $L^q$ bounds as $\delta\to 0$, for all $1\le q<\infty$. This implies that  one can find a subsequence, $Z u_{\delta_k}$  converging to a $L^q_{{\rm loc}}$ function, which in view of the definition of weak derivative, is also a derivative along the center of the uniform limit of the $u_\delta$. Since such limit is the original solution of \eqref{maineq-zero}, this proves the local  integrability of $Zu$ in every $L^q$ class as stated in Theorem \ref{main1}.
\end{rmrk}

\begin{proof}
We apply the inequality \eqref{last-eq} in the previous lemma to the solution $u$ and invoke Lemma \ref{stimaZu} to estimate the integral
\[
M=\int_{t_1} ^{t_2}\int_\Omega
|\nabla_0u|^{p-2}|Zu|^{\beta} | \nabla_0Zu|^2\eta^{4+\beta},
\]
  obtaining
\begin{equation}\label{ML4}
\begin{aligned}
M &\le C
\int_{t_1} ^{t_2}\int_\Omega
(\delta+|\nabla_0u|^2)^{\frac{p-2}{2}}|Zu|^{\beta+2}|\nabla_0\eta|^2 \eta^{2+\beta}
+ C\int_{t_1}^{t_2}\int_\Omega |Zu|^{\beta+2} |\partial_t \eta |  \eta^{3+\beta}\\
&\le C ||\nabla_0 \eta||_{L^\infty}^2 R^{\frac{p-2}{p+\beta}}L^{\frac{\beta+2}{p+\beta}}
+C||\eta\partial_t\eta||_{L^\infty}|spt(\eta)|^{\frac{p-2}{p+\beta}}L^{\frac{\beta+2}{p+\beta}},
\end{aligned}
\end{equation}
where $C=C(\lambda, \Lambda)>0$, and $L$ is as in \eqref{ML1}. In the second inequality of \eqref{ML4}, we used H\"older's inequality.
Combining \eqref{ML3}, \eqref{last-eq} and \eqref{ML4}, we obtain the estimate for $I_1$,
\begin{equation}\label{ML5}
I_1\le C(p+\beta) ||\nabla_0 \eta||_{L^\infty} R^{\frac{1}{p+\beta}}L^{\frac{p-1+\beta}{p+\beta}}
+C (p+\beta) ||\eta\partial_t\eta||_{L^\infty}^{\frac 1 2}|spt(\eta)|^{\frac{p-2}{2(p+\beta)}} R^{\frac{4-p}{2(p+\beta)}}
L^{\frac{p-1+\beta}{2(p+\beta)}}.
\end{equation}
Next, we substitute the latter in  the estimate \eqref{ML3} for $I_1$ and the estimate \eqref{ML2} for $I_2$ to \eqref{ML1}, and conclude
\[
L\le C(p+\beta) ||\nabla_0 \eta||_{L^\infty} R^{\frac{1}{p+\beta}}L^{\frac{p-1+\beta}{p+\beta}}
+C (p+\beta) ||\eta\partial_t\eta||_{L^\infty}^{\frac 1 2}|spt(\eta)|^{\frac{p-2}{2(p+\beta)}} R^{\frac{4-p}{2(p+\beta)}}
L^{\frac{p-1+\beta}{2(p+\beta)}},
\]
which yields immediately  \eqref{mainestimate}. \end{proof}

The following result follows from  Lemma \ref{mainlemma}, and the energy estimate in Lemma \ref{lemma3.4}. It yields a Caccippoli inequality for the horizontal derivatives of weak solutions, which extends to the subRiemannian setting the analogue  Euclidean estimate   proved in \cite[Proposition 3.2 (3.7), page 225]{DB}.

\begin{prop}\label{cor1}
Let $u$ be a weak solution of \eqref{maineq} in $Q$, with $\delta>0$ and $2\leq p \leq 4$.
Then for every $ \beta \geq 0$ and non-negative $\eta\in C^{1}([0,T], C^\infty_0 (\Omega))$ vanishing on the parabolic boundary of $Q$, we have
\begin{equation}\label{maincacci}
\begin{aligned}
&\sup_{t_1<t<t_2}\int_\Omega (\delta+ |\nabla_0u|^2)^{\frac{\beta+2}{2}}\eta^2 +\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{(p-2+\beta)/2} |\nabla_0^2u|^2 \eta^2\\
\leq & C(p+\beta)^7\big( ||\nabla_0\eta||_{L^\infty}^2+||\eta Z\eta||_{L^\infty}\big)\int\int_{spt(\eta)} (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}}\\
& +
C(p+\beta)^7||\eta\partial_t\eta||_{L^\infty}|spt(\eta)|^{\frac{p-2}{p+\beta}}\Big(\int\int_{spt(\eta)}  (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}} \Big)^{\frac{\beta+2}{p+\beta}},
\end{aligned}
\end{equation}
where $C=C(n,p,\lambda,\Lambda)>0$.
\end{prop}
\begin{rmrk} Although the statement addresses the approximate solution $u_\delta$, in view of   arguments analogue to those in  Remark \ref{limit-remark}, the same estimate holds for weak solutions of \eqref{maineq-zero}.
\end{rmrk}
\begin{proof}
Lemma \ref{lemma3.4} gives us the following estimate for the left hand side of \eqref{maincacci}
\begin{equation}\label{NH1}
\begin{aligned}
&\sup_{t_1<t<t_2}\int_\Omega (\delta+ |\nabla_0u|^2)^{\frac{\beta+2}{2}}\eta^2 +\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{(p-2+\beta)/2} |\nabla_0^2u|^2 \eta^2\\
\leq & \, C (p+\beta)
\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{\frac{p+\beta}{2}} ( |\nabla_0\eta|^2 + |Z\eta| \eta)
+
C
\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{\frac{\beta+2}{2}}|\partial_t \eta| \eta \\
& +
C(p+\beta )^5
\int_{t_1} ^{t_2} \int_\Omega
(
\delta + |\nabla_0u|^2)^{\frac{p-2+\beta}{
2}} |Zu|^2 \eta^2.
\end{aligned}
\end{equation}

To obtain the desired estimate \eqref{NH1}, we will show that each integral on the right hand side of \eqref{NH1} can be bounded from above by the right hand side of \eqref{maincacci}. For the first integral on the right hand of \eqref{NH1}, it is obviously bounded from above by the first item on the right hand side of \eqref{maincacci}. For the second integral, H\"older's inequality gives us
\[
\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{\frac{\beta+2}{2}}|\partial_t \eta| \eta
\le ||\eta\partial_t\eta||_{L^\infty} |spt(\eta)|^{\frac{p-2}{p+\beta}}\Big(\int\int_{spt(\eta)}  (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}} \Big)^{\frac{\beta+2}{p+\beta}},
\]
which shows that it is bounded from above by the second item on the right hand side of \eqref{maincacci}.

For the third integral on the right hand side of \eqref{NH1}, we use H\"older's inequality and our main lemma, Lemma \ref{mainlemma},  and we have
\begin{equation*}
\begin{aligned}
\int_{t_1} ^{t_2} \int_\Omega
(\delta + |\nabla_0u|^2)^{\frac{p-2+\beta}{2}} |Zu|^2 \eta^2 &\le \Big(\int\int_{spt(\eta)}  (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}}\Big)^{\frac{p-2+\beta}{p+\beta}}\Big(\int_{t_1} ^{t_2}\int_\Omega
|Zu|^{p+\beta} \eta^{p+\beta} \Big)^{\frac{2}{p+\beta}}\\
&\le C(p+\beta)^2 ||\nabla_0\eta||_{L^\infty}^2\int\int_{spt(\eta)} (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}}\\
&\ +
C(p+\beta)^2||\eta\partial_t\eta||_{L^\infty}|spt(\eta)|^{\frac{p-2}{p+\beta}}\Big(\int\int_{spt(\eta)}  (\delta+ |\nabla_0u|^2)^{\frac{p+\beta}{2}} \Big)^{\frac{\beta+2}{p+\beta}},
\end{aligned}
\end{equation*}
which concludes the proof of the lemma.
\end{proof}

\section{Boundedness of the horizontal gradient}\label{lipschitz estimate}
In this section we conclude the proof of Theorem \ref{main1}, i.e. we establish that weak solutions of the approximating equation \eqref{maineq} with $\delta>0$ are Lipschitz continuous with respect to the subRiemannian distance, uniformly in the parameter $\delta$.  The proof follows immediately from Proposition \ref{cor1} and from the Moser type iteration in Theorem \ref{main2} below. The proof of Theorem \ref{main2} should be known, but we can not find the precise reference in the literature. It is similar to the proof of Theorem 4 in \cite{Ch} for the case $1<p<2$. The proof is included for the reader's convenience.

First, we recall a few definitions needed in the proof. We will denote by  $d_0(x,y)=||y^{-1}x||$ the subRiemannian distance, where
$$||x||^4=(\sum_{i=1}^{2n} x_i^2)^2+ 16 x_{2n+1}^2,$$
is the Koranyi gauge.
The corresponding parabolic metric is $d_0 ((x,t), (y,s))= d_0(x,y)+|t-s|^2$.

\begin{dfn}\label{cylinder} A parabolic cylinder $Q_{r}(x_0, t_0)\subset Q$
is a set of the form
$Q_{r}(x_0, t_0)= B(x_0, r) \times (t_0-r^2, t_0 ).$
where $r>0$, $B(x_0, r)=\{y|\ ||yx_0^{-1}||<r\} \subset \Omega$ denotes the gauge ball of center $x_0$.
The parabolic boundary of the cylinder $Q_{ r}(x_0, t_0)\subset Q$ is
the  set $B(x_0, r) \times \{t_0-r^2\} \cup \partial B(x_0, r) \times [t_0-r^2, t_0).$
For $r$, $\mu>0$ we also define the cylinders
$$ Q_{\mu, r} : = B(x, R) \times [t_0 - \mu R, t_0] $$
\end{dfn}

\begin{thrm}\label{main2}
Let $u\in C^\infty(Q)$, with $Q=\Omega\times (0,T)$. If $u$ satisfies the Caccioppoli type inequality \eqref{maincacci}, then  for every $p\ge 2$, and  for any $Q_{\mu, 2r}\subset Q$, we have
\begin{equation}\label{desired}
\sup_{Q_{\mu, r}}  |\nabla_0u|
\le C \mu^{\frac 1 2} \max\Big( \Big(\frac{1}{\mu r^{N+2}}\int\int_{Q_{\mu,2r}} (\delta+|\nabla_0u|^2)^{\frac{p}{2}}\Big)^{\frac{1}{2}}, \mu^{\frac{p}{2(2-p)}}\Big),
\end{equation}
where $C=C(n,p,\lambda, \Lambda)>0$.
\end{thrm}
\begin{rmrk} Suspending temporarily the notation established in Remark \ref{notation}, we denote by $u_\delta$ the solutions of the approximating equation \eqref{maineq}.  As mentioned earlier, there is a subsequence $u_\delta\to u$ converging uniformly in compact subsets of $Q$ to the  weak solution $u$ of \eqref{maineq-zero}. In view of the uniform  bound on the Lipschitz constant of $u_\delta$ in \eqref{desired}, then the Lipschitz regularity of $u$ follows immediately.
\end{rmrk}
\begin{proof}
Let  $\eta\in C^{1}([0,T], C^\infty_0 (\Omega))$ be a non-negative cut-off function vanishing on the parabolic boundary of $Q$ such that
$|\eta|\le 1$ in $Q$. For $\beta\ge 0$, we set
\[
v=(\delta+|\nabla_0u|^2)^{\frac{p+\beta}{4}}\eta^2.
\]
Then the Caccioppoli inequality \eqref{maincacci} gives us
\begin{equation}\label{maincacci2}
\begin{aligned}
\sup_{t_1<t<t_2}\int_\Omega v^m +\int_{t_1} ^{t_2} \int_\Omega
|\nabla v|^2
\leq &  C(p+\beta)^7\big( ||\nabla_0\eta||_{L^\infty}^2+||\eta Z\eta||_{L^\infty}\big)\int\int_{spt(\eta)} v^2\\
&+
C(p+\beta)^7||\eta\partial_t\eta||_{L^\infty}|spt(\eta)|^{\frac{p-2}{p+\beta}}\Big(\int\int_{spt(\eta)} v^2\Big)^{\frac{\beta+2}{p+\beta}},
\end{aligned}
\end{equation}
where $C=C(n,p,\lambda,\Lambda)>0$. Here $m=2(\beta+2)/(p+\beta)$. Note that $4/p<m \le 2$. Now let $q=2(m+N)/N$. We have
\[
\int_{t_1} ^{t_2} \int_\Omega v^q\le  \int_{t_1}^{t_2}\Big(
\int_\Omega v^m\Big)^{\frac 2 N}\Big(\int_\Omega v^{\frac{2N}{N-2}}\Big)^{\frac{N-2}{N}}\le C \Big(\sup_{t_1<t<t_2}\int_\Omega v^m\Big)^{\frac 2 N} \Big(\int_{t_1} ^{t_2} \int_\Omega
|\nabla v|^2\Big),
\]
where $C=C(n)>0$. Here in the second inequality, we used the Sobolev inequality in the space variables. Now we plug the estimate \eqref{maincacci2} into the above inequality and we obtain that
\begin{equation}\label{MM1}
\begin{aligned}
\Big(\int_{t_1} ^{t_2} \int_\Omega v^q\Big)^{\frac{N}{N+2}}\le &  C(p+\beta)^7\big( ||\nabla_0\eta||_{L^\infty}^2+||\eta Z\eta||_{L^\infty}\big)\int\int_{spt(\eta)} v^2\\
&+
C(p+\beta)^7||\eta\partial_t\eta||_{L^\infty}|spt(\eta)|^{\frac{p-2}{p+\beta}}\Big(\int\int_{spt(\eta)} v^2\Big)^{\frac{\beta+2}{p+\beta}},
\end{aligned}
\end{equation}
where $C=C(n,p,\lambda,\Lambda)>0$.  Here $q=2+4(\beta+2)/(N(p+\beta))$. This is the inequality on which our iteration is based.

Let $Q_{\mu, 2r}\subset Q$. We define, for $i=0,1,2,...$, a sequence of radius
$r_i=(1+2^{-i})r$ and a  sequence of exponent $\beta_i$ such  that $\beta_0=0$ and
\[
p+\beta_{i+1}=(p+\beta_i) \big(1+\frac{2(\beta_i+2)}{N(p+\beta_i)}\big),
\]
that is,
\[
\beta_i=2(\kappa^i-1), \quad \kappa=\frac{N+2}{N}.
\]
We denote $Q_i=Q_{\mu, r_i}$. Note that $Q_0=Q_{\mu, 2r}$ and $Q_\infty=Q_{\mu, r}$. The we choose a standard parabolic cut-off function $\eta_i\in C^\infty(Q_i)$ such that $\eta_i=1$ in $Q_{i+1}$ with
\[
|\nabla_0 \eta_i|\le \frac{2^{i+8}}{r}, \quad |Z\eta_i|\le  \frac{2^{2i+8}}{r^2}, \quad |\partial_t \eta_i|\le \frac{2^{2i+8}}{\mu r^2} \quad \text{in }Q_i.
\]
Now we let $\eta=\eta_i$ and $\beta=\beta_i$ in \eqref{MM1} and we obtain that for $i=0,1,...$
\begin{equation}\label{MM2}
\begin{aligned}
\Big(\int\int_{Q_{i+1}} (\delta+|\nabla_0u|^2)^{\frac{\alpha_{i+1}}{2}}\Big)^{\frac{N}{N+2}}\le & C2^{2i}\alpha_i^7r^{-2}
\Big[ \Big(\int\int_{Q_i}(\delta+|\nabla_0u|^2)^{\frac{\alpha_{i}}{2}}\Big)^{\frac{p-2}{\alpha_i}}+\\
&+
\mu^{-1}\big(\mu r^{N+2}\big)^{\frac{p-2}{\alpha_i}}
\Big] \Big( \int\int_{Q_i}  (\delta+|\nabla_0u|^2)^{\frac{\alpha_{i}}{2}}\Big)^{\frac{\alpha_i-p+2}{\alpha_i}},
\end{aligned}
\end{equation}
where $C=C(n,p,\lambda,\Lambda)>0$ and $\alpha_i=p+\beta_i=p-2+2\kappa^i$. We denote
\[
M_i=\Big(\frac{1}{\mu r^{N+2}}\int\int_{Q_i} (\delta+|\nabla_0u|^2)^{\frac{\alpha_{i}}{2}}\Big)^{\frac{1}{\alpha_i}}.
\]
Then  we can write \eqref{MM2} as
\[
M_{i+1}^{\frac{\alpha_{i+1}}{\kappa}}\le C \mu^{\frac{2}{N+2}}2^{2i}\alpha_i^7\big(M_i^{p-2}+\mu^{-1}\big)M_i^{\alpha_i-p+2}.
\]
We set
\[
{\overline M}_i=\max\big( M_i, \mu^{\frac{1}{2-p}}\big).
\]
Then it follows from the above inequality that
\begin{equation}\label{iter}
{\overline M}_{i+1}^{\frac{\alpha_{i+1}}{\kappa}}\le C \mu^{\frac{2}{N+2}}2^{2i}\alpha_i^7\overline{M}_i^{\alpha_i},
\end{equation}
since we may assume that $C=C(n,p, \lambda,\Lambda)\ge 1$.
Iterating (\ref{iter}), we obtain that
\[
{\overline M}_{i+1}\le \Big(\prod_{j=0}^i K_j^{\frac{\kappa^{i+1-j}}{\alpha_{i+1}}}\Big){\overline M}_0^{\frac{\alpha_0\kappa^{i+1}}{\alpha_{i+1}}},
\]
where
\[
K_i=C \mu^{\frac{2}{N+2}}2^{2i}\alpha_i^7.
\]
Recall that $\alpha_i=p-2+2\kappa^i$ and $\kappa= (N+2)/N$. Let $i$ go to infinity. We obtain that
\begin{equation}\label{MM3}
{\overline M}_\infty=\limsup_{i\to \infty} {\overline M}_i\le C \mu^{\frac{1}{2}}{\overline M}_0^{\frac p 2},
\end{equation}
where $C=(n,p,\lambda,\Lambda)>0$. Note that
\[
{\overline M}_\infty \ge \sup_{Q_{\mu,r}}|\nabla_0 u|, \quad {\overline M}_0=\max \Big( \Big(\frac{1}{\mu r^{N+2}}\int\int_{Q_{\mu,2r}} (\delta+|\nabla_0u|^2)^{\frac{p}{2}}\Big)^{\frac{1}{p}}, \mu^{\frac{1}{2-p}}\Big).
\]
Thus \eqref{MM3} gives us the desired inequality \eqref{desired}, completing the proof.
\end{proof}

\section{Higher integrability of $\partial_t u$}\label{time-derivative}

In this section, we prove that the time derivative $\partial_tu$ of weak solutions of \eqref{maineq} in the range $2\le p\le 4$ belongs to
$L^q_{{\rm loc}}(\Omega\times (0,T))$ for every $q\ge 1$.  As observed before, once we  establish uniform estimates for $\p_t u_\delta$,  then arguing as in Remark \ref{limit-remark},  one can readily  conclude the integrability of $\p_t u$.

\begin{lemma} Let $u$ be a solution of equation \eqref{maineq} in $Q=\Omega\times (0,T)$. Then we have
\[
\partial_t u\in L^q_{{\rm loc}}(\Omega\times (0,T))
\]
for every $q\ge 1$. Moreover, for every $\beta\ge 0$, and non-negative $\eta\in C^1([0,T], C^\infty_0(\Omega))$ vanishing on the
parabolic boundary, we have
\begin{equation}\label{tderivative}
\int_{t_1}^{t_2}\int_\Omega|\partial_tu|^{\beta+2}\eta^{\beta+2}\le  C\big(M^{2p-2}||\nabla_0\eta||_{L^\infty}^2+M^p || \eta\partial_t\eta||_{L^\infty}\big)^{\frac{\beta+2}{2}}\vert spt(\eta)\vert,
\end{equation}
where $C=C(p,\lambda, \Lambda, \beta)>0$ and $M=\sup_{spt(\eta)}(\delta+|\nabla_0u|^2)^{\frac{1}{2}}$.
\end{lemma}

\begin{proof}
Let $\beta\ge 0$. Since $u$ is a solution of \eqref{maineq}, we can write
\[
|\partial_tu|^{\beta+2}=|\partial_tu|^\beta \partial_tuX_i(A_{i,\delta}(x,\nabla_0u)).
\]
We denote by $L$ the integral on the left hand side of \eqref{tderivative}, which is the object we will estimate. Let $\eta\in C^1([0,T], C^\infty_0(\Omega))$
be a non-negative cut-off function, vanishing on the
parabolic boundary. Since $\eta(\cdot, t)\in C^\infty_0(\Omega)$ for every $t\in [0,T]$,  integration by parts gives us that
\begin{equation}\label{test1}
\begin{aligned}
L=&\int_{t_1}^{t_2}\int_\Omega|\partial_tu|^{\beta+2}\eta^{\beta+2}=\int_{t_1}^{t_2}\int_\Omega|\partial_tu|^{\beta}\partial_tuX_i(A_{i,\delta}(x,\nabla_0u))\eta^{\beta+2}\\
=&-(\beta+2)\int_{t_1}^{t_2}\int_\Omega|\partial_tu|^\beta\partial_tuA_{i,\delta}(x,\nabla_0u)\eta^{\beta+1}X_i\eta\\
&-(\beta+1)\int_{t_1}^{t_2}\int_\Omega|\partial_tu|^\beta X_i(\partial_tu)A_{i,\delta}(x,\nabla_0u)\eta^{\beta+2}=I_1+I_2.
\end{aligned}
\end{equation}
We will estimate the integrals $I_1,I_2$ in the right hand side of the above equality as follows. First, we use the structure condition and H\"older's inequality to estimate $I_1$. We have
\begin{equation}\label{test2}
\begin{aligned}
|I_1|\le & (\beta+2)\Lambda\int_{t_1}^{t_2}\int_\Omega (\delta+|\nabla_0u|^2)^{\frac{p-1}{2}}|\partial_tu|^{\beta+1}\eta^{\beta+1}|\nabla_0\eta|\\
\le & (\beta+2)\Lambda \Big( \int_{t_1}^{t_2}\int_\Omega |\partial_tu|^{\beta+2}\eta^{\beta+2}\Big)^{\frac{\beta+1}{\beta+2}}\times\\
& \times\Big(\int_{t_1}^{t_2}\int_\Omega
(\delta+|\nabla_0u|^2)^{\frac{p-1}{2}(\beta+2)}|\nabla_0\eta|^{\beta+2}\Big)^{\frac{1}{\beta+2}}\\
=&(\beta+2)\Lambda
||\nabla_0\eta||_{L^\infty} |spt(\eta)|^{\frac{1}{\beta+2}} M^{p-1} L^{\frac{\beta+1}{\beta+2}},
\end{aligned}
\end{equation}
where $M=\sup_{spt(\eta)}(\delta+|\nabla_0u|^2)^{\frac{1}{2}}$.

Second, we also use the structure condition and H\"older's inequality to estimate $I_2$.
We have
\begin{equation}\label{test3}
\begin{aligned}
|I_2|\le & (\beta+1)\Lambda\int_{t_1}^{t_2}\int_\Omega  (\delta+|\nabla_0u|^2)^{\frac{p-1}{2}}|\partial_tu|^{\beta}|\nabla_0\partial_tu|\eta^{\beta+2}\\
\le & (\beta+1)\Lambda
\Big(\int_{t_1}^{t_2}\int_\Omega |\partial_tu|^{\beta+2}\eta^{\beta+2}\Big)^{\frac{\beta}{2(\beta+2)}}\Big(\int\int_{spt(\eta)}(\delta+|\nabla_0u|^2)^{\frac{p}{4}(\beta+2)}\Big)^{\frac{1}{\beta+2}} J^{\frac 1 2}\\
\le & (\beta+1)\Lambda  |spt(\eta)|^{\frac{1}{\beta+2}}
M^{\frac{p}{2}}L^{\frac{\beta}{2(\beta+2)}}J^{\frac 1 2},
\end{aligned}
\end{equation}
where
\[
J=\int_{t_1}^{t_2}\int_\Omega (\delta+|\nabla_0u|^2)^{\frac{p-2}{2}}|\partial_tu|^\beta|\nabla_0\partial_tu|^2\eta^{\beta+4}.
\]
To estimate the integral $J$, we differentiate equation \eqref{maineq} with respect to $t$ and we obtain that
\[
\partial_t(\partial_tu)=X_i(\partial_{\xi_j}A_{i,\delta}(x,\nabla_0u)X_j(\partial_tu)).
\]
Then we use $\varphi=|\partial_tu|^\beta\partial_tu\eta^{\beta+4}$ as a test function to the above equation and we obtain the following
Caccioppoli inequality by the structure condition and the Cauchy-Schwarz inequality
\begin{equation}\label{test4}
\begin{aligned}
J\le &C \int_{t_1}^{t_2}\int_\Omega (\delta+|\nabla_0u|^2)^{\frac{p-2}{2}}|\partial_tu|^{\beta+2}\eta^{\beta+2}|\nabla_0\eta|^2
+C \int_{t_1}^{t_2}\int_\Omega|\partial_tu|^{\beta+2}\eta^{\beta+3}|\partial_t\eta|\\
\le & C (M^{p-2} ||\nabla_0\eta||_{L^\infty}^2+||\eta\partial_t\eta||_{L^\infty})L,
\end{aligned}
\end{equation}
where $C=C(p,\lambda,\Lambda,\beta)>0$. Here we used the fact that $\eta$ vanishes on the parabolic boundary.
Combining \eqref{test3} and \eqref{test4}, we obtain the following estimate for $I_2$.
\begin{equation}\label{test5}
|I_2|\le C  |spt(\eta)|^{\frac{1}{\beta+2}}
M^{\frac{p}{2}}L^{\frac{\beta+1}{\beta+2}}(M^{p-2} ||\nabla_0\eta||_{L^\infty}^2+||\eta\partial_t\eta||_{L^\infty})^{\frac 1 2}.
\end{equation}
Now we combine \eqref{test1} with the estimates \eqref{test2} and \eqref{test5} and we end up with
\[
L\le C M^{p-1}||\nabla_0\eta||_{L^\infty} |spt(\eta)|^{\frac{1}{\beta+2}} L^{\frac{\beta+1}{\beta+2}}+C  |spt(\eta)|^{\frac{1}{\beta+2}}
M^{\frac{p}{2}}L^{\frac{\beta+1}{\beta+2}}(M^{p-2} ||\nabla_0\eta||_{L^\infty}^2+||\eta\partial_t\eta||_{L^\infty})^{\frac 1 2},
\]
from which \eqref{tderivative} follows. This completes the proof.
\end{proof}

\section{Concluding remarks and some open problems}

There are a number of immediate extensions which we want to highlight, as well as some more involved,  plausible extensions which we listen as open problems.

First of all, the prototype for the class of operators in \eqref{maineq-zero}  is the regularized $p$-Laplacian operator $$L_p u= \text{div}_{g_0,\mu_0}( (1+|\nabla_0 u|_{g_0}^2)^{\frac{p-2}{2}} \nabla_0 u)$$
in a sub-Riemannian contact manifold $(M,\omega,g_0)$, where $M$ is the underlying differentiable manifold, $\omega$ is the contact form and $g_0$ is a Riemannian metric on the contact distribution. The measure $\mu_0$ is the corresponding Popp measure.
Since the structure conditions \eqref{structure-zero} and equation \eqref{maineq-zero} are invariant by contact diffeomorphisms, then invoking Darboux coordinates one can pull-back the PDE from the setting of contact subRiemannian manifolds to that of the Heisenberg group. Consequently all our results extends to the more general contact subRiemannian setting.
For a more detailed description, see \cite[Section 6.1]{CCLDO}.  As an immediate corollary of Theorem \ref{main1} one has the following.

\begin{thrm}\label{main3}  Let $(M,\omega,g_0)$ be a contact, sub-Riemannian manifold and let $\Omega\subset M$ be an open set. For $2\le p\le 4$,  $\delta\ge 0$,
consider $u\in L^p((0,T), W_0^{1,p}(\Omega))$ be a weak solution of
$$\p_t u =\operatorname{div}_{g_0,\mu_0}( (1+|\nabla_0 u|_{g_0}^2)^{\frac{p-2}{2}} \nabla_0 u) ,$$ in $Q=\Omega\times (0,T)$.
 For any open ball $B\subset \subset \Omega$ and $T>t_2\ge t_1\ge 0$,  and $q\ge 1$,
there exist  constants $C=C(n, p, d(B, \p\Omega), T-t_2, 1)>0$  and $C_q=C(n, p,q, d(B, \p\Omega), T-t_2, 1)>0$ such that
$$
||\nabla_0 u||_{L^\infty(B\times (t_1,t_2))} \le C \text{ and } ||\p_t u||_{L^q(B\times (t_1,t_2))} + ||Z u||_{L^q(B\times (t_1,t_2))} \le C_q.$$
\end{thrm}
Of course, if $\delta>0$ then in view of the results in \cite{CCG}, the solutions are smooth in $Q$.

Some of the following extensions seem   challenging, and we list them as open problems in increasing order of their perceived difficulty.

\begin{enumerate}
\item Standard, but technically involved, modifications should allow to extend our work to the case of equations of the type
$$\p_t u - X_i A_i (x,t, u, \nabla_0 u) = B(x,t,u,\nabla_0 u)$$
with structure conditions similar to those in \cite[Section 1, Chapter VIII]{DB}.
\item We  feel it should possible to weaken the bounds in the structure conditions for $\p_{x_k} A_i$ and request instead  only horizontal derivatives bounds, bounds  on $X_k A_i$, although this would require some additional work in the proof of Lemma \ref{SobXU}.
\item This paper only deals with scalar equations, however in the Euclidean case the results continue to hold also for systems of equations with additional structure (see \cite{DB}). The extension in the subelliptic setting would involve first extending the results of \cite{CCG}, and all the regularity theory literature that is used there.
\item Because our argument rests in a crucial way on Lemma \ref{SobXU}, the Lipschitz regularity for the range $4< p <\infty$ is currently beyond our reach. We conjecture that our main Caccioppoli inequality \eqref{maincacci} still holds with exactly the same statement in this extended range.
\item Proof of the  H\"older regularity of horizontal derivatives, in any range of $p\neq 2$.
\item Just as in the Euclidean case, the regularity problem in the range $1<p<2$ is more challenging, and would require completely different arguments. In the stationary case this has been solved by Mukherjee and Zhong in \cite{MZ}.
\item Our work extends easily to any step two Carnot group. Although there is promising work by Domokos and Manfredi \cite{DomMan} in the stationary case, the extension of our result to the higher step setting seems quite challenging.
\end{enumerate}


\section{Summary}\label{summary}

\textsc{lenstronomy} is an Astropy-affiliated (\hyperlink{astropy:2013}{Astropy Collaboration 2013}, \hyperlink{astropy:2018}{2018}) Python package for gravitational lensing simulations and analyses.
\textsc{lenstronomy} was introduced by \hyperlink{lenstronomy1}{Birrer \& Amara (2018)} and is based on the linear basis set approach by \hyperlink{Birrer:2015}{Birrer et a. (2015)}.
The user and developer base of \textsc{lenstronomy} has substantially grown since then, and the software has become an integral part of a wide range of recent analyses, such as measuring the Hubble constant with time-delay strong lensing or constraining the nature of dark matter from resolved and unresolved small scale lensing distortion statistics.
The modular design has allowed the community to incorporate innovative new methods, as well as to develop enhanced software and wrappers with more specific aims on top of the \textsc{lenstronomy} API.
Through community engagement and involvement, \textsc{lenstronomy} has become a foundation of an ecosystem of affiliated packages extending the original scope of the software and proving its robustness and applicability
at the forefront of the strong gravitational lensing community in an open source and reproducible manner.

\begin{figure}
\centering
\includegraphics{paper_fig.png}
\caption{Illustration of the strong gravitational lensing phenomenology and the capabilities of \textsc{lenstronomy} in performing realistic simulations as well as reconstructing lensing and source properties from a given data set. Top row from left to right along the green arrow:
A galaxy is lensed around a foreground massive object, becomes highly distorted, and has components appearing multiple times. Observations of this phenomena are limited in resolution (convolution), depending on the detector (pixelation), and are subject to noise.
Bottom row from right to left along the red arrow: The inverse problem is solved with a linear basis set in the source morphology maximizing the likelihood of the model given the data.}
\label{fig:example}
\end{figure}

\section{Background}

Gravitational lensing displaces the observed positions and distorts the shapes of apparent objects on the sky due to intervening inhomogeneous matter along the line of sight. Strong gravitational lensing describes the regime where the background source, such as a galaxy or quasar, is lensed by a massive foreground object, such as another galaxy or cluster of galaxies, to produce multiple images of the source in a highly distorted manner.
The top row of Figure~\ref{fig:example} illustrates such a process from the intrinsic galaxy to the data product at hand, including the lensing distortions, effects of the instrument, observational conditions, and noise.

Analyses of strong gravitational lensing have provided a wealth of key insights into cosmology and astrophysics.
For example, relative time delays of multiply imaged variable sources provided precision measurements on the expansion rate of the Universe (\hyperlink{Wong:2020}{Wong et al., 2020}; \hyperlink{Shajib:2020strides}{Shajib et al., 2020}; \hyperlink{Birrer:2020tdcosmoiv}{Birrer et al., 2020}). Small scale distortions in the lensing signal of resolved sources
(\hyperlink{Vegetti:2012}{Vegetti et al., 2012}; \hyperlink{Hezaveh:2016}{Hezaveh et al., 2016}; \hyperlink{Birrer:2017}{Birrer et al., 2017})
and unresolved flux ratios
(\hyperlink{Gilman:2020}{Gilman et al. 2020}; \hyperlink{Hsueh:2020}{Hsueh et al. 2020}) constrain the nature of dark matter. Combined strong lensing and kinematic observables constrain the formation and evolution of galaxies
(\hyperlink{Sonnenfeld:2015}{Sonnenfeld et al., 2015}; \hyperlink{Shajib:2021slacs}{Shajib et al. 2021a})
, and the lensing magnification effect provides an otherwise inaccessible angle on the early Universe (\hyperlink{Zheng:2012}{Zheng et al. 2012}; \hyperlink{Cava:2018}{Cava et al. 2018}).

\section{Statement of need}\label{statement-of-need}

Strong lensing studies have significantly enhanced, and sometimes challenged, our current fundamental understanding of the Universe.
In the near future, with the onset of the next-generation ground and space-based wide and deep astronomical imaging
(Rubin, Roman, Euclid observatories; \hyperlink{LSST}{Ivezić et al., 2019}; \hyperlink{Roman}{Spergel et al., 2013}; \hyperlink{Euclid}{Laureijs et al., 2011} and interferometric
(SKA; \hyperlink{SKA}{Dewdney et al., 2009}) surveys, the number of discovered lenses of different types will be growing by more than an order of magnitude (\hyperlink{Collett:2015}{Collett 2015}; \hyperlink{OM10}{Oguri \& Marshall, 2010}).
Such large samples can provide unprecedented statistical precision to stress-test our current understanding and exploit discovery potential.
It is key that these demanding studies, at present and in the future, are conducted by reliable software and supported by reproducible and open-source analysis products to provide the most compelling and transparent evidence required to further our physical understanding.

The primary design goal of \textsc{lenstronomy} is to facilitate scientific investigations into the outstanding and most pressing questions in the cosmology and astrophysics community.
\textsc{lenstronomy} has been applied throughout its development to the most demanding modeling and inference problems in strong lensing and the software has evolved around the requirements of the scientific applications to facilitate robust analyses. The modular API of the original design of \textsc{lenstronomy} (\hyperlink{lenstronomy1}{Birrer \& Amara 2018}) has accommodated the addition of new features. Code review processes in the development phase have led to additional benefits for the user community at large beyond the specific needs of the developer.

\textsc{lenstronomy} provides reliable and well-tested specific functionalities, as well as top-level interfaces, which allow for adaptive and innovative usage in control by the scientific investigator.
Guidance for the user community is provided on multiple levels. First, source code is well documented and provided through
\href{http://lenstronomy.readthedocs.org}{readthedocs.org}. Second, a set of \textsc{jupyter} notebooks are provided in an \href{https://github.com/sibirrer/lenstronomy_extensions}{extension repository}. These notebooks demonstrate simplified example use cases, each notebook individually highlighting different specific functionalities of \textsc{lenstronomy}, including a \href{https://github.com/sibirrer/lenstronomy_extensions/blob/v1.8.1/lenstronomy_extensions/Notebooks/starting_guide.ipynb}{starting guide notebook}  to introduce the modular design structure of \textsc{lenstronomy}. Third, end-to-end analysis pipelines of some of the published work are publicly available, providing ‘real-life’ examples at advanced levels.

\section{Track-record of applications}

\textsc{lenstronomy} has been applied in and contributed to more than 30 peer reviewed publications since its first public release in 2018.
In particular, \textsc{lenstronomy} has been used to provide state-of-the-art measurements on real data sets, such as:
(i) Hubble constant measurements from three quadruly lensed quasars with Hubble Space Telescope (HST) imaging (\hyperlink{Birrer:2016}{Birrer et al., 2016}; \hyperlink{Birrer:2019}{Birrer et al., 2019}; \hyperlink{Shajib:2020strides}{Shajib et al., 2020}),
dynamical modeling in the hierarchical analysis by \hyperlink{Birrer:2020tdcosmoiv}{Birrer et al., (2020)}, and modeling of lensed supernovae (\hyperlink{Moertsell:2020}{Mörtsell et al. 2020});
(ii) inference of small scale dark matter properties from detailed studies of both, resolved imaging (\hyperlink{Birrer:2017}{Birrer et al., 2017}), and unresolved flux ratio statistics (\hyperlink{Gilman:2020}{Gilman et al., 2020});
(iii) decomposition of quasar and host galaxy light in both, lensed and unlensed cases (\hyperlink{Ding:2020}{Ding et al., 2020}; \hyperlink{Bennert:2021}{Bennert et al., 2021});
(iv) morphological studies of high-redshift sources in the cluster environment
(\hyperlink{Yang:2020}{Yang et al., 2020}; \hyperlink{Yang:2021}{2021});
(v) internal structure of galaxies (\hyperlink{Shajib:2021slacs}{Shajib et al., 2021a}; \hyperlink{Shajib:2021AO}{2021b});
(vi) measurements of the weak lensing effect imprinted in Einstein rings
(\hyperlink{Birrer:2017cosmos}{Birrer et al., 2017}; \hyperlink{Kuhn:2021}{Kuhn et al., 2021}).
Among the studies, some of them have applied a pipeline to uniformly analyse dozens of lenses of different types
(\hyperlink{Shajib:2019strides}{Shajib et al., 2019}; \hyperlink{Shajib:2021slacs}{2021a}; \hyperlink{Shajib:2021AO}{2021b}),
a milestone in moving towards utilizing thousands of lenses in the near future.

Beyond analyzing data, many theoretical studies have been conducted using \textsc{lenstronomy} to investigate statistical robustness in present and anticipated future analyses
(\hyperlink{BirrerTreu:2019}{Birrer \& Treu, 2019}; \hyperlink{Millon:2020}{Millon et al., 2020}; \hyperlink{vdVyvere:2020}{Van de Vyvere et al., 2020}; \hyperlink{Li:2021}{Li et al., 2021}; \hyperlink{Ding:2021transient}{Ding et al. 2021}),
as well as to provide forecasts for anticipated future constraints for different science cases (\hyperlink{Gilman:2019}{Gilman et al., 2019}; \hyperlink{Sengul:2020}{Çaǧan Şengül et al. 2020}; \hyperlink{BirrerTreu:2021}{Birrer \& Treu, 2021}).
Particularly, three separate teams participated in the blind time-delay lens modeling challenge (\hyperlink{Ding:2021tdlmc}{Ding et al. 2021}) using \textsc{lenstronomy}.

\textsc{lenstronomy} has seen a substantial development and incorporation of innovations and numerical recipes
(\hyperlink{Tessore:2015}{Tessore \& Metcalf, 2015}; \hyperlink{Shajib:2019unified}{Shajib, 2019}; \hyperlink{Joseph:2019}{Joseph et al., 2019}; \hyperlink{Galan:2021}{Galan et al., 2021}; \hyperlink{Birrer:2021arcs}{Birrer, 2021}),
and has found applications beyond its original aim due to the robust and high-standard design requirements.

\section{Ecosystem of affiliated packages}

\textsc{lenstronomy} has allowed the community to develop third-party analysis products and software products utilizing its core functionalities to provide more targeted and integrated software solutions for a wide range of scientific analyses.
These open-source \href{https://github.com/sibirrer/lenstronomy/blob/1.8.1/AFFILIATEDPACKAGES.rst}{affiliated packages} effectively create an ecosystem enhancing the capability of \textsc{lenstronomy}.
They provide specified and tested solution for specific scientific investigations, such as plug-ins and direct implementation for innovative source reconstruction algorithms (\href{https://github.com/aymgal/SLITronomy}{\textsc{SLITronomy}}; \hyperlink{Joseph:2019}{Joseph et al., 2019}; \hyperlink{Galan:2021}{Galan et al., 2021}),
gravitational wave lensing computations (\href{https://gitlab.com/gpagano/lensinggw}{\textsc{lensingGW}}; \hyperlink{Pagano:2020}{Pagano et al., 2020}),
automated pipelines for gravitational lensing reconstruction
(\href{https://github.com/ajshajib/dolphin}{\textsc{dolphin}}; \hyperlink{Shajib:2021slacs}{Shajib et al., 2021a}),
cluster source reconstruction and local perturbative lens modeling
(\href{https://github.com/ylilan/lenstruction}{\textsc{lenstruction}}; \hyperlink{Yang:2020}{Yang et al., 2020}),
enhancement in large-scale structure imaging survey simulations
(\href{https://github.com/LSSTDESC/SLSprinkler}{DESC \textsc{SLSprinkler}}; \hyperlink{LSSTDESC:2021}{Dark Energy Science Collaboration (LSST DESC) et al., 2021}),
rendering of sub-halos and line-of-sight halos
(\href{https://github.com/dangilman/pyHalo}{\textsc{pyHalo}}; \hyperlink{Gilman:2020}{Gilman et al., 2020}),
galaxy morphology analysis (\href{https://github.com/dartoon/galight}{\textsc{galight}}; \hyperlink{Ding:2020}{Ding et al., 2020}),
and hierarchical analyses to measure the Hubble constant
(\href{https://github.com/sibirrer/hierarc}{\textsc{hierArc}}; \hyperlink{Birrer:2020tdcosmoiv}{Birrer et al., 2020}).
With the rise in popularity and the promises in dealing with ever complex data problems with fast deep-learning methods,
dedicated tools for simulating large datasets for applying such methods to strong gravitational lensing
(\href{https://github.com/deepskies/deeplenstronomy}{\textsc{deeplenstronomy}}; \hyperlink{Morgan:2021}{Morgan et al., 2021}),
(\href{https://github.com/jiwoncpark/baobab}{\textsc{baobab}}; \hyperlink{Park:2021}{Park et al., 2021}),
as well as end-to-end Bayesian Neural Network training and validation packages for Hubble constant measurements
(\href{https://github.com/jiwoncpark/h0rton}{\textsc{h0rton}}; \hyperlink{Park:2021}{Park et al., 2021}),
and for a hierarchical analysis of galaxy-galaxy lenses
(\href{https://github.com/swagnercarena/ovejero}{\textsc{ovejero}}; \hyperlink{Wagner-Carena:2021}{Wagner-Carena et al., 2021}) have been developed.
The affiliated packages make best use of the \textsc{lenstronomy} modules without duplicating source code and make it possible to combine aspects of multiple affiliated packages in one single analysis.

\section{Related open source software}
\begin{itemize}
    \item \href{https://github.com/sibirrer/lenstronomy}{\textsc{lenstronomy}}(\hyperlink{Birrer:2015}{Birrer et al., 2015}; \hyperlink{lenstronomy1}{Birrer \& Amara 2018})
    \item \href{https://github.com/Jammy2211/PyAutoLens}{\textsc{PyAutoLens}}(\hyperlink{Nightingale:2018}{Nightingale et al., 2018}; \hyperlink{Nightingale:2021}{2021})
    \item \href{http://www.physics.rutgers.edu/~keeton/gravlens/}{\textsc{gravlens}}(\hyperlink{Keeton:2011}{Keeton, 2011})
    \item \href{https://www.slac.stanford.edu/~oguri/glafic/}{\textsc{glafic}}(\hyperlink{Oguri:2010}{Oguri, 2010})
    \item \href{https://github.com/jspilker/visilens}{\textsc{visilens}} (\hyperlink{Spilker:2016}{Spilker et al., 2016})
    \item \href{https://www.physik.uzh.ch/~psaha/lens/pixelens.php}{\textsc{PixeLens}} (\hyperlink{PixeLens}{Saha \& Williams, 2011})
    \item \href{https://github.com/j0r1/GRALE2}{\textsc{grale}} (\hyperlink{GRALE}{Liesenborgs et al., 2006})
    \item \href{http://projets.lam.fr/projects/lenstool/wiki}{\textsc{lenstool}} (\hyperlink{Jullo:2009}{Jullo \& Kneib, 2009})
\end{itemize}

\section{Acknowledgements}\label{acknowledgements}

Support for this work was provided by the National Science Foundation through NSF AST-1716527.
AJS was supported by NASA through the STScI grant HST-GO-15320 and by a Dissertation Year Fellowship from the UCLA Graduate Division.
This research was supported by the U.S. Department of Energy (DOE) Office of Science Distinguished Scientist Fellow Program.
DG is supported by NASA HST-GO-15177.
AG, MM LvdV, DS are supported by COSMICLENS: ERC grant agreement No 787886.
LT is supported by International Helmholtz-Weizmann Research School for Multimessenger Astronomy.
MU is supported by KIPAC and the Stanford Summer Research Program.
XD is supported by NASA HST-GO-15115.
TS is supported by NASA grant HST-GO-15320 and HST-GO-15652.
MZ is supported by the National Science Foundation of China.
AA is supported by a Royal Society Wolfson Fellowship.
We are grateful to the user community for valuable feedback and encouragement in continuing the development.

\section*{References}\label{references}
\addcontentsline{toc}{section}{References}

\hypertarget{refs}{}

\hypertarget{astropy:2013}{Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., Greenfield, P., Droettboom, M., Bray,
E., Aldcroft, T., Davis, M., Ginsburg, A., Price-Whelan, A. M., Kerzendorf, W. E., Conley,
A., Crighton, N., Barbary, K., Muna, D., Ferguson, H., Grollier, F., Parikh, M. M., Nair,
P. H., … Streicher, O. (2013). Astropy: A community Python package for astronomy.
\textit{Astronomy \& Astrophysics, 558, A33}. \href{https://doi.org/10.1051/0004-6361/201322068}{https://doi.org/10.1051/0004-6361/201322068}}

\hypertarget{astropy:2018}{
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., Lim, P. L.,
Crawford, S. M., Conseil, S., Shupe, D. L., Craig, M. W., Dencheva, N., Ginsburg, A.,
Vand erPlas, J. T., Bradley, L. D., Pérez-Suárez, D., de Val-Borro, M., Aldcroft, T. L.,
Cruz, K. L., Robitaille, T. P., Tollerud, E. J., … Astropy Contributors. (2018). The
Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package.
\textit{The Astronomical Journal, 156(3), 123}. \href{https://doi.org/10.3847/1538-3881/aabc4f}{https://doi.org/10.3847/1538-3881/aabc4f}}

\hypertarget{Bennert:2021}{
Bennert, V. N., Treu, T., Ding, X., Stomberg, I., Birrer, S., Snyder, T., Malkan, M. A.,
Stephens, A. W., \& Auger, M. W. (2021). A local baseline of the black hole mass scaling relations for active galaxies. IV. Correlations between $M_{\rm BH}$ and host galaxy $\sigma$, stellar
mass, and luminosity. \textit{arXiv e-Prints}, arXiv:2101.10355. \href{http://arxiv.org/abs/2101.10355}{http://arxiv.org/abs/2101.10355}}

\hypertarget{Birrer:2021arcs}{Birrer, S. (2021). Gravitational lensing formalism in a curved arc basis: A continuous description
of observables and degeneracies from the weak to the strong lensing regime. \textit{arXiv e-Prints}, arXiv:2104.09522. \href{http://arxiv.org/abs/2104.09522}{http://arxiv.org/abs/2104.09522}}

\hypertarget{lenstronomy1}{Birrer, S., \& Amara, A. (2018). lenstronomy: Multi-purpose gravitational lens modelling
software package. \textit{Physics of the Dark Universe, 22, 189–201}. \href{https://doi.org/10.1016/j.dark.2018.11.002}{https://doi.org/10.1016/j.dark.2018.11.002}}

\hypertarget{Birrer:2017}{
Birrer, S., Amara, A., \& Refregier, A. (2017). Lensing substructure quantification in RXJ1131-
1231: a 2 keV lower bound on dark matter thermal relic mass. \textit{Journal of Cosmology and
Astroparticle Physics, 2017(5), 037}. \href{https://doi.org/10.1088/1475-7516/2017/05/037}{https://doi.org/10.1088/1475-7516/2017/05/037}}

\hypertarget{Birrer:2016}{
Birrer, S., Amara, A., \& Refregier, A. (2016). The mass-sheet degeneracy and time-delay
cosmography: analysis of the strong lens RXJ1131-1231.
\textit{Journal of Cosmology and
Astroparticle Physics, 2016(8), 020}. \href{https://doi.org/10.1088/1475-7516/2016/08/020}{https://doi.org/10.1088/1475-7516/2016/08/020}}

\hypertarget{Birrer:2015}{
Birrer, S., Amara, A., \& Refregier, A. (2015). Gravitational Lens Modeling with Basis Sets.
\textit{The Astrophysical Journal, 813(2), 102}. \href{https://doi.org/10.1088/0004-637X/813/2/102}{https://doi.org/10.1088/0004-637X/813/2/102}}

\hypertarget{Birrer:2020tdcosmoiv}{
Birrer, S., Shajib, A. J., Galan, A., Millon, M., Treu, T., Agnello, A., Auger, M., Chen, G.
C.-F., Christensen, L., Collett, T., Courbin, F., Fassnacht, C. D., Koopmans, L. V. E.,
Marshall, P. J., Park, J.-W., Rusu, C. E., Sluse, D., Spiniello, C., Suyu, S. H., … Van de
Vyvere, L. (2020). TDCOSMO. IV. Hierarchical time-delay cosmography joint inference
of the Hubble constant and galaxy density profiles.
\textit{Astronomy \& Astrophysics, 643, A165}.
\href{https://doi.org/10.1051/0004-6361/202038861}{https://doi.org/10.1051/0004-6361/202038861}}

\hypertarget{BirrerTreu:2021}{
Birrer, S., \& Treu, T. (2021). TDCOSMO V: strategies for precise and accurate measurements of the Hubble constant with strong lensing.
\textit{Astronomy \& Astrophysics, 649, A61, 6}.
\href{https://doi.org/10.1051/0004-6361/202039179}{https://doi.org/10.1051/0004-6361/202039179}}

\hypertarget{BirrerTreu:2019}{
Birrer, S., \& Treu, T. (2019). Astrometric requirements for strong lensing time-delay cosmography.
\textit{Monthly Notices of the Royal Astronomical Society, 489(2), 2097–2103}.
\href{https://doi.org/10.1093/mnras/stz2254}{https://doi.org/10.1093/mnras/stz2254}}

\hypertarget{Birrer:2019}{
Birrer, S., Treu, T., Rusu, C. E., Bonvin, V., Fassnacht, C. D., Chan, J. H. H., Agnello,
A., Shajib, A. J., Chen, G. C.-F., Auger, M., Courbin, F., Hilbert, S., Sluse, D., Suyu,
S. H., Wong, K. C., Marshall, P., Lemaux, B. C., \& Meylan, G. (2019). H0LiCOW
- IX. Cosmographic analysis of the doubly imaged quasar SDSS 1206+4332 and a new measurement of the Hubble constant.
\textit{Monthly Notices of the Royal Astronomical Society,
484(4), 4726–4753}. \href{https://doi.org/10.1093/mnras/stz200}{https://doi.org/10.1093/mnras/stz200}}

\hypertarget{Birrer:2017cosmos}{
Birrer, S., Welschen, C., Amara, A., \& Refregier, A. (2017). Line-of-sight effects in strong lensing: putting theory into practice.
\textit{Journal of Cosmology and Astroparticle Physics,
2017(4), 049}. \href{https://doi.org/10.1088/1475-7516/2017/04/049}{https://doi.org/10.1088/1475-7516/2017/04/049}}

\hypertarget{Cava:2018}{
Cava, A., Schaerer, D., Richard, J., Pérez-González, P. G., Dessauges-Zavadsky, M., Mayer,
L., \& Tamburello, V. (2018). The nature of giant clumps in distant galaxies probed by the anatomy of the cosmic snake.
\textit{Nature Astronomy, 2, 76–82}. \href{https://doi.org/10.1038/
s41550-017-0295-x}{https://doi.org/10.1038/
s41550-017-0295-x}}

\hypertarget{Collett:2015}{
Collett, T. E. (2015). The Population of Galaxy-Galaxy Strong Lenses in Forthcoming Optical Imaging Surveys.
\textit{The Astrophysical Journal, 811(1), 20}. \href{https://doi.org/10.1088/0004-637X/811/1/20}{https://doi.org/10.1088/0004-637X/811/1/20}}

\hypertarget{Sengul:2020}{
Çaǧan Şengül, A., Tsang, A., Diaz Rivero, A., Dvorkin, C., Zhu, H.-M., \& Seljak, U. (2020). Quantifying the line-of-sight halo contribution to the dark matter convergence power spectrum from strong gravitational lenses. \textit{Physical Review D, 102(6), 063502}.
\href{https://doi.org/10.1103/PhysRevD.102.063502}{https://doi.org/10.1103/PhysRevD.102.063502}}

\hypertarget{SKA}{
Dewdney, P. E., Hall, P. J., Schilizzi, R. T.,
\& Lazio, T. J. L. W. (2009). The Square
Kilometre Array.
\textit{IEEE Proceedings, 97(8), 1482–1496}. \href{https://doi.org/10.1109/JPROC.2009.2021005}{https://doi.org/10.1109/JPROC.2009.2021005}}

\hypertarget{Ding:2021transient}{
Ding, X., Liao, K., Birrer, S., Shajib, A. J., Treu, T., \& Yang, L. (2021). Improved time delay lens modelling and H0 inference with transient sources. \textit{Monthly Notices of the Royal Astronomical Society 504(4), 5621–5628}. \href{https://doi.org/10.1093/mnras/stab1240}{https://doi.org/10.1093/mnras/stab1240}}

\hypertarget{Ding:2020}{
Ding, X., Silverman, J., Treu, T., Schulze, A., Schramm, M., Birrer, S., Park, D., Jahnke,
K., Bennert, V. N., Kartaltepe, J. S., Koekemoer, A. M., Malkan, M. A., \& Sanders, D. (2020). The Mass Relations between Supermassive Black Holes and Their Host Galaxies at 1 < z < 2 HST-WFC3.
\textit{The Astrophysical Journal, 888(1), 37}. \href{https://doi.org/10.3847/1538-4357/ab5b90}{https://doi.org/10.3847/1538-4357/ab5b90}}

\hypertarget{Ding:2021tdlmc}{
Ding, X., Treu, T., Birrer, S., Chen, G. C.-F., Coles, J., Denzel, P., Frigo, M., Galan, A.,
Marshall, P. J., Millon, M., More, A., Shajib, A. J., Sluse, D., Tak, H., Xu, D., Auger,
M. W., Bonvin, V., Chand, H., Courbin, F., … Williams, L. L. R. (2021). Time delay lens modelling challenge.
\textit{Monthly Notices of the Royal Astronomical Society, 503(1),
1096–1123}.
\href{https://doi.org/10.1093/mnras/stab484}{https://doi.org/10.1093/mnras/stab484}}

\hypertarget{Galan:2021}{
Galan, A., Peel, A., Joseph, R., Courbin, F., \& Starck, J.-L. (2021). SLITRONOMY: Towards a fully wavelet-based strong lensing inversion technique.
\textit{Astronomy \& Astrophysics, 647, A176}. \href{https://doi.org/10.1051/0004-6361/202039363}{https://doi.org/10.1051/0004-6361/202039363}}

\hypertarget{Gilman:2020}{
Gilman, D., Birrer, S., Nierenberg, A., Treu, T., Du, X., \& Benson, A. (2020). Warm dark matter chills out: constraints on the halo mass function and the free-streaming length of dark matter with eight quadruple-image strong gravitational lenses.
\textit{Monthly Notices of the
Royal Astronomical Society, 491(4), 6077–6101}. \href{https://doi.org/10.1093/mnras/stz3480}{https://doi.org/10.1093/mnras/stz3480}}

\hypertarget{Gilman:2019}{
Gilman, D., Birrer, S., Treu, T., Nierenberg, A., \& Benson, A. (2019). Probing dark matter
structure down to 107 solar masses: flux ratio statistics in gravitational lenses with line-of-sight haloes.
\textit{Monthly Notices of the Royal Astronomical Society, 487(4), 5721–5738}.
\href{https://doi.org/10.1093/mnras/stz1593}{https://doi.org/10.1093/mnras/stz1593}}

\hypertarget{Hezaveh:2016}{
Hezaveh, Y. D., Dalal, N., Marrone, D. P., Mao, Y.-Y., Morningstar, W., Wen, D., Blandford,
R. D., Carlstrom, J. E., Fassnacht, C. D., Holder, G. P., Kemball, A., Marshall, P. J., Murray, N., Perreault Levasseur, L., Vieira, J. D., \& Wechsler, R. H. (2016). Detection of Lensing Substructure Using ALMA Observations of the Dusty Galaxy SDP.81.
\textit{The Astrophysical Journal, 823(1), 37}. \href{https://doi.org/10.3847/0004-637X/823/1/37}{https://doi.org/10.3847/0004-637X/823/1/37}}

\hypertarget{Hsueh:2020}{
Hsueh, J.-W., Enzi, W., Vegetti, S., Auger, M. W., Fassnacht, C. D., Despali, G., Koopmans,
L. V. E., \& McKean, J. P. (2020). SHARP - VII. New constraints on the dark matter free-streaming properties and substructure abundance from gravitationally lensed quasars.
\textit{Monthly Notices of the Royal Astronomical Society, 492(2), 3047–3059}.
\href{https://doi.org/10.1093/mnras/stz3177}{https://doi.org/10.1093/mnras/stz3177}}

\hypertarget{LSST}{
Ivezić, Ž., Kahn, S. M., Tyson, J. A., Abel, B., Acosta, E., Allsman, R., Alonso, D., AlSayyad,
Y., Anderson, S. F., Andrew, J., Angel, J. R. P., Angeli, G. Z., Ansari, R., Antilogus,
P., Araujo, C., Armstrong, R., Arndt, K. T., Astier, P., Aubourg, É., … Zhan, H. (2019).
LSST: From Science Drivers to Reference Design and Anticipated Data Products.
\textit{The Astrophysical Journal, 873(2), 111}. \href{https://doi.org/10.3847/1538-4357/ab042c}{https://doi.org/10.3847/1538-4357/ab042c}}

\hypertarget{Joseph:2019}{
Joseph, R., Courbin, F., Starck, J.-L., \& Birrer, S. (2019). Sparse Lens Inversion Technique (SLIT): lens and source separability from linear inversion of the source reconstruction problem.
\textit{Astronomy \& Astrophysics, 623, A14}. \href{https://doi.org/10.1051/0004-6361/201731042}{https://doi.org/10.1051/0004-6361/201731042}}

\hypertarget{Jullo:2009}{
Jullo, E., \& Kneib, J.-P. (2009). Multiscale cluster lens mass mapping - I. Strong lensing modelling.
\textit{Monthly Notices of the Royal Astronomical Society 395(3), 1319–1332}. \href{https://doi.org/10.1111/j.1365-2966.2009.14654.x}{https://doi.org/10.1111/j.1365-2966.2009.14654.x}}

\hypertarget{Keeton:2011}{
Keeton, C. R. (2011). GRAVLENS: Computational Methods for Gravitational Lensing (ascl:1102.003).}

\hypertarget{Kuhn:2021}{
Kuhn, F. A., Birrer, S., Bruderer, C., Amara, A., \& Refregier, A. (2021). Combining strong and weak lensing estimates in the Cosmos field.
\textit{Journal of Cosmology and Astroparticle
Physics, 2021(4), 010}. \href{https://doi.org/10.1088/1475-7516/2021/04/010}{https://doi.org/10.1088/1475-7516/2021/04/010}}

\hypertarget{Euclid}{
Laureijs, R., Amiaux, J., Arduini, S., Auguères, J.-L., Brinchmann, J., Cole, R., Cropper, M.,
Dabin, C., Duvet, L., Ealet, A., Garilli, B., Gondoin, P., Guzzo, L., Hoar, J., Hoekstra, H.,
Holmes, R., Kitching, T., Maciaszek, T., Mellier, Y., … Zucca, E. (2011). Euclid Definition Study Report.
\textit{arXiv e-Prints}, arXiv:1110.3193.
\href{http://arxiv.org/abs/1110.3193}{http://arxiv.org/abs/1110.3193}}

\hypertarget{Li:2021}{
Li, N., Becker, C., \& Dye, S. (2021). The impact of line-of-sight structures on measuring H0 with strong lensing time-delays.
\textit{Monthly Notices of the Royal Astronomical Society}.
\href{https://doi.org/10.1093/mnras/stab984}{https://doi.org/10.1093/mnras/stab984}}

\hypertarget{GRALE}{
Liesenborgs, J., De Rijcke, S., \& Dejonghe, H. (2006). A genetic algorithm for the nonparametric inversion of strong lensing systems.
\textit{Monthly Notices of the Royal Astronomical Society, 367(3), 1209–1216}.
\href{https://doi.org/10.1111/j.1365-2966.2006.10040.x}{https://doi.org/10.1111/j.1365-2966.2006.10040.x}}

\hypertarget{LSSTDESC:2021}{
LSST Dark Energy Science Collaboration (LSST DESC), Abolfathi, B., Alonso, D., Armstrong,
R., Aubourg, É., Awan, H., Babuji, Y. N., Bauer, F. E., Bean, R., Beckett, G., Biswas,
R., Bogart, J. R., Boutigny, D., Chard, K., Chiang, J., Claver, C. F., Cohen-Tanugi, J.,
Combet, C., Connolly, A. J., … Zuntz, J. (2021). The LSST DESC DC2 Simulated Sky Survey.
\textit{The Astrophysical Journal Supplement Series, 253(1), 31}. \href{https://doi.org/10.3847/1538-4365/abd62c}{https://doi.org/10.3847/1538-4365/abd62c}}

\hypertarget{Millon:2020}{
Millon, M., Galan, A., Courbin, F., Treu, T., Suyu, S. H., Ding, X., Birrer, S., Chen, G.
C.-F., Shajib, A. J., Sluse, D., Wong, K. C., Agnello, A., Auger, M. W., Buckley-Geer,
E. J., Chan, J. H. H., Collett, T., Fassnacht, C. D., Hilbert, S., Koopmans, L. V. E., …
Van de Vyvere, L. (2020). TDCOSMO. I. An exploration of systematic uncertainties in the inference of H0 from time-delay cosmography.
\textit{Astronomy \& Astrophysics, 639, A101}.
\href{https://doi.org/10.1051/0004-6361/201937351}{https://doi.org/10.1051/0004-6361/201937351}}

\hypertarget{Morgan:2021}{
Morgan, R., Nord, B., Birrer, S., Lin, J., \& Poh, J. (2021). deeplenstronomy: A dataset simulation package for strong gravitational lensing.
\textit{The Journal of Open Source Software,
6(58), 2854}. \href{https://doi.org/10.21105/joss.02854}{https://doi.org/10.21105/joss.02854}}

\hypertarget{Moertsell:2020}{
Mörtsell, E., Johansson, J., Dhawan, S., Goobar, A., Amanullah, R., \& Goldstein, D. A.
(2020). Lens modelling of the strongly lensed Type Ia supernova iPTF16geu. \textit{Monthly Notices of the Royal Astronomical Society, 496(3), 3270–3280}. \href{https://doi.org/10.1093/mnras/staa1600}{https://doi.org/10.1093/mnras/staa1600}}

\hypertarget{Nightingale:2018}{
Nightingale, J. W., Dye, S., \& Massey, R. J. (2018). AutoLens: automated modeling of a strong lens’s light, mass, and source.
\textit{Monthly Notices of the Royal Astronomical Society, 478(4), 4738–4784}. \href{https://doi.org/10.1093/mnras/sty1264}{https://doi.org/10.1093/mnras/sty1264}}

\hypertarget{Nightingale:2021}{
Nightingale, J. W., Hayes, R., Kelly, A., Amvrosiadis, A., Etherington, A., He, Q., Li, N., Cao,
X., Frawley, J., Cole, S., Enia, A., Frenk, C., Harvey, D., Li, R., Massey, R., Negrello, M., \& Robertson, A. (2021). PyAutoLens: Open-Source Strong Gravitational Lensing.
\textit{The Journal of Open Source Software, 6(58), 2825}. \href{https://doi.org/10.21105/joss.02825}{https://doi.org/10.21105/joss.02825}}

\hypertarget{Oguri:2010}{
Oguri, M. (2010). glafic: Software Package for Analyzing Gravitational Lensing (ascl:1010.012).}

\hypertarget{OM10}{
Oguri, M., \& Marshall, P. J. (2010). Gravitationally lensed quasars and supernovae in future wide-field optical imaging surveys.
\textit{Monthly Notices of the Royal Astronomical Society,
405(4), 2579–2593}. \href{https://doi.org/10.1111/j.1365-2966.2010.16639.x}{https://doi.org/10.1111/j.1365-2966.2010.16639.x}}

\hypertarget{Pagano:2020}{
Pagano, G., Hannuksela, O. A., \& Li, T. G. F. (2020). LENSINGGW: a PYTHON package for lensing of gravitational waves.
\textit{Astronomy \& Astrophysics, 643, A167}.
\href{https://doi.org/10.1051/0004-6361/202038730}{https://doi.org/10.1051/0004-6361/202038730}}

\hypertarget{Park:2021}{
Park, J. W., Wagner-Carena, S., Birrer, S., Marshall, P. J., Lin, J. Y.-Y., Roodman, A., \&
LSST Dark Energy Science Collaboration. (2021). Large-scale Gravitational Lens Modeling with Bayesian Neural Networks for Accurate and Precise Inference of the Hubble Constant.
\textit{The Astrophysical Journal, 910(1), 39}. \href{https://doi.org/10.3847/1538-4357/abdfc4}{https://doi.org/10.3847/1538-4357/abdfc4}}

\hypertarget{PixeLens}{
Saha, P., \& Williams, L. L. R. (2011). PixeLens: A Portable Modeler of Lensed Quasars (ascl:1102.007).}

\hypertarget{Shajib:2019unified}{
Shajib, A. J. (2019). Unified lensing and kinematic analysis for any elliptical mass profile.
\textit{Monthly Notices of the Royal Astronomical Society, 488(1), 1387–1400}.
\href{https://doi.org/10.1093/mnras/stz1796}{https://doi.org/10.1093/mnras/stz1796}}

\hypertarget{Shajib:2020strides}{
Shajib, A. J., Birrer, S., Treu, T., Agnello, A., Buckley-Geer, E. J., Chan, J. H. H., Christensen,
L., Lemon, C., Lin, H., Millon, M., Poh, J., Rusu, C. E., Sluse, D., Spiniello, C., Chen,
G. C.-F., Collett, T., Courbin, F., Fassnacht, C. D., Frieman, J., … Zhang, Y. (2020).
STRIDES: a 3.9 per cent measurement of the Hubble constant from the strong lens system DES J0408-5354.
\textit{Monthly Notices of the Royal Astronomical Society, 494(4),
6072–6102}. \href{https://doi.org/10.1093/mnras/staa828}{https://doi.org/10.1093/mnras/staa828}}

\hypertarget{Shajib:2019strides}{
Shajib, A. J., Birrer, S., Treu, T., Auger, M. W., Agnello, A., Anguita, T., Buckley-Geer, E.
J., Chan, J. H. H., Collett, T. E., Courbin, F., Fassnacht, C. D., Frieman, J., Kayo, I.,
Lemon, C., Lin, H., Marshall, P. J., McMahon, R., More, A., Morgan, N. D., … Walker,
A. R. (2019). Is every strong lens model unhappy in its own way? Uniform modelling of a sample of 13 quadruply+ imaged quasars.
\textit{Monthly Notices of the Royal Astronomical
Society, 483(4), 5649–5671}.
\href{https://doi.org/10.1093/mnras/sty3397}{https://doi.org/10.1093/mnras/sty3397}}

\hypertarget{Shajib:2021AO}{
Shajib, A. J., Molina, E., Agnello, A., Williams, P. R., Birrer, S., Treu, T., Fassnacht, C.
D., Morishita, T., Abramson, L., Schechter, P. L., \& Wisotzki, L. (2021). High-resolution imaging follow-up of doubly imaged quasars.
\textit{Monthly Notices of the Royal Astronomical
Society, 503(2), 1557–1567}. \href{https://doi.org/10.1093/mnras/stab532}{https://doi.org/10.1093/mnras/stab532}}

\hypertarget{Shajib:2021slacs}{
Shajib, A. J., Treu, T., Birrer, S., \& Sonnenfeld, A. (2021). Dark matter haloes of massive elliptical galaxies at z $\sim$ 0.2 are well described by the Navarro-Frenk-White profile.
\textit{Monthly Notices of the Royal Astronomical Society, 503(2), 2380–2405}. \href{https://doi.org/10.1093/mnras/stab536}{https://doi.org/10.1093/mnras/stab536}}

\hypertarget{Sonnenfeld:2015}{
Sonnenfeld, A., Treu, T., Marshall, P. J., Suyu, S. H., Gavazzi, R., Auger, M. W., \& Nipoti,
C. (2015). The SL2S Galaxy-scale Lens Sample. V. Dark Matter Halos and Stellar IMF of Massive Early-type Galaxies Out to Redshift 0.8.
\textit{The Astrophysical Journal, 800(2),
94}. \href{https://doi.org/10.1088/0004-637X/800/2/94}{https://doi.org/10.1088/0004-637X/800/2/94}}

\hypertarget{Roman}{
Spergel, D., Gehrels, N., Breckinridge, J., Donahue, M., Dressler, A., Gaudi, B. S., Greene,
T., Guyon, O., Hirata, C., Kalirai, J., Kasdin, N. J., Moos, W., Perlmutter, S., Postman,
M., Rauscher, B., Rhodes, J., Wang, Y., Weinberg, D., Centrella, J., … Shaklan, S. (2013).
Wide-Field InfraRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRSTAFTA
Final Report.
\textit{arXiv e-Prints}, arXiv:1305.5422.
\href{http://arxiv.org/abs/1305.5422}{http://arxiv.org/abs/1305.5422}}

\hypertarget{Spilker:2016}{
Spilker, J. S., Marrone, D. P., Aravena, M., Béthermin, M., Bothwell, M. S., Carlstrom, J.
E., Chapman, S. C., Crawford, T. M., de Breuck, C., Fassnacht, C. D., Gonzalez, A. H.,
Greve, T. R., Hezaveh, Y., Litke, K., Ma, J., Malkan, M., Rotermund, K. M., Strandet,
M., Vieira, J. D., … Welikala, N. (2016). ALMA Imaging and Gravitational Lens Models of South Pole TelescopeSelected Dusty, Star-Forming Galaxies at High Redshifts.
\textit{The Astrophysical Journal, 826, 112}. \href{https://doi.org/10.3847/0004-637X/826/2/112}{https://doi.org/10.3847/0004-637X/826/2/112}}

\hypertarget{Tessore:2015}{
Tessore, N., \& Metcalf, R. B. (2015). The elliptical power law profile lens.
\textit{Astronomy \& Astrophysics, 580, A79}. \href{https://doi.org/10.1051/0004-6361/201526773}{https://doi.org/10.1051/0004-6361/201526773}}

\hypertarget{vdVyvere:2020}{
Van de Vyvere, L., Sluse, D., Mukherjee, S., Xu, D., \& Birrer, S. (2020). The impact of mass
map truncation on strong lensing simulations.
\textit{Astronomy \& Astrophysics, 644, A108}.
\href{https://doi.org/10.1051/0004-6361/202038942}{https://doi.org/10.1051/0004-6361/202038942}}

\hypertarget{Vegetti:2012}{
Vegetti, S., Lagattuta, D. J., McKean, J. P., Auger, M. W., Fassnacht, C. D., \& Koopmans,
L. V. E. (2012). Gravitational detection of a low-mass dark satellite galaxy at cosmological distance.
\textit{Nature, 481(7381), 341–343}.
\href{https://doi.org/10.1038/nature10669}{https://doi.org/10.1038/nature10669}}

\hypertarget{Wagner-Carena:2021}{
Wagner-Carena, S., Park, J. W., Birrer, S., Marshall, P. J., Roodman, A., Wechsler, R. H., \& LSST Dark Energy Science Collaboration. (2021). Hierarchical Inference with Bayesian
Neural Networks: An Application to Strong Gravitational Lensing.
\textit{The Astrophysical Journal, 909(2), 187}. \href{https://doi.org/10.3847/1538-4357/abdf59}{https://doi.org/10.3847/1538-4357/abdf59}}

\hypertarget{Wong:2020}{
Wong, K. C., Suyu, S. H., Chen, G. C.-F., Rusu, C. E., Millon, M., Sluse, D., Bonvin, V.,
Fassnacht, C. D., Taubenberger, S., Auger, M. W., Birrer, S., Chan, J. H. H., Courbin, F.,
Hilbert, S., Tihhonova, O., Treu, T., Agnello, A., Ding, X., Jee, I., … Meylan, G. (2020).
H0LiCOW XIII. A 2.4 per cent measurement of H0 from lensed quasars: 5.3$\sigma$ tension between early- and late-Universe probes.
\textit{Monthly Notices of the Royal Astronomical
Society, 498(1), 1420–1439}. \href{https://doi.org/10.1093/mnras/stz3094}{https://doi.org/10.1093/mnras/stz3094}}

\hypertarget{Yang:2020}{
Yang, L., Birrer, S., \& Treu, T. (2020). A versatile tool for cluster lensing source reconstruction - I. Methodology and illustration on sources in the Hubble Frontier Field Cluster MACS
J0717.5+3745.
\textit{Monthly Notices of the Royal Astronomical Society, 496(3), 2648–2662}.
\href{https://doi.org/10.1093/mnras/staa1649}{https://doi.org/10.1093/mnras/staa1649}}

\hypertarget{Yang:2021}{
Yang, L., Roberts-Borsani, G., Treu, T., Birrer, S., Morishita, T., \& Bradač, M. (2021). The
evolution of the size-mass relation at z = 1-3 derived from the complete Hubble Frontier Fields data set.
\textit{Monthly Notices of the Royal Astronomical Society, 501(1), 1028–1037}.
\href{https://doi.org/10.1093/mnras/staa3713}{https://doi.org/10.1093/mnras/staa3713}}

\hypertarget{Zheng:2012}{
Zheng, W., Postman, M., Zitrin, A., Moustakas, J., Shu, X., Jouvel, S., H{\o}st, O., Molino,
A., Bradley, L., Coe, D., Moustakas, L. A., Carrasco, M., Ford, H., Ben{\'\i}tez, N., Lauer,
T. R., Seitz, S., Bouwens, R., Koekemoer, A., Medezinski, E., … van der Wel, A. (2012).
A magnified young galaxy from about 500 million years after the Big Bang. \textit{Nature, 489(7416), 406–408}. \href{https://doi.org/10.1038/nature11446}{https://doi.org/10.1038/nature11446}}

\end{document}