# A Quasi-Newton method for physically-admissible simulation of Poiseuille
flow under fracture propagation.

Guotong Ren
Department of Petroleum Engineering
University of Tulsa
Tulsa, OK 740104
<EMAIL_ADDRESS>
&Rami M. Younis
Department of Petroleum Engineering
University of Tulsa
Tulsa, OK 74104
<EMAIL_ADDRESS>

###### Abstract

Coupled hydro-mechanical processes are of great importance to numerous
engineering systems, e.g., hydraulic fracturing, geothermal energy, and carbon
sequestration. Fluid flow in fractures is modeled after a Poiseuille law that
relates the conductivity to the aperture by a cubic relation. Newton’s method
is commonly employed to solve the resulting discrete, nonlinear algebraic
systems. It is demonstrated, however, that Newton’s method will likely
converge to nonphysical numerical solutions, resulting in estimates with a
negative fracture aperture. A Quasi-Newton approach is developed to ensure
global convergence to the physical solution. A fixed-point stability analysis
demonstrates that both physical and nonphysical solutions are stable for
Newton’s method, whereas only physical solutions are stable for the proposed
Quasi-Newton method. Additionally, it is also demonstrated that the Quasi-
Newton method offers a contraction mapping along the iteration path. Numerical
examples of fluid-driven fracture propagation demonstrate that the proposed
solution method results in robust and computationally efficient performance.

_K_ eywords Quasi-Newton, Fracture propagation, Coupled hydro-mechanics,
Extended finite element method, Finite volume method.

##  1 Introduction

In hydro-mechanical processes, the interplay between rock deformation and
fluid pressure in fractures can dictate first-order effects in several
engineering systems such as hydraulic fracturing [2, 14], geothermal
utilization [3], and $\text{CO}_{2}$ sequestration [5, 18]. In the limits of
linearity, and assuming a constant positive-definite permeability tensor in
fracture, existence and uniqueness of solution to the governing continuity
equations has been established [7, 6]. The physically-accepted limit of
Poiseuille flow introduces a nonlinear relation for fracture conductivity as a
cubic function of local aperture. While formal conditions for the existence
and uniqueness of solutions are well-established, numerical and semi-
analytical approximations are widely applied. Newton’s method is commonly
employed to solve resulting algebraic equations whether it be in the contexts
of fully- or sequentially-coupled solution processes. Numerical convergence
difficulties can arise, as is generally the case with black-box nonlinear
solution processes. More notably however, it is observed that the Newton
iteration can converge to nonphysical solutions (i.e. aperture is negative).
This is a frequently occurring issue and numerical evidence is easily
obtained. Consider for example, a simple two-dimensional model with a fracture
located at the center of a domain, and where fluid is injected into the mid-
point of the fracture. Initially, the scaled fluid pressure and displacement
fields are zero. Two alternate aperture profiles along the fracture are
presented in fig. 1(b); while both solutions satisfy the same nonlinear
convergence criteria for the same model, one profile is oscillatory and
includes negative aperture, whereas the other does not. Negative fracture
(aperture) conductivity violates the positive-definite tensor requirement of
[7, 6], and permits flow from low fluid pressure to high fluid pressure. This
simple example may support the hypothesis that the hydro-mechanical
formulation possesses multiple solutions and that the nonphysical one can be
obtained by application of Newton’s method. We omit the specifics of the
numerical approximation and physical parameters employed in this illustrative
example; in a subsequent section, it is demonstrated that the nonphysical
solution exists regardless of the discretization employed.

(a)

(b)

Figure 1: Example of multiple solutions: (a) case schematic; (b) aperture
profile. The results are generated using coupled XFEM-FVM [15]

While significant to robust numerical simulation, the literature on the
aforementioned issues is scarce. In this article, our focus is on linear
elastic fracture propagation problems without fluid lag; i.e., fluid is
assumed to occupy the entire fracture space at each instance in which the
fracture tip reaches a failure state. Two classes of numerical fracture
propagation simulation algorithm are: (1) Given a fixed tip advancement length
step, determine the time increment such that tip is in an equilibrium failure
state, and the solution satisfies the continuity equations (e.g., [9, 13,
15]); and (2) given a time step, determine the advancement length step for
each tip such that the solution fulfills equilibrium failure (e.g., [8, 10]).
In both classes of algorithm, the nonphysical solution may be encountered. [9,
15] propose initialization schemes for unknowns at every time that new
fracture segments are introduced based on surrogate models. Nevertheless,
while these initial guesses may improve convergence, it is observed that
convergence to a nonphysical solution remains likely to occur in viscosity-
dominated fracture propagation. [7] propose a safeguarding strategy by
applying a relaxation to $\omega$ within the iterative coupled scheme. As
$\omega<0$ is detected during nonlinear iterations, the $\omega$ iterate is
manually reassigned to a non-negative value. Although such an approach ensures
a positive fracture conductivity should the method converge, an open question
pertains to the effects on convergence rate and reliability.

A Quasi-Newton method is proposed to address global convergence to positive
aperture solutions. The method proposes a Jacobian matrix modification that is
inspired by that used in nonlinear two-point-flux finite volume methods (e.g.
[12, 17]). For example in [12, 17], in order to maintain a monotone Jacobian
matrix, derivatives of the transmissibility with respect to pressure are
neglected. As a result, the converged pressure solutions are assured to be
non-negative. For the fracture propagation problem and coupled hydro-
mechanical problems at hand, monotonicity of the linearized operator is not
guaranteed. For instance, the physical pressure solution may take on negative
values near the tip region in the viscosity-dominated fracture propagation
regime (fluid lag zones for example). Subsequently, the proposed approach will
demonstrate that: (1) the fixed-point stability for the Newton’s and proposed
Quasi-Newton methods are different; (2) the proposed Quasi-Newton method
provides a contraction mapping property with Lipschitz constant between $0$
and $1$ on a designed iteration path; and (3) a non-negative aperture
$\omega>=0$ is ensured at every nonlinear iteration.

Next, the mathematical formulation of the problem is detailed in section 2.
Then, in section 3, the Quasi-Newton method is proposed and analyzed under the
setting of an infinite domain where the mechanics equation can be formulated
using a closed Green’s function between $\omega$ and pressure $p$. In section
3.2, the development is extended to more general fracture configurations.
Specifically, we extend the algorithm to discretization using a coupled
extended-finite-element method (XFEM) and finite volume method (FVM) for
mechanics and flow, respectively [15].

##  2 Problem statement

Consider a spatial domain $\Omega\subset\mathbb{R}^{2}$ with external boundary
$\Gamma$ and its associated outward-oriented unit-normal $\bm{n}_{\Gamma}$.
Dirichlet and Neumann boundaries for mechanics are $\Gamma_{u}$ and
$\Gamma_{t}$. The boundary segments are disjoint
($\Gamma_{u}\bigcap\Gamma_{t}=\varnothing$), and
$\Gamma_{u}\cup\Gamma_{t}=\Gamma$.

A one-dimensional lower representation of the fracture, $\mathcal{C}$, is
adopted. Fracture aperture, $\omega$, is defined as,

$\omega=\llbracket\bm{u}\rrbracket\cdot\bm{n}_{c}.$ (1)

where $\llbracket\bm{u}\rrbracket$ is the jump of the displacement over the
fracture, $\bm{n}_{c}$ is the fracture unit-normal vector.

A single-phase incompressible fluid occupies fracture, and linear elastic
mechanical deformation model is assumed. The continuity equation in the
fracture is

$\frac{\partial\omega}{\partial
t}-\nabla_{c}\cdot\left(\frac{\omega^{3}}{12\mu}\nabla_{c}p\right)=0\quad\text{on}\quad\mathcal{C},$
(2)

where $\nabla_{c}$ is the gradient operator defined on the fracture path. $p$
is the fluid pressure in the fracture and $\mu$ is the fluid viscosity.
Poiseuille’s law posits that the fracture conductivity is
$\frac{\omega^{3}}{12}$. The boundary condition at the inlet of the fracture
is

$\displaystyle-\frac{\omega^{3}}{12\mu}\frac{\partial p}{\partial x}=Q,$ (3)

where $Q$ is the fluid injection rate, and at the tip

$\displaystyle-\frac{\omega^{3}}{12\mu}\frac{\partial p}{\partial x}=0.$ (4)

The quasi-static momentum equation reads

$\nabla\cdot\bm{\sigma}=\bm{0},\quad\text{on}\quad\Omega,$ (5)

where $\bm{\sigma}$ is the second order tensor. On the external boundary
$\Gamma$, Neumann (force) and Dirichlet (displacement) conditions are

$\displaystyle\bm{\sigma}\cdot\bm{n}_{\Gamma}=\bm{t}\quad$
$\displaystyle\text{on}\quad\Gamma_{u},$ (6a)
$\displaystyle\bm{u}=\bm{\hat{u}}\quad$
$\displaystyle\text{on}\quad\Gamma_{t},$ (6b)

while on immersed fracture boundaries, $p$ is imposed onto the oriented
surfaces of the fracture:

$\displaystyle\bm{\sigma}\cdot\bm{n_{c}}=p\bm{I}\cdot\bm{n_{c}}\quad$
$\displaystyle\text{on}\quad\mathcal{C}.$ (7a)

The stress $\bm{\sigma}$ is modeled using linear elastic theory:

$\bm{\bm{\sigma}}=\frac{E\nu}{(1+\nu)(1-2\nu)}(\nabla\cdot\bm{u})\bm{I}+\frac{E}{1+\nu}\bm{\varepsilon},$
(8)

where $\bm{I}$ is the identity matrix, $E$ is Young’s modulus, $\nu$ is
Poisson’s ratio and the strain $\bm{\varepsilon}$ is a second order tensor.
Under infinitesimal deformation, the strain tensor is a function of
displacement as

$\bm{\varepsilon}=\frac{1}{2}(\nabla^{\text{T}}\bm{u}+\nabla\bm{u}).$ (9)

###  2.1 Numerical Discretization

Two discretization schemes will be employed in computational examples which
are referred to as DS1 and DS2. DS1 considers a single fracture in an infinite
domain while DS2 deals with rather general fracture and domain geometry. The
simulation test performed only considers the static fracture. The propagation
scenario will be investigated in the next section. Both schemes are briefly
introduced.

####  2.1.1 DS1

A single fracture is modeled in an infinite domain as shown in fig. 1(a). Due
to symmetry, only half of the domain is modeled. According to [16], the
aperture field $\omega$ is obtained by solving eqs. 5, 7, 8 and 9, and can be
written as an explicit function of $p$ by the Green’s function $\mathcal{G}$,

$\begin{split}\omega(x)&=-\frac{2(1-\nu^{2})}{\pi
E}\int_{0}^{a}\mathcal{G}(s;x)p(s)ds\\\ &=-\frac{2(1-\nu^{2})}{\pi
E}\int_{0}^{a}\ln{\left|\frac{(a^{2}-x^{2})^{\frac{1}{2}}-(a^{2}-s^{2})^{\frac{1}{2}}}{(a^{2}-x^{2})^{\frac{1}{2}}+(a^{2}-s^{2})^{\frac{1}{2}}}\right|}p(s)ds\end{split}$
(10)

where $a$ is the fracture half length. Consequently, an integro-differential
equation system, eqs. 2 and 10, is formulated. Next, a discrete approximation
of this formulation is described.

The domain $\mathcal{C}$ is approximated with $n_{c}$ equally-spaced grid
cells using length equal to $\Delta x=\frac{a}{n_{c}}$.
$\bm{\omega}=\\{\omega_{i},i=1,...,n_{c}\\}$ and
$\bm{p}=\\{p_{i},i=1,...,n_{c}\\}$, both of which are co-located at the center
of $n_{c}$ cells. A Gaussian quadrature integration rule is used to evaluate
$\omega_{i}$,

$\begin{split}\omega_{i}\big{(}x_{i}=\frac{2i-1}{2}\Delta
x\big{)}&=-\frac{2(1-\nu^{2})}{\pi E}\sum_{j=1}^{n_{c}}\int_{(j-1)\Delta
x}^{j\Delta x}\mathcal{G}(s;x_{i})p_{j}ds\\\ &=-\frac{2(1-\nu^{2})}{\pi
E}\sum_{j=1}^{n_{c}}\sum_{k=1}^{n_{g}}\mathcal{G}(s^{*}_{k};x_{i})\mathcal{W}_{k}p_{j}\end{split}$
(11)

where $n_{g}$ is the number of quadrature points; $s_{k}^{*}$ is the
quadrature point; $\mathcal{W}_{k}$ is the weight at the point $s_{j}^{*}$.
eq. 11 can be simplified using a matrix-vector notation,

$\bm{\omega}=\bm{A}\bm{p}$ (12)

where each item in full matrix $\bm{A}$ stores the numerical integration
result of each cell. Since for all $s\text{ and }x\subseteq(0,a)$,
$\mathcal{G}(s;x)<0$, then  $A_{ij}>0$.

Pertaining to fluid flow, a first-order backward Euler temporal
discretization, and a second-order central-difference spatial discretization
are applied. The discretized form reads,

$(\omega^{n+1}_{i}-\omega_{i}^{n})-T\sum_{\pm}\left(\omega^{n+1}_{i\pm\frac{1}{2}}\right)^{3}(p^{n+1}_{i\pm
1}-p^{n+1}_{i})=0$ (13)

where superscripts $n$ or $n+1$ indicate time level, and $\Delta
t=t^{n+1}-t^{n}$; $T=\frac{\Delta t}{12\mu\Delta x^{2}}$ is the static
transmissibility while $\omega$ at the interface is evaluated by the
arithmetic average of its neighboring cells’, i.e.
$\omega_{i\pm\frac{1}{2}}=\frac{\omega_{i}+\omega_{i\pm 1}}{2}$. Combining
eqs. 12 and 13, a matrix-vector form of the nonlinear system is formulated

$\big{(}\bm{A}+\bm{F}(\bm{p}^{n+1})\big{)}\bm{p}^{n+1}=\bm{q}+\bm{\omega}^{n}$
(14)

where $\bm{q}=(\frac{\Delta t}{\Delta x}Q,0,...,0)^{T}$ is a $n_{c}$ by $1$
vector, $\bm{F}$ represents the flux term and is a tri-diagonal sparse matrix
under a 2D setting. The detailed expression of $\bm{F}$ is given in A.

####  2.1.2 DS2

In general scenarios with multiple fractures, as well as complex-fracture
geometries, the original problem is approximated using a mixed discretization;
coupled extended FEM and embedded FVM for mechanics deformation and fluid
flow, respectively. The numerical details of the method is referred to the
appendix C.

##  3 Proposed Quasi-Newton Solver

Newton’s method is typically applied to solve nonlinear algebraic systems that
arise from discretization. It is observed that these methods may converge to
nonphysical solutions. In the following subsections, this is analyzed
theoretically and empirically in the contexts of models DS1 and DS2 introduced
above. Furthermore, a Quasi-Newton approach is proposed, and demonstrated to
provide robust global convergence to the physical solution.

###  3.1 Analysis and development using model DS1

The nonlinear system eq. 14 is to be solved, where the residual vector
$\bm{R}$ and Jacobian matrix $\bm{J}$ in the context of DS1 become,

$\displaystyle
R^{n+1,v+1}_{i}=(\omega^{n+1,v}_{i}-\omega_{i}^{n})-T\sum_{\pm}\left(\omega^{n+1,v}_{i\pm\frac{1}{2}}\right)^{3}(p^{n+1,v}_{i\pm
1}-p^{n+1,v}_{i})-q_{i}$ (15a)
$\displaystyle\begin{split}&J^{n+1,v+1}_{ij}=\frac{\partial\omega^{n+1,v}_{i}}{\partial
p^{n+1,v}_{j}}-T\sum_{\pm}\left(\omega^{n+1,v}_{i\pm\frac{1}{2}}\right)^{3}\frac{\partial(p^{n+1,v}_{i\pm
1}-p^{n+1,v}_{i})}{\partial p_{j}^{n+1,v}}\\\
&-T\sum_{\pm}\frac{\partial\left(\omega^{n+1,v}_{i\pm\frac{1}{2}}\right)^{3}}{\partial
p^{n+1,v}_{j}}(p^{n+1,v}_{i\pm 1}-p^{n+1,v}_{i})\end{split}$ (15b)
$\displaystyle\bm{J}^{n+1,v+1}\delta\bm{p}^{n+1,v+1}=-\bm{R}^{n+1,v+1}$ (15c)

where superscripts $v$ or $v+1$ represent the iteration level. The first term
on the right-hand side (RHS) of eq. 15b,
$\frac{\partial\omega^{n+1,v}_{i}}{\partial p^{n+1,v}_{j}}$, is equal to
$A_{ij}$; the second term on the RHS of eq. 15b forms a diagonally dominant
sparse matrix whose diagonal terms are positive and off-diagonal terms are
negative. The last term on the RHS of eq. 15b produces a full matrix with
elements consisting of derivatives of $\omega^{3}$ with respect to $p$. The
subsequent numerical study demonstrates that the presence of this term can
lead to convergence of Newton’s method to a nonphysical solution. In the
proposed Quasi-Newton approach, the third term is neglected in the Jacobian
calculation. Note that in [11, 17], a similar strategy is adopted but for the
purpose of guaranteeing a monotone linear matrix. The QN method now reads,

$\displaystyle\tilde{J}^{n+1,v+1}_{ij}=\frac{\partial\omega^{n+1,v}_{i}}{\partial
p^{n+1,v}_{j}}-T\sum_{\pm}\left(\omega^{n+1,v}_{i\pm\frac{1}{2}}\right)^{3}\frac{\partial(p^{n+1,v}_{i\pm
1}-p^{n+1,v}_{i})}{\partial p_{j}^{n+1,v}}$ (16a)
$\displaystyle\tilde{\bm{J}}^{n+1,v+1}\bm{p}^{n+1,v+1}=\big{(}\bm{A}+\bm{F}(\bm{p}^{n+1,v})\big{)}\bm{p}^{n+1,v+1}=\bm{q}+\bm{\omega}^{n}$
(16b)

From eq. 16b, a mapping $K:X\rightarrow X,X\subset\mathbb{R}^{n_{c}}$ is
defined as

$K=\big{(}\bm{A}+\bm{F}(\bm{p})\big{)}^{-1}(\bm{q}+\bm{\omega}^{n})$ (17)

Consequently, our Quasi-Newton eq. 16b searches for fixed points over set $X$.
The following proposition proposes a property on set $X$,

######  Proposition 1.

$\forall\bm{p}\in X$, mapping $K:X\rightarrow X$, and
$X\subset\mathbb{R}^{n_{c}}$ is a vector space that satisfies
$\sum_{i}(\bm{A}\bm{p}^{*})_{i}=\sum_{i}(\bm{q}+\bm{\omega}^{n})_{i},\text{
where }\bm{p}^{*}=K(\bm{p})$.

Proposition 1 states that the total mass balance is satisfied during each
iteration. The proof of the Proposition 1 is in Appendix B.

####  3.1.1 Analysis of fixed-point stability

There may exist multiple fixed-points in $X$. The stability of such fixed-
points for nonlinear mappings $K$ corresponding to both Newton’s and Quasi-
Newton methods is analyzed. The stability of fixed points is defined by
Definition 1.

######  Definition 1.

Fixed-points can be classified as,

  * •

The fixed point $\bm{p}_{0}$ is stable if there exists an open set $U\subset
X$ containing $\bm{p}_{0}$ such that
$\|K(\bm{p})-\bm{p}_{0}\|\leq\|\bm{p}-\bm{p}_{0}\|$ for all $\bm{p}\in U$

  * •

The fixed point $\bm{p}_{0}$ is unstable if there exists an open set $U\subset
X$ containing $\bm{p}_{0}$ such that
$\|K(\bm{p})-\bm{p}_{0}\|\geq\|\bm{p}-\bm{p}_{0}\|$ for all $\bm{p}\in U$

A well-studied result in discrete dynamics characterizes the stability of
fixed-points for multi-variable nonlinear operators.

######  Theorem 1.

(see [4]) Suppose map $K:X\rightarrow X$ is differentiable at a fixed point
$\bm{p}_{0}$ and then let $\lambda_{1},\lambda_{2},...,\lambda_{n}$ be the
eigenvalues of the Jacobian matrix $K^{\prime}$ evaluated at $\bm{p}_{0}$.

  * •

if all the eigenvalues $\lambda_{j}$ have $|\lambda_{j}|<1$, then $\bm{p}_{0}$
is stable or attracting.

  * •

if one eigenvalue $\lambda_{j0}$ has $|\lambda_{j0}|>1$, $\bm{p}_{0}$ is
unstable, which can be either saddle or repelling.

  * •

if all the eigenvalues $\lambda_{j}$ have $|\lambda_{j}|>1$, then $\bm{p}_{0}$
is repelling.

An empirical analysis of the stability of the fixed-points is conducted. The
Buckingham $\pi$ theorem is applied in a dimensional analysis to yield the
following dimensionless groups,

$\Pi_{1}=\frac{\mu}{Et},\quad\Pi_{2}=\frac{qt}{a^{2}},\quad\Pi_{3}=\frac{p}{E}$
(18)

Based on the Buckingham $\pi$ theorem, the solution $\Pi_{3}$ only depends on
the magnitude of $\Pi_{1}$ and $\Pi_{2}$. The domain setup is shown in fig.
1(a). A single fracture is in the center of an infinite domain. Initially, the
fracture is void of fluid, i.e. $\bm{\omega}^{n}=0$. The simulation runs for
one time step $\Delta t=t$. The fracture is represented using four grid cells.
The input parameter $\Pi_{1}$ is varied across simulations from $1e-17$ to
$1.5e-1$ and $\Pi_{2}$ from $1e-5$ to $2e-2$. As a result, there are 8000
simulation cases in total. The physical and nonphysical solutions are
generated for each one of 8000 simulation cases. Subsequently, then spectral
radius of $K^{\prime}$, $\rho(K^{\prime})$ is evaluated for each case at the
physical and nonphysical solutions. These results are illustrated in fig. 2.
All nonphysical solutions are unstable fixed-points of the Quasi-Newton
mapping, while all physical solutions are stable. In other words, if the
Quasi-Newton method converges, it will not do so to a nonphysical solution
since any small perturbation will push the next iterate away from the unstable
fixed point. To demonstrate this, we apply the nonphysical solution perturbed
by $0.01\%$ as an initial guess for the Quasi-Newton method. All 8000 cases
converge to the physical solution, with a maximum of 12 nonlinear iterations.
The number of nonlinear iterations required to convergence over the entire
test-set appears in fig. 3.

(a)

(b)

Figure 2: $\log(\rho(K^{\prime}))$ of Quasi-Newton method v.s
$\left(\log(\Pi_{1}),\log(\Pi_{2})\right)$. (a) physical solution; (b)
nonphysical solution.  Figure 3: nonlinear iterations of each case using the
slightly perturbed nonphysical solution as the initial guess

Similarly, the stability of fixed-points for the Newton operator is also
investigated. A mapping $K$ for Newton’s method can be defined as,

$K(\bm{p})=\bm{p}-\bm{J}^{-1}(\bm{p})\bm{R}(\bm{p})$ (19)

The computed spectral radius $\rho(K^{\prime})$ at physical and nonphysical
solutions over the investigation space are presented in fig. 4. Clearly, both
physical and nonphysical solutions are stable fixed points in $X$. This
explains why Newton’s method may converge to the negative aperture solution
when applied to any of these problems.

(a)

(b)

Figure 4: $\log(\rho(K^{\prime}))$ of Newton’s method v.s
$\left(\log(\Pi_{1}),\log(\Pi_{2})\right)$. (a) physical solution; (b)
nonphysical solution.

####  3.1.2 contraction mapping

######  Definition 2.

Let $(D,d)$ define a metric space. Mapping $G:D\rightarrow D$ is a
contraction, if there exists a constant $c$, with $0<c<1$, such that

$d\big{(}G(y_{1})-G(y_{2})\big{)}<cd(y_{1},y_{2})$ (20)

for all $y_{1},y_{2}\in D$. Note that $c$ in eq. 20 is the Lipschitz constant.

Next, it is demonstrated that a mapping $G$ for $\omega$ preserves the
contraction mapping property on a certain iteration path. The RHS of eq. 17 is
left multiplied by matrix $A$, to produce $G$ as,

$G=\bm{A}(\bm{A}+\bm{F}(\bm{\omega}))^{-1}(\bm{q}+\bm{\omega}^{n})$ (21)

Firstly, we define a set $D$ that only contains physical solutions.

######  Definition 3.

Set $D$ is defined such that dimensionless aperture
$\frac{\bm{\omega}}{\sqrt{Qt}}\geq\varepsilon_{0}$. Here $\varepsilon_{0}$ is
a fixed tolerance whose absolute value, $\mid\varepsilon_{0}\mid$, is close to
$0$.

In order to show there exists a $0<c<1$ such that mapping $G:D\rightarrow D$
is a contraction, we design an iteration path for mapping $G$

$G=\bm{A}\big{(}\bm{A}+\bm{F}(\bm{\omega}^{n+1,v})\big{)}^{-1}(\bm{q}+\bm{\omega}^{n})$
(22)

where $v$ stands for applying $G$ for $v$ iterations, and $\omega^{n+1,v}$ is

$\omega^{n+1,v}=\left\\{\begin{matrix}&\omega^{n}&v=1\\\
&\omega^{n+1,v-1}&v\geq 1\end{matrix}\right.$ (23)

where in the first iteration, the previous time-step solution is applied as
the initial guess.

There are two things that need to be demonstrated for a contraction mapping:
(1) if $\frac{\bm{\omega}}{\sqrt{Qt}}\in D$,
$\frac{K(\bm{\omega})}{\sqrt{Qt}}\in D$; (2) there exists $0<c<1$ on the
designed iteration path, eq. 22 and eq. 23. To present the value of $c$ during
the numerical simulation, the following calculation is adopted

$c=\frac{d(G^{v}(\omega^{v}),G^{v}(\omega^{v-1}))}{d(\omega^{v},\omega^{v-1})}=\frac{d(\omega^{v+1},\omega^{v})}{d(\omega^{v},\omega^{v-1})}$
(24)

where operator $d$ is chosen as $L_{2}$ norm in this manuscript.

The same simulation cases used the analysis of fixed-point stability are used
except with $n_{c}=15$. The results appear in fig. 5. In fig. 5(a), the
Lipschitz constant $c$ is below $1$ for all simulation cases. The minimum
dimensionless aperture $\frac{\omega}{\sqrt{qt}}$ is recorded during the
nonlinear iteration. The dark blue area in fig. 5(b) indicates that all
iterations produce an $\omega>0$. In the yellow and orange regions,
$\varepsilon_{0}$ is on the order of $-1e-5$, which is approximately three
orders of magnitude less than the positive values of
$\frac{\omega}{\sqrt{qt}}$. In summary, on the designed iteration path,
$G(\bm{\omega})$ produces aperture $\omega\in D$. fig. 5(c) shows that the
maximum number of iterations needed is below 20.

(a)

(b)

(c)

Figure 5: Quasi-Newton Contraction Verification: (a) the maximum $c$ during
nonlinear iterations; (b) the minimum dimensionless $\frac{\omega}{\sqrt{qt}}$
during nonlinear iterations and $\log(-\frac{\omega}{\sqrt{Qt}})$ is
presented; (c) nonlinear iterations for each case.

######  Remark 1.

Contraction mapping is verified empirically under uniform grid size. In order
to maintain the contraction in non-uniform grid systems, a constraint on the
time-step size may be enforced.

##### Physical interpretations

A simulation case using one of the dimensional groups $(\Pi_{1},\Pi_{2})$ in
the contraction-mapping test is used to motivate the evolution of $\omega$ and
$p$ during the proposed Quasi-Newton iteration. In fig. 6, $\omega$ and $p$
profiles are captured at different nonlinear iteration numbers, $\nu$. A
moving fluid front in the fracture is captured for each iteration. Aperture
profiles are split over two sub-regions: A fluid-filled region where
$\omega>0$ and a fluid-void region where $\omega=0$. As the solver iterates,
$\omega$ is gradually opened by the fluid front. Similarly, on the pressure
profile, the fluid-filled region is represented by positive pressure, whereas
the fluid-void space is indicated by negative pressure.

In terms of the structure of $\bm{F}$, the number of nonzero entries in matrix
$\bm{F}$ increases as the iteration grows. For example, in the first iteration
($\nu=1)$, only the aperture of the first grid cell is positive and,
therefore, the only nonzero entries are first and second rows of $\bm{F}$,
which means there exists flux between the first and second cells. In the
second iteration, since aperture of the first two cells is positive, a third
row will be added into nonzero entries of $\bm{F}$. Consequently, the fluid
front is moving exactly one grid cell after each iteration. As a result, the
number of iterations needed for convergence approximately depends on the
number of fracture cells, as well as the fluid front location at the
convergence. This in turn, may be interpreted to scale with the target time-
size under continued injection.

(a)

(b)

Figure 6: $\omega$ and $p$ profiles during Quasi-Newton’s iteration. Iteration
number, $v$, increases from left to right.

###  3.2 Extension and application to model DS2

The proposed Quasi-Newton formulation is extended to the co-solution of $p$
and $u$ using a general XFEM-FVM discretization. As motivated in the previous
section, the contraction mapping property may be verified computationally in
terms of $\omega$ as an independent-variable. In order to translate to a
similar effect for aperture maintained in the coupled XFEM-FVM system, we
first derive a nonlinear mapping function similar to eq. 21, and multiple
cases are investigated. In the first case, a group of injection tests are
conducted in which the Lipschitz constant $c$ is evaluated and analyzed for
contraction in aperture. In the second case, the viscosity–dominated KGD
problem is studied to show the robustness of the solver. In the final case,
multiple fracture propagation is considered and solver performance is
reported.

In the XFEM-FVM scheme, XFEM offers a relationship between $\bm{\omega}$ and
$\bm{u}$ as,

$\bm{\omega}=\bm{B}\bm{u}$ (25)

where $\bm{B}$ is a rank-deficient linear matrix. Substituting eq. 25 into
eqs. 15a and 15b yields the discretized system as well as Jacobian matrix
using Quasi-Newton method. A similar form, such as eq. 16b, can be written
here for the XFEM–FVM discretized system,

$\begin{bmatrix}\tilde{\bm{J}}_{ff}&\tilde{\bm{J}}_{fm}\\\
\tilde{\bm{J}}_{mf}&\tilde{\bm{J}}_{mm}\end{bmatrix}^{v}\begin{bmatrix}\bm{p}\\\
\bm{u}\end{bmatrix}^{v+1}=\begin{bmatrix}\bm{F}_{f}\\\
\bm{F}_{m}\end{bmatrix}$ (26)

where $f,m$ indicates flow and mechanics, respectively. During iterations,
$\tilde{\bm{J}}_{ff}$ is a function of $\bm{\omega}^{v}$ whose values are
taken from the previous iteration using eq. 25. In the Quasi-Newton approach,
$\tilde{\bm{J}}_{fm}$ is the same as $\bm{B}$ defined in eq. 25. The rest of
the terms in eq. 26 are constant and do not change across iterations.
$\bm{F}_{f}=\bm{q}+\bm{\omega}^{n}$ contains the source or sink terms while
$\bm{F}_{m}$ contains the Neumann boundary information.

In order to derive a nonlinear mapping for $\bm{\omega}$ in the current
system, the Schur complement is constructed to demonstrate that the solution
of $\omega$ by eq. 26 is equivalent to $\bm{\omega}$ solution obtained by
using the following mapping $G$,

$\bm{\omega}=G(\bm{\omega})=\bm{B}\tilde{\bm{J}}^{-1}_{mm}\left[\bm{F}_{m}-\tilde{\bm{J}}_{mf}\bm{S}^{-1}(\bm{F}_{f}-\tilde{\bm{J}}_{fm}\tilde{\bm{J}}^{-1}_{mm}\bm{F}_{m})\right]$
(27)

where all of terms in eq. 27 are constant except that the Schur complement
$\bm{S}$ is a function of $\bm{\omega}$, which is defined as

$\bm{S}=\tilde{\bm{J}}_{ff}-\tilde{\bm{J}}_{fm}\tilde{\bm{J}}^{-1}_{mm}\tilde{\bm{J}}_{mf}$
(28)

Starting from eq. 26, the deformation solution,

$\bm{u}=\tilde{\bm{J}}^{-1}_{mm}(\bm{F}_{m}-\tilde{\bm{J}}_{mf}\bm{p})$ (29)

and the pore-pressure solution is,

$\bm{p}=\bm{S}^{-1}(\bm{F}_{f}-\tilde{\bm{J}}_{fm}\tilde{\bm{J}}^{-1}_{mm}\bm{F}_{m})$
(30)

Note that substituting eq. 30 into eq. 29 and left multiplication of the RHS
by $\bm{B}$ results in eq. 27.

####  3.2.1 Injection test

A schematic of the case is illustrated in fig. 1(a). The length of the domain
is $100\text{ m}$ by $100\text{ m}$. The fracture half-length is taken as
$40\text{ m}$. The displacement at the mid point of four edges is fixed and
stresses on the boundaries are assumed to be zero. Initially, the fracture is
void of fluid, i.e., $\omega^{0}=0\text{ m}$. An injection rate of
$0.001\text{ m}^{2}/\text{s}$ is applied and a time step $\Delta t$ is set at
$85\text{ s}$. The Young’s modulus and Poisson’s ratio for the rock are
$8.3\text{ GPa}$ and $0.25$, respectively. Plane strain conditions are
assumed.

Fluid is injected at the mid-point of the fracture for a single-time step
$\Delta t$. The mesh sizes and fluid viscosity are varied to verify the
designed iteration path offers a contraction mapping using XFEM–FVM, where $c$
with respect to $\omega$ is calculated, and the result is shown in table 1.
First, $c$ does not exceed $1$ in any test case. Secondly, as $\mu$ is
increased or mesh size is decreased, $c$ decreases. These observations agree
with the results in fig. 5(a).

Next, time is marched forward with a fixed $\Delta t=85\text{ sec}$ until the
injected fluid fills the entire fracture. Four different mesh sizes in table 1
with $\mu=20\text{ Pa}\cdot\text{s}$ are evaluated. The maximum $c$ in each
time step solution and snapshots of $\omega$ profiles are shown in fig. 7. In
fig. 7(b), during injection, $\omega$ is gradually opened while the fluid
front is moving from the injection point to the fracture tip. Before the
injected fluid fills the fracture, $c$ remains below 1 under the four-tested
meshes presented in fig. 7(a). A higher value of $c$ is a result of
subsequently opened fracture segments with each nonlinear iteration. Once the
fluid front reaches the fracture tip, $c$ drops significantly below $0.1$. A
significant drop of $c$ indicates an improvement of local convergence speed.

Table 1: Maximum $c$ during nonlinear iterations (first number in the bracket)
and the number of nonlinear iterations (second number in the bracket)

mesh sizes $\mu(\text{Pa}\cdot\text{s})$ 20 200 2000 20000 $157\times 107$
$(0.83,34)$ $(0.81,25)$ $(0.69,25)$ $(0.53,17)$ $257\times 107$ $(0.89,44)$
$(0.86,39)$ $(0.82,26)$ $(0.75,26)$ $357\times 107$ $(0.92,51)$ $(0.89,40)$
$(0.87,42)$ $(0.76,35)$ $405\times 107$ $(0.95,65)$ $(0.91,59)$ $(0.88,37)$
$(0.83,41)$

(a)

(b)

Figure 7: (a) $\text{max}(c)$ of the case $\mu=20\text{ Pa}\cdot\text{s}$
recorded at each time step solve; (b) snapshots of $\omega$ profiles over half
fracture length at different times. Time increases from right to left. Figure
8: Comparison of three different nonlinear solvers. ”failed” above the ”Newton
with Penalty” bar means the solver fails after a certain number of iterations.

In computational mechanics, the negative aperture is addressed by applying
supporting forces on fracture surfaces, also known as contact force. The
penalty method or Lagrangian multipliers are widely-used numerical techniques.
We briefly describe how the contact force is added into the coupled system in
appendix C. However, a reckless application of contact constraints to hydro-
mechanics system may not overcome the non-physical solutions. We test various
viscosity scenarios under the mesh size $357\times 107$. Two numerical schemes
are considered here: Newton with the penalty method and Quasi-Newton with the
penalty method. The results show in fig. 8. First of all, Newton’s solver with
the penalty method fails for all of test cases here. Newton’s method yields
giant updates which destroy the nonlinear solver system. Since both physical
and nonphysical solutions are stable for Newton’s method, direct application
of the penalty method will not guide Newton’s path to the correct one and give
rise to the failure of the nonlinear solver. On the other hand, all Quasi-
Newton with the penalty method converges to the physical solution at a slight
higher cost than the pure Quasi-Newton method. The Quasi-Newton could possibly
produce very small negative aperture which will introduce the penalty method
into calculation. Nevertheless, it converges to the correct solution at the
end.

####  3.2.2 KGD fracture propagation

The proposed method is applied to the simulation of linear-elastic fluid-
driven fracture propagation. KGD fracture is considered under two propagation
regimes; viscosity-dominated and toughness-dominated. These are are controlled
by the dimensionless parameter $\mathcal{K}_{m}$ ([19])

$\mathcal{K}_{m}=\frac{8K_{c}}{(2E^{{}^{\prime}3}\mu^{{}^{\prime}}q)^{1/4}},$
(31)

where $K_{c}$ is the critical stress intensity factor (SIF);
$\mu^{{}^{\prime}}=12\mu$, $E^{\prime}=E/(1-\nu^{2})$ is the equivalent
Young’s modulus in the plain strain condition; and $q$ is the flow rate into
two wings of the fracture. If $\mathcal{K}_{m}<1$, the flow lies in the
viscosity-dominated regime where rock is very brittle ($K_{c}\rightarrow 0$)
and energy dissipation is dominant in viscous flow. On the other hand, if
$\mathcal{K}_{m}>4$, the process becomes toughness-dominated regime where
energy is mostly used to break the rock and factors from either small aperture
or highly-viscous fluid could be neglected.

The rock-failure criterion is based on Irwin’s law of linear elastic fracture
mechanics and SIF, $K_{eq}$, is used as an indicator. $K_{eq}$ is never
allowed to be larger than the critical $K_{c}$ during fracture propagation.
When $K_{eq}<K_{c}$, fracture is static. Otherwise, fracture is extended by a
certain length $\Delta a$ defined by the user. In terms of fracture-mesh
updates, multiple fracture segments within the same background mesh grid are
prohibited. In this way, multiple small fracture segments will be eliminated
and large grid-size contrast in the numerical calculation can be avoided. For
example, in fig. 9, the fracture segment 3 that partially cuts the grid (fig.
9(a)) will extend to reach the boundary of the grid (fig. 9(b)), when fracture
propagation is triggered. More details of the fracture propagation algorithm
could be referred to appendix D.

(a)

(b)

Figure 9: fracture grids and the numbering; (a) before propagation (b) after
propagation.

The results suggest that Newton’s method is more likely to converge to the
nonphysical solution when fracture propagation lies in the viscosity-dominated
regime. To further investigate this, consider application of the Quasi-Newton
method using the KGD analytical solution in the viscosity-dominated regime.

A square domain of dimensions $100\times 100\text{ m}^{2}$ is modeled with a
Cartesian mesh of size $317\times 117$. The initial fracture position is at
the center of the domain with a half length of $2\text{ m}$. The rock and
fluid parameters are listed in table 2 and the maximum time-step size is
$0.5\text{ sec}$. The total simulation time is $90\text{ sec}$. Two
advancement lengths, $\Delta a=2,4\text{ m}$, are applied. The convergence
criterion is controlled as,

$\sqrt{\frac{\sum_{i=1}^{N}(\omega^{v+1}-\omega^{v})^{2}}{N}}<1e-8$ (32)

where $N$ is the number of fracture grids.

Table 2: Input parameters for KGD model verification  $E$(GPa) | $\nu$ | $K_{c}(\text{MPa}\sqrt{\text{m}})$ | $\mu(\text{Pa}\cdot\text{s})$ | $q(\text{m}^{2}/\text{s})$ | $\mathcal{K}_{m}$
---|---|---|---|---|---
8.3 | 0.25 | 0.5 | 2e-3 | 0.001 | 0.78

(a)

(b)

Figure 10: (a) fracture half length evolution with time (b) $p$ profile when
half fracture length is equal to $16$ m.

The viscosity-dominated analytical solution is generated using the form given
in [1]. A comparison of numerical and analytical results are illustrated in
fig. 10 for both time-step sizes considered. The pressure profile of
viscosity-dominated fracture propagation exhibits negative values near the tip
and asymptotically decays. The singular behavior of pressure at the tip adds
additional difficulty to Newton’s method in obtaining this physical solution.

During the simulation, the Quasi-Newton approach always converges to the
physical solution and the performance is reported in fig. 10. The constant $c$
computed in each iteration over the entire simulation for both step-length
sizes are plotted in fig. 11(a) and fig. 11(b). $c$ is below $1$, confirming
the contraction mapping property on the designed iteration path throughout the
course of entire simulation. The number of iterations for each time step is
illustrated in fig. 11(c) and fig. 11(b).

(a)

(b)

(c)

(d)

Figure 11: (a) and (b) $c$ over the simulation history (c) and (d) the number
of newton iterations in each time step solve. First column corresponds to
results with $\Delta a=2$; Second column corresponds to results with $\Delta
a=4$.

####  3.2.3 multiple fracture propagation

Consider a pre-existing fracture network as illustrated in fig. 12(a). In this
test case, point source fluid injection into the network and fracture
propagation are simulated under various conditions. The input parameters are
listed in table 3. The upper and lower boundaries are under the maximum
horizontal stress $\sigma_{H}$ while the left and right boundaries are under
the minimum one $\sigma_{h}$. The point sources are located on the two
horizontal fractures at $(50,48)$ and $(50,52)$, respectively. The injection
rates that equal $0.001\text{ m}^{2}/\text{s}$ are the same for two point
sources. At $t=0$, fractures are void of fluid. The simulation time lasts
$122$ secs.

Table 3: Input parameters for the case of multiple fracture propagation | value | unit
---|---|---
$E$ | $8.3$ | GPa
$\nu$ | 0.25 |
$\mu$ | 2e-3 | $\text{Pa}\cdot\text{s}$
$K_{c}$ | 2e6 | $\text{Pa}\sqrt{\text{m}}$
$\Delta a$ | 5 | m
$\sigma_{H}$ | 0.4 | MPa
$\sigma_{h}$ | 0.2 | MPa

(a)

(b)

(c)

(d)

Figure 12: (a) initial fracture configuration; (b) (c) and (d) fracture
configurations at three different times in chronological order. Figure 13:
nonlinear iterations for each time step solve Figure 14: fracture geometry
when $\sigma_{H}$ is in the horizontal direction and $\sigma_{h}$ is in the
vertical direction.

The plots in figs. 12(b), 12(c) and 12(d) capture three snapshots during
fracture propagation. At the earliest time, injected fluid first opens two
horizontal fractures while the other three vertical fractures are closed due
to stresses on the boundaries. In figs. 12(b) and 12(c), as three vertical
fractures have the least resistance from far field stresses than the two
horizontal ones, vertical fractures are gradually opened by fluid and start to
propagate in the vertical direction. In the meantime, only one horizontal
fracture starts to propagate. In the latest time shown in fig. 12(d), only
propagation of the vertical fractures is observed in the vertical direction
perpendicular to the minimum principal stress. The number of nonlinear
iterations is shown in fig. 13. There are a few times when the iteration
number is above 30 since nonuniform-spaced fracture grids are inevitably
introduced during fracture propagation . On the other hand, Newton’s method is
also applied to this case and fails in the first time when a fracture starts
to propagate.

We also rotate $\sigma_{H}$ and $\sigma_{h}$ on the boundaries and run the
same simulation case. In this regard, only two horizontal fractures are
advancing in the x direction while three vertical fractures remain static in
fig. 14.

##  4 Discussion and Conclusion

A Quasi-Newton approach is proposed to avoid convergence to nonphysical
solutions and to improve solver robustness for coupled hydro-mechanics and
fracture propagation simulation. Supported by thorough empirical observation
and some analysis, the following findings are listed:

  * •

multiple solutions may exist for coupled hydro-mechanical models of Poiseiulle
flow in fracture, and this has been observed in fully discrete and semi-
analytical models.

  * •

the derivative of the flux function with respect to aperture in the Jacobian
matrix can influence the Newton path towards nonphysical solution.

  * •

empirical stability analysis suggests that the physical solution is the only
stable fixed point for the proposed Quasi-Newton method, whereas both physical
and nonphysical solutions are stable fixed points for Newton’s method.

  * •

the proposed Quasi-Newton method is a contraction for aperture on uniform
grids.

In fracture propagation, the simulation meshing is adjusted with simulation
time. An inappropriate initialization of unknowns in the newly-formed fracture
space can possibly steer the nonlinear Newton solution process for state
(displacement and pressure) to a nonphysical fixed-point attractor. This more
frequently observed in viscosity-dominated propagation problems. The proposed
Quasi-Newton method is demonstrated to overcome these issues in viscosity-
dominated KGD fracture propagation, as well as under the simultaneous
propagation of multiple fractures. The proposed method has been tested under
the setting of a 2D domain, but can be readily extended to 3D domain.

##  5 Acknowledgements

This material is based upon work supported by the U.S. Department of Energy
under Award Number DE-FE-0031777. The authors also acknowledge partial funding
from the members of the TU Future Reservoir Simulation Systems & Technology
(FuRSST) Industry-University Consortium.

## References

  * [1] José Ignacio Adachi.  Fluid-driven fracture in permeable rock.  PhD thesis, 2002.
  * [2] Zdeněk P Bažant, Marco Salviato, Viet T Chau, Hari Viswanathan, and Aleksander Zubelewicz.  Why fracking works.  Journal of Applied Mechanics, 81(10), 2014.
  * [3] Katrin Breede, Khatia Dzebisashvili, Xiaolei Liu, and Gioia Falcone.  A systematic review of enhanced (or engineered) geothermal systems: past, present and future.  Geothermal Energy, 1(1):1–27, 2013.
  * [4] Beatriz Campos, Alicia Cordero, Juan R Torregrosa, and Pura Vindel.  Stability of king’s family of iterative methods with memory.  Journal of Computational and Applied Mathematics, 318:504–514, 2017\.
  * [5] Ryan WJ Edwards, Michael A Celia, Karl W Bandilla, Florian Doster, and Cynthia M Kanno.  A model to estimate carbon dioxide injectivity and storage capacity for geological sequestration in shale gas wells.  Environmental science & technology, 49(15):9222–9229, 2015.
  * [6] Vivette Girault, Kundan Kumar, and Mary F Wheeler.  Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium.  Computational Geosciences, 20(5):997–1011, 2016.
  * [7] Vivette Girault, Mary F Wheeler, Benjamin Ganis, and Mark E Mear.  A lubrication fracture model in a poro-elastic medium.  Mathematical Models and Methods in Applied Sciences, 25(04):587–645, 2015.
  * [8] Elizaveta Gordeliy and Anthony Peirce.  Implicit level set schemes for modeling hydraulic fractures using the xfem.  Computer Methods in Applied Mechanics and Engineering, 266:125–143, 2013.
  * [9] P Gupta and Carlos Armando Duarte.  Coupled hydromechanical-fracture simulations of nonplanar three-dimensional hydraulic fracture propagation.  International Journal for Numerical and Analytical Methods in Geomechanics, 42(1):143–180, 2018.
  * [10] Michael J Hunsweck, Yongxing Shen, and Adrián J Lew.  A finite element approach to the simulation of hydraulic fractures with lag.  International Journal for Numerical and Analytical Methods in Geomechanics, 37(9):993–1015, 2013.
  * [11] Konstantin Lipnikov, Mikhail Shashkov, Daniil Svyatskiy, and Yu Vassilevski.  Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes.  Journal of Computational Physics, 227(1):492–512, 2007.
  * [12] Konstantin Lipnikov, Daniil Svyatskiy, and Yuri Vassilevski.  Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes.  Journal of Computational Physics, 228(3):703–716, 2009.
  * [13] Fushen Liu.  Modeling hydraulic fracture propagation in permeable media with an embedded strong discontinuity approach.  International Journal for Numerical and Analytical Methods in Geomechanics, 44(12):1634–1655, 2020.
  * [14] Zhe Liu and Albert C Reynolds.  History matching an unconventional reservoir with a complex fracture network.  In SPE Reservoir Simulation Conference. OnePetro, 2019.
  * [15] Guotong Ren and Rami M Younis.  An integrated numerical model for coupled poro-hydro-mechanics and fracture propagation using embedded meshes.  Computer Methods in Applied Mechanics and Engineering, 376:113606, 2021.
  * [16] DA Spence and P Sharp.  Self-similar solutions for elastohydrodynamic cavity flow.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 400(1819):289–313, 1985.
  * [17] Kirill M Terekhov, Bradley T Mallison, and Hamdi A Tchelepi.  Cell-centered nonlinear finite-volume methods for the heterogeneous anisotropic diffusion problem.  Journal of Computational Physics, 330:245–267, 2017.
  * [18] Shiqian Xu, Guotong Ren, Rami M Younis, and Qihong Feng.  Revisiting field estimates for carbon dioxide storage in depleted shale gas reservoirs: The role of geomechanics.  International Journal of Greenhouse Gas Control, 105:103222, 2021\.
  * [19] Qinglei Zeng, Zhanli Liu, Tao Wang, Yue Gao, and Zhuo Zhuang.  Fully coupled simulation of multiple hydraulic fractures to propagate simultaneously from a perforated horizontal wellbore.  Computational Mechanics, 61(1):137–155, 2018.

## Appendix A

Matrix $\bm{F}$ from eq. 14 is

$\bm{F}=\frac{\Delta t}{12\mu\Delta
x^{2}}\begin{pmatrix}\omega^{3}_{1+\frac{1}{2}}&-\omega^{3}_{1+\frac{1}{2}}&0&0&\ldots&0\\\
-\omega^{3}_{1+\frac{1}{2}}&\omega^{3}_{1+\frac{1}{2}}+\omega^{3}_{2+\frac{1}{2}}&-\omega^{3}_{2+\frac{1}{2}}&0&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&0&-\omega^{3}_{n-\frac{1}{2}}&\omega^{3}_{n-\frac{1}{2}}\end{pmatrix}$
(33)

where each row and column sum of $\bm{F}$ is $0$.

## Appendix B

Proof of Proposition 1

###### Proof.

According to eq. 16b

$\bm{A}\bm{p}^{*}+\bm{F}(\bm{p})\bm{p}^{*}=\bm{q}+\bm{\omega}^{n}$ (34)

Take the column sum of eq. 34, we get

$\sum_{i}(\bm{A}\bm{p}^{*})_{i}+\sum_{i}(\bm{F}(\bm{p})\bm{p}^{*})_{i}=\sum_{i}(\bm{q}+\bm{\omega}^{n})_{i}$
(35)

Since $\forall\bm{p}\in\mathbb{R}^{n_{c}}$,
$\sum_{i}\big{(}\bm{F}(\bm{p})\big{)}_{ij}=\sum_{j}\big{(}\bm{F}(\bm{p})\big{)}_{ij}=0$,
then $\sum_{i}(\bm{F}(\bm{p})\bm{p}^{*})_{i}=0$. Suppose the linear solver can
yield the exact solution of $\bm{p}^{*}$, then Proof 1 is valid. ∎

## Appendix C

The discretization of the proposed the coupled FVM-XFEM will be briefly
described in this section. First, flow in the fracture is handled by the
finite difference method and its discretized form is the same as eq. 15. On
the other side, the displacement approximation of XFEM is

$\displaystyle\bm{u}$ $\displaystyle=\sum_{i\in
I}N_{i}\bar{\bm{u}}_{i}+\sum_{i\in
L}N_{i}(H_{c}-H^{i}_{c})\bar{\bm{a}}_{i}+\sum_{i\in
K}\sum_{l=1}^{4}N_{i}(F_{l}-F^{i}_{l})\bar{\bm{b}}^{l}_{i}$ (36)

where $I,L,K$ is the set of standard nodes, Heaviside-enriched nodes and tip-
enriched nodes, respectively. $N_{i}$ is the shape function, $H_{c}$ is the
Heaviside function, and $F_{l}$ is the tip enrichment function. Expressions
for the function above are listed in [15]. $\bar{\bm{u}}_{i}$ is the standard
nodal displacement, $\bar{\bm{a}}_{i}$ is Heaviside-enriched nodal
displacement, and $\bar{\bm{b}}^{l}_{i}$ is tip-enriched nodal displacement.
Starting from the weak form of eq. 5

$\int_{\Omega}\delta\bm{\varepsilon}:\bm{\sigma}d\bm{x}=\int_{\Gamma_{t}}\delta\bar{\bm{u}}\cdot\bm{t}d\bm{x}+\int_{\mathcal{C}}\llbracket\delta\bar{\bm{u}}\rrbracket\cdot
p_{F}\bm{n}_{c}d\bm{x}$ (37)

where $\delta\bm{\varepsilon}$ and $\delta\bar{\bm{u}}$ are the trial
functions. Substitution of eq. 36 into the weak form obtains

$(\delta\bar{\bm{u}})^{T}\sum_{e=1}^{N}\int_{\Omega^{e}}(\bm{B}_{i}^{r})^{\rm{T}}D\bm{B}_{j}^{s}d\bm{x}\bar{\bm{u}}=(\delta\bar{\bm{u}})^{T}\sum_{e=1}^{N_{c}}\bigg{(}\int_{\mathcal{C}^{e}}N_{i}\llbracket
F_{l}\rrbracket p_{F}\bm{n}_{c}d\bm{x}+\int_{\mathcal{C}^{e}}N_{i}\llbracket
H_{\bm{\gamma}_{c}}\rrbracket
p_{F}\bm{n}_{c}d\bm{x}\bigg{)}+(\delta\bar{\bm{u}})^{T}\sum_{e=1}^{N_{\gamma}}\int_{\Gamma^{e}_{t}}N_{i}\bm{t}d\bm{x}$
(38)

where $r,s=u,a,b^{l}$, $B^{r,s}$ are the (enriched) shape functions
derivatives with respect to coordinates. $N$ is the total number of grid
blocks, $N_{c}$ is the number of fracture grids and $N_{\gamma}$ is the number
of grids along the domain outer boundaries.

The aperture $\omega$ can be derived from eq. 36

$\omega=2\sum_{i\in L}N_{i}\bar{\bm{a}}_{i}+\sum_{i\in
K}\sum_{l=1}^{4}N_{i}(F^{+}_{l}-F^{-}_{l})\bar{\bm{b}}^{l}_{i}$ (39)

eq. 26 can be obtained using the discretized form of XFEM and FVM. The
resultant Newton’s linearized system reads

$\begin{pmatrix}J_{ff}&J_{fm}\\\ J_{mf}&J_{mm}\\\
\end{pmatrix}\begin{pmatrix}\delta p\\\ \delta\bm{u}\\\
\end{pmatrix}=-\begin{pmatrix}R_{f}\\\ R_{m}\\\ \end{pmatrix},$ (40)

Where $J_{ff},J_{fm},J_{mf},J_{mm}$ are residuals’ derivatives with respect to
pressure and displacement. The quasi-newton method can be achieved by
modification of $J_{fm}$. To add contact force into the system, we adopt the
penalty method and the weak form now becomes

$\int_{\Omega}\delta\bm{\varepsilon}:\bm{\sigma}d\bm{x}=\int_{\Gamma_{t}}\delta\bar{\bm{u}}\cdot\bm{t}d\bm{x}+\int_{\mathcal{C}}\llbracket\delta\bar{\bm{u}}\rrbracket\cdot
p_{F}\bm{n}_{c}d\bm{x}+\int_{\mathcal{C}}\llbracket\delta\bar{\bm{u}}\rrbracket\cdot(-K_{ct}\omega)\bm{n}_{c}d\bm{x}$
(41)

where $K_{ct}$ is the stiffness of the contact. Hence, the discretized form
becomes

$\displaystyle(\delta\bar{\bm{u}})^{T}\sum_{e=1}^{N}\int_{\Omega^{e}}(\bm{B}_{i}^{r})^{\rm{T}}D\bm{B}_{j}^{s}d\bm{x}\bar{\bm{u}}=(\delta\bar{\bm{u}})^{T}\sum_{e=1}^{N_{c}}\bigg{(}\int_{\mathcal{C}^{e}}N_{i}\llbracket
F_{l}\rrbracket p_{F}\bm{n}_{c}d\bm{x}+\int_{\mathcal{C}^{e}}N_{i}\llbracket
H_{\bm{\gamma}_{c}}\rrbracket
p_{F}\bm{n}_{c}d\bm{x}\bigg{)}+(\delta\bar{\bm{u}})^{T}\sum_{e=1}^{N_{\gamma}}\int_{\Gamma^{e}_{t}}N_{i}\bm{t}d\bm{x}$
(42)
$\displaystyle+(\delta\bar{\bm{u}})^{T}\sum_{e=1}^{N_{ct}}\bigg{(}\int_{\mathcal{C}^{e}}N_{i}\llbracket
F_{l}\rrbracket(-K_{ct}\omega)\bm{n}_{c}d\bm{x}+\int_{\mathcal{C}^{e}}N_{i}\llbracket
H_{\bm{\gamma}_{c}}\rrbracket(-K_{ct}\omega)\bm{n}_{c}d\bm{x}\bigg{)}$ (43)

where $N_{ct}$ denotes the number of fracture elements whose $\omega$ are
negative.

## Appendix D

The fracture propagation algorithm of [15] will be briefly described in this
section. The goal is to capture the time spot $t_{c}$ when the stress
intensity factor reaches the critical value,

$K_{eq}=K_{c}$ (44)

Two time-step adjustment mechanisms are involved here,

$\Delta t^{n+1}:=\left\\{\begin{matrix}\alpha\Delta t^{n}&K_{eq}<K_{c}\\\
\frac{K_{c}-K_{eq}^{old}}{K_{eq}-K_{eq}^{old}}\Delta
t^{n}&K_{eq}>K_{c}\end{matrix}\right.$ (45)

where $\alpha$ is the growth factor. $t^{n+1}=t^{n}+\Delta t^{n}$.
$K_{eq}^{old}$ and $K_{eq}$ is recorded at $t^{n}$ and $t^{n+1}$,
respectively. For each time-step solve, the coupled system eq. 40 is solved
till eq. 44 is achieved.


# Non-singular word maps for linear groups

Henry Bradford  H. Bradford, Univ. of Cambridge, Cambridge CB3 0WB, Unites
Kingdom<EMAIL_ADDRESS>,  Jakob Schneider  J. Schneider, TU Dresden, 01062
Dresden, Germany<EMAIL_ADDRESS>and  Andreas Thom  A. Thom, TU
Dresden, 01062 Dresden, Germany<EMAIL_ADDRESS>
###### Abstract.

We study the word image of words with constants in $\operatorname{GL}(V)$ and
show that it is large provided the word satisfies some natural conditions on
its length and its critical constants.

There are various consequences: We prove that for every $l\geq 1$, there are
only finitely many pairs $(n,q)$ such that the length of the shortest non-
singular mixed identity $\operatorname{PSL}_{n}(q)$ is bounded by $l$. We
generalize the Hull–Osin dichotomy for highly transitive permutation groups to
linear groups over finite fields. Finally, we show that the rank limit of
$\operatorname{GL}_{n}(q)$ for $q$ fixed and $n\to\infty$ is mixed identity
free.

###### Contents

  1. 1 Basic notions
  2. 2 Multiply linearly transitive groups
  3. 3 Diameter and dimension of the word image
    1. 3.1 Diameter estimates
    2. 3.2 Esimates of the dimension
  4. 4 A Hull–Osin type result for linear groups
  5. 5 No mixed identities for $A(q)$

## Introduction

Denote by $\mathbf{F}_{r}=\langle x_{1},\ldots,x_{r}\rangle$ the free group on
$r$ generators. Let $G$ be a group and let $w\in G\ast\mathbf{F}_{r}$ be a
word with $r$ variables $x_{1},\dots,x_{r}$ and constants in $G$. We say that
$w$ is a _mixed identity_ for $G$ if $w$ is non-trivial and
$w(h_{1},\dots,h_{r})=1_{G}$ for all $h_{1},\dots,h_{r}\in G.$ Here we write
$w(h_{1},\ldots,h_{r})$ for the image $\varphi(w)$ of $w$ under the unique
homomorphism $\varphi\colon G\ast\mathbf{F}_{r}\to G$ which maps $x_{i}\mapsto
h_{i}$ for all $i\in\\{1,\ldots,r\\}$ and fixes $G$ elementwise. On
$G\ast\mathbf{F}_{r}$ we consider the word length, which assigns length one to
the generators $x_{1}^{\pm 1},\dots,x_{r}^{\pm 1}$ and length zero to elements
of $G$. In this setting, the _augmentation map_ $\varepsilon\colon
G\ast\mathbf{F}_{r}\to\mathbf{F}_{r}$ denotes the unique homomorphism which
maps all $G\ni g\mapsto 1_{\mathbf{F}_{r}}$ to the identity and fixes each
element of $\mathbf{F}_{r}$. The elements in the kernel of $\varepsilon$ are
called _singular_ and $\varepsilon(w)$ is called the _content_ of $w$.

In this article we continue the study of non-singular word maps on finite
groups, which we started in [schneiderthom2022word], where the second and
third named authors treated the case of symmetric and alternating groups. One
of the main results in [schneiderthom2022word] was that any non-singular mixed
identity of $S_{n}$ is of length at least $\Omega(\log(n)/\log(\log(n))).$
This is in contrast to the length of singular mixed identities, which can be
of bounded length. See [bradfordschneiderthom2023length] for more information
on mixed identities in this context. Non-singular word maps for compact Lie
groups where studied in [klyachkothom2017new] and for linear algebraic groups
in [gordeevkunyavskiiplotkin2018wordmaps].

A _law_ or an _identity_ of a group $G$ is a non-trivial word
$w\in\mathbf{F}_{r}$ such that $w(h_{1},\dots,h_{r})=1_{G}$ for all
$h_{1},\dots,h_{r}\in G$. We say that $w\in\mathbf{F}_{r}$ is a _coset
identity_ if there exists a subgroup of finite index $H$ and cosets
$g_{1}H,\dots,g_{r}H$, such that $w(h_{1},\dots,h_{r})=1_{G}$ for all
$h_{i}\in g_{i}H,i\in\\{1,\dots,r\\}$. It is a major open problem to decide if
every finitely generated pro-finite group that admits a coset identity also
admits an identity, see for example [mazurovkhukhro2014kourovka]*Problem 12.95
and the work of Larsen–Shalev [larsenshalev2016probabilistic]*Problem 3.1.
Note that every coset identity yields a non-singular mixed identity

$w^{\prime}(x_{1},\dots,x_{r})\coloneqq
w(g_{1}x_{1}^{d},\dots,g_{r}x_{r}^{d})$

for $d=[G:H]!$. Hence, a potential strategy to answer Larsen and Shalev’s
question is to study the length of the shortest non-singular mixed identity of
a finite group $G$ and compare it to the length of the shortest identity. It
can be shown that the answer to Larsen and Shalev’s question is positive if
these quantities can be bounded in terms of each other.

######  Question 1.

Does the existence of a non-singular mixed identity $w\in G\ast\mathbf{F}_{r}$
for a finite group $G$ imply the existence of a law for $G$ of length bounded
only in terms of the length of $w$?

We are able to prove that the length of the shortest non-singular mixed
identity for $\operatorname{PSL}_{n}(q)$ tends to infinity as $n$ or $q$
increase. Note, that the minimal length of laws for groups
$\operatorname{PSL}_{n}(q)$ has been studied in [bradfordthom2018laws]. Hence,
this contributes another class of cases where this question has a positive
answer, see Corollary 3.

In view of the corresponding result of Jones [jones1974varieties] on laws for
non-abelian finite simple groups, we put forward the following conjecture:

######  Conjecture 1.

Let $l\in\mathbb{N}$. The number of isomorphism classes of non-abelian finite
simple groups satisfying a non-singular mixed identity of length $\leq l$ is
finite.

Similarly to our results for symmetric groups, our analysis of word maps
allows also for detailed information on the image of word maps, singular or
non-singular, as long as the lengths of the critical constants (see Section 1
for a definition) are under certain control. Our main result is stated as
Theorem 1. We have various applications of this result in analogy to results
on permutation actions. In particular, we have estimates on the diameter and
dimension of the image.

It was shown by Gordeev–Kunyavskii–Plotkin in
[gordeevkunyavskiiplotkin2018wordmaps] that non-singular word maps on linear
algebraic groups with generic constants are dominant. As a consequence of our
results, we can prove a non-trivial lower bound on the dimension of the image
in case of ${\rm SL}_{n}$ and non-singular word maps with arbitrary constants,
see Theorem 2.

In Section 4, we can prove an analogue of the Hull–Osin dichotomy
[hullosin2016transitivity]*Theorem 1.6 for highly linearly transitive groups,
which in its original form says that a highly transitive permutation group
either contains a locally finite, highly transitive, normal subgroup or it is
mixed identity free. We prove a similar result for highly linearly transitive
groups acting on a countably infinite dimensional vector space over a finite
field, see Theorem 5. There is also a quantitative version of this,
generalizing a result of Le Boudec–Matte Bon
[leboudecmattebon2022triple]*Proposition A.1, see Theorem 6.

It is a consequence of the original Hull–Osin dichotomy that full groups of
ergodic, probability measure preserving, measurable equivalence relations on a
non-atomic probability space are mixed identity free, see the discussion in
Section 5. Moreover, it follows from Popa’s work [popa1995free] that the
unitary group of a II1-factor does not satisfy a mixed identity. No such
approach can be taken to understand the mod-$p$ analogue of the unitary group
of the hyperfinite II1-factor as introduced by Carderi and the third author in
[carderithom2018exotic]. However, we give a direct proof that this group is
also mixed identity free, see Theorem 7.

Finally, this paper might be the starting point of a more systematic study of
mixed identities for oligomophic permutation groups such as ${\rm
Sym}(\mathbb{N})$ of $\operatorname{GL}(\mathbb{F}_{q}^{\oplus\omega})$. Note
that ${\rm Aut}(\mathbb{Q},<)$ satisfies a mixed identity as consequence of
work of Zarzycki [zarzycki2010limits], while the automorphism group of the
Rado graph and the universal Urysohn metric space does not, see
Etedadialiabadi–Gao–Le Maître-Melleray
[etedadialiabadigaolemaitremelleray2021dense]. See also Remark 3 for further
discussion.

##  1\. Basic notions

Let’s recall some basics where we follow the notation introduced in
[schneiderthom2022word]. Let $G$ be a group. Throughout the entire article,
let

$w=c_{0}x_{\iota(1)}^{\varepsilon(1)}c_{1}\cdots
c_{l-1}x_{\iota(l)}^{\varepsilon(l)}c_{l}\in G\ast\mathbf{F}_{r}=G\ast\langle
x_{1},\ldots,x_{r}\rangle$

be a fixed word with constants $c_{j}\in G$ (for $j\in\\{0,\ldots,l\\}$) with
$\varepsilon(j)\in\\{\pm 1\\}$ and $\iota(j)\in\\{1,\ldots,r\\}$ (for
$j\in\\{1,\ldots,l\\}$). Subsequently, $H\leq G$ will be the group whose
elements we plug into $w$. Set
$J_{0}(w)\coloneqq\\{j\in\\{1,\ldots,l-1\\}\,|\,\iota(j)\neq\iota(j+1)\\}$,
$J_{+}(w)\coloneqq\\{j\in\\{1,\ldots,l-1\\}\,|\,\iota(j)=\iota(j+1)\text{ and
}\varepsilon(j)=\varepsilon(j+1)\\}$, and
$J_{-}(w)\coloneqq\\{j\in\\{1,\ldots,l-1\\}\,|\,\iota(j)=\iota(j+1)\text{ and
}\varepsilon(j)=-\varepsilon(j+1)\\}$. We call the elements of $J_{-}(w)$ the
_critical indices_ and the group elements $c_{j}$ for $j\in J_{-}(w)$ the
_critical constants_ of $w$. Furthermore, we assume that the above expression
for $w$ is _reduced_ , i.e. $c_{j}\notin\mathbf{C}_{G}(H)\setminus\mathbf{1}$
for $j\in J_{0}(w)\cup J_{+}(w)$ and $c_{j}\notin\mathbf{C}_{G}(H)$ for all
$j\in J_{-}(w)$. Here $\mathbf{C}_{G}(H)$ denotes the _centralizer_ of $H$ in
$G$. We call the word $w$ _strong_ if it has no critical constants, i.e.
$J_{-}(w)=\emptyset$. And we call it _cyclically strong_ if $w$ is strong and,
in addition, it does not hold that $\iota(1)=\iota(l)$ and
$\varepsilon(1)=-\varepsilon(l)$. Finally, we write $\lvert w\rvert\coloneqq
l$ for the _length_ of the word $w$.

Throughout, we will need some measure of distance on our groups. That is why
we define the notion of a _seminorm_ on a group. For a group $G$, this is a
mapping $\lvert\bullet\rvert\colon G\to I$, where $I$ is either the unit
interval $[0,1]\subseteq\mathbb{R}$ or the set $\mathbb{N}\cup\\{\infty\\}$,
such that $\lvert g\rvert\geq 0$ (_non-negativity_), $\lvert
gh\rvert\leq\lvert g\rvert+\lvert h\rvert$ (_subadditivity_), and $\lvert
g\rvert=\lvert g^{-1}\rvert$ (_homogeneity_) for all $g,h\in G$. Each seminorm
induces a left invariant _pseudometric_ $\operatorname{dist}\colon G\times
G\to I$ by $\operatorname{dist}(g,h)\coloneqq\lvert g^{-1}h\rvert.$ A seminorm
becomes a norm if $\lvert g\rvert=0$ implies $g=1_{G}$ for all $g\in G$
(_positivity_). In this case, it induces a metric. For each seminorm
$\lvert\bullet\rvert$, the set $G_{0}\coloneqq\lvert\bullet\rvert^{-1}(0)\leq
G$ is a subgroup. The seminorm is called _invariant_ if it is invariant under
conjugation, i.e. $\lvert g^{h}\rvert=\lvert h^{-1}gh\rvert=\lvert g\rvert$
for all $g,h\in G$. In this case, $G_{0}\trianglelefteq G$ is a normal
subgroup; then $\lvert\bullet\rvert$ descends to a norm on $G/G_{0}$. Then,
the corresponding pseudometric $\operatorname{dist}\colon G\times G\to I$ will
be _bi-invariant_ , i.e.
$\operatorname{dist}(fg,fh)=\operatorname{dist}(g,h)=\operatorname{dist}(gf,hf)$
for all $f,g,h\in G$. In this article, all seminorms will be invariant, whence
all pseudometrics will be bi-invariant. For every seminorm
$\lvert\bullet\rvert\colon G\to I$ and a normal subgroup $N\trianglelefteq G$,
we define the _quotient seminorm_ $\lvert\bullet\rvert\colon G/N\to I$ by
$\lvert\overline{g}\rvert\coloneqq\operatorname{dist}(g,N)=\inf\\{\operatorname{dist}(g,n)\,|\,n\in
N\\}$ for all $g\in G$.

For $w\in G\ast\mathbf{F}_{r}$ and a fixed seminorm on $G$, we define the
_critical length_

$\lvert w\rvert_{\rm crit}\coloneqq\min\\{\operatorname{diam}(G),\lvert
c_{j}\rvert\,|\,j\in J_{-}(w)\\},$

where the diameter is measured with respect to the fixed seminorm. Besides the
word length of a word, the critical length turns out to be the most important
quantity controlling lower bounds on the size of the word image.

##  2\. Multiply linearly transitive groups

Let $V$ be a vector space over the field $K$ and of finite or countably
infinite dimension. We only consider the cases when $K=\mathbb{F}_{q}$ is a
finite field, or $K=\mathbb{C}$ are the complex numbers. Write
$\operatorname{GL}(V)$ for the general linear group of $V$. Denote by
$\lvert\bullet\rvert\colon\operatorname{GL}(V)\to\mathbb{N}\cup\\{\infty\\}$
the _projective rank seminorm_ , which is defined by

$\lvert g\rvert\coloneqq\min_{\lambda\in
K^{\times}}\operatorname{rk}(g-\lambda 1_{V})\in\mathbb{N}\cup\\{\infty\\}$

for $g\in\operatorname{GL}(V)$. It descends to a _norm_ , which we also denote
by

$\lvert\bullet\rvert\colon\operatorname{PGL}(V)\to\mathbb{N}\cup\\{\infty\\}.$

This induces the _projective rank pseudometric_ on $\operatorname{GL}(V)$ and
a _projective rank metric_ on $\operatorname{PGL}(V)$, both denoted by
$\operatorname{dist}$.

######  Definition 1.

Let $V$ be a vector space of finite or countably infinite dimension and let
$k\in\mathbb{N}$. A linear group $H\leq\operatorname{GL}(V)$ is called
_linearly $k$-transitive_ if for all $k$-element sets of linearly independent
vectors $\\{u_{1},\ldots,u_{k}\\}$ and $\\{w_{1},\ldots,w_{k}\\}$ there is a
map $h\in H$ such that $u_{i}.h=w_{i}$ for all $i\in\\{1,\ldots,k\\}$.

Here is our main technical result:

######  Theorem 1.

Let $V$ be some vector space of finite or countably infinite dimension over
the field $K$. If $K$ is finite, we allow $\dim(V)$ to be infinite, whereas
for $K=\mathbb{C}$, we require it to be finite. Let $D,d,l,r\in\mathbb{N}$ and
let the word

$w=x_{\iota(1)}^{\varepsilon(1)}c_{1}\cdots
c_{l-1}x_{\iota(l)}^{\varepsilon(l)}\in\operatorname{GL}(V)\ast\mathbf{F}_{r}$

be reduced of length $l\geq 2$. Let $H\leq\operatorname{GL}(V)$ be a linearly
$D$-transitive group. Let $\\{u_{1},\ldots,u_{d}\\}$ and
$\\{w_{1},\ldots,w_{d}\\}$ be $d$-element sets of linearly independent
vectors. If $\iota(1)=\iota(l)$ and $\varepsilon(1)=-\varepsilon(l)$, then
assume in addition that the linear spans of $\\{u_{1},\ldots,u_{d}\\}$ and
$\\{w_{1},\ldots,w_{d}\\}$ intersect trivially. Then, if $d$ satisfies

$d\leq\frac{1}{l}\cdot\min\\{\lvert w\rvert_{\rm crit}-1,D\\},$

there exist $h_{1},\ldots,h_{r}\in H$ such that
$u_{i}.w(h_{1},\ldots,h_{r})=w_{i}$ for $i\in\\{1,\ldots,d\\}$.

###### Proof.

We set $U=\langle u_{1},\ldots,u_{d}\rangle$ and $W\coloneqq\langle
w_{1},\ldots,w_{d}\rangle.$ In order to prove the theorem, we will define
vectors $v_{i,j}^{-\varepsilon(j)}\in V$ for
$(i,j)\in\\{1,\dots,d\\}\times\\{1,\dots,l\\}$ and elements
$h_{1},\ldots,h_{r}\in H$ such that

$\displaystyle u_{1}=v_{1,1}^{\varepsilon(1)}\stackrel{{\scriptstyle
h^{\varepsilon(1)}_{\iota(1)}}}{{\mapsto}}v_{1,1}^{-\varepsilon(1)}\stackrel{{\scriptstyle
c_{1}}}{{\mapsto}}v_{1,2}^{\varepsilon(2)}\stackrel{{\scriptstyle
h^{\varepsilon(2)}_{\iota(2)}}}{{\mapsto}}\cdots\stackrel{{\scriptstyle
h^{\varepsilon(l-1)}_{\iota(l-1)}}}{{\mapsto}}v_{1,l-1}^{-\varepsilon(l-1)}\stackrel{{\scriptstyle
c_{l-1}}}{{\mapsto}}v_{1,l}^{\varepsilon(l)}\stackrel{{\scriptstyle
h^{\varepsilon(l)}_{\iota(l)}}}{{\mapsto}}v_{1,l}^{-\varepsilon(l)}=w_{1}$
$\displaystyle\vdots$ $\displaystyle
u_{d}=v_{d,1}^{\varepsilon(1)}\stackrel{{\scriptstyle
h^{\varepsilon(1)}_{\iota(1)}}}{{\mapsto}}v_{d,1}^{-\varepsilon(1)}\stackrel{{\scriptstyle
c_{1}}}{{\mapsto}}v_{d,2}^{\varepsilon(2)}\stackrel{{\scriptstyle
h^{\varepsilon(2)}_{\iota(2)}}}{{\mapsto}}\cdots\stackrel{{\scriptstyle
h^{\varepsilon(l-1)}_{\iota(l-1)}}}{{\mapsto}}v_{d,l-1}^{-\varepsilon(l-1)}\stackrel{{\scriptstyle
c_{l-1}}}{{\mapsto}}v_{d,l}^{\varepsilon(l)}\stackrel{{\scriptstyle
h^{\varepsilon(l)}_{\iota(l)}}}{{\mapsto}}v_{d,l}^{-\varepsilon(l)}=w_{d}.$

This requirement will determine vectors $v^{\varepsilon(j)}_{i,j}\in V$ by the
equation $v_{i,j}^{\varepsilon(j)}\coloneqq
v_{i,j-1}^{-\varepsilon(j-1)}.c_{j-1}$ for $2\leq j\leq l$ and
$v^{\varepsilon(1)}_{i,1}\coloneqq u_{i}.$ The choice of the vectors
$v^{-\varepsilon(j)}_{i,j}\in V$ proceeds by induction on the lexicographic
order $\leq$ on the set $\\{1,\ldots,d\\}\times\\{1,\ldots,l\\}$. During the
procedure, we want to ensure that the sets
$\\{v_{i,j}^{\varepsilon}\,|\,i\in\\{1,\ldots,d\\},j\in\\{1,\ldots,l\\},\iota(j)=k\\}$
are linearly independent for all $\varepsilon\in\\{\pm 1\\}$ and
$k\in\\{1,\dots,r\\}$ and we can put $v_{i,l}^{-\varepsilon(l)}\coloneqq
w_{i}$ at the end of the procedure.

Once this is done, we can choose $h_{k}\in H$ such that
$v^{+}_{i,j}.h_{k}=v^{-}_{i,j}$ whenever $\iota(j)=k$ and the picture is
complete. For this, it is enough that $H$ is linearly $dl$-transitive which we
assumed to be the case. Thus, in order to finish the proof, it only remains to
carry out the inductive process and keep track of the linear independence.

Before we start with the induction, let’s first describe some auxiliary
objects which are needed in the inductive procedure. Suppose that all
$v^{\varepsilon}_{i^{\prime},j^{\prime}}$ with
$(i^{\prime},j^{\prime})\leq(i,j)$ are already defined. We define subspaces
for $i\in\\{1,\ldots,d\\}$, $j\in\\{1,\ldots,l\\}$, $k\in\\{1,\ldots,r\\}$,
and $\varepsilon\in\\{\pm\\}$.

$\displaystyle V_{i,j,k}^{\varepsilon}$ $\displaystyle\coloneqq\langle
v_{i^{\prime},j^{\prime}}^{\varepsilon}\,|\,(i^{\prime},j^{\prime})\leq(i,j),\iota(j^{\prime})=k\rangle$
$\displaystyle\leq\sum_{k=1}^{r}{V_{i,j,k}^{\varepsilon}}=\langle
v_{1,1}^{\varepsilon},\ldots,v_{1,l}^{\varepsilon},\ldots,v_{i,1}^{\varepsilon},\ldots,v_{i,j}^{\varepsilon}\rangle\eqqcolon
V_{i,j}^{\varepsilon}\leq V,$

so that $\dim(V_{i,j}^{\varepsilon})\leq(i-1)l+j$. We also set
$V_{i,0,k}^{\varepsilon}\coloneqq V_{i-1,l,k}^{\varepsilon}$ for
$i\in\\{2,\ldots,d\\}$, $k\in\\{1,\ldots,r\\}$, $\varepsilon\in\\{\pm\\}$. For
$\varepsilon\in\\{\pm\\}$ and $\iota\in\\{1,\ldots,r\\}$, we need to introduce
the subspaces:

$\displaystyle V_{\iota}^{\varepsilon}\coloneqq\begin{cases}U&\text{if
}(\iota,\varepsilon)=(\iota(1),\varepsilon(1))\text{ and
}(\iota,\varepsilon)\neq(\iota(l),-\varepsilon(l))\\\ W&\text{if
}(\iota,\varepsilon)\neq(\iota(1),\varepsilon(1))\text{ and
}(\iota,\varepsilon)=(\iota(l),-\varepsilon(l))\\\ U+W&\text{if
}(\iota,\varepsilon)=(\iota(1),\varepsilon(1))\text{ and
}(\iota,\varepsilon)=(\iota(l),-\varepsilon(l))\\\
\mathbf{0}&\text{otherwise.}\end{cases}$

which we will need later in the proof.

We will now describe how to choose the vectors $v_{i,j}^{-\varepsilon(j)}$. In
the $(i,j)$th step of the induction with $j<l$, all the vectors
$v_{i^{\prime},j^{\prime}}^{\varepsilon}$ for $(i^{\prime},j^{\prime})<(i,j)$
and $\varepsilon\in\\{\pm\\}$, and the vector $v_{i,j}^{\varepsilon(j)}$, are
already defined. Hence the subspace
$V_{i^{\prime},j^{\prime},k}^{\varepsilon}$ for
$(i^{\prime},j^{\prime})<(i,j)$ for $k\in\\{1,\ldots,r\\}$ and
$\varepsilon\in\\{\pm\\}$ and the subspace
$V_{i,j,\iota(j)}^{\varepsilon(j)}=V_{i,j-1,\iota(j)}^{\varepsilon(j)}+\langle
v_{i,j}^{\varepsilon(j)}\rangle$ are already defined. The vector
$v_{i,j}^{-\varepsilon(j)}$ is not yet defined and our aim is to construct it
now.

In order to retain the linear independence condition, we need to ensure that

$v_{i,j}^{-\varepsilon(j)}\notin
V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+V_{\iota(j)}^{-\varepsilon(j)}\quad\mbox{and}\quad
v_{i,j+1}^{\varepsilon(j+1)}=v_{i,j}^{-\varepsilon(j)}.c_{j}\notin
V_{i,j,\iota(j+1)}^{\varepsilon(j+1)}+V_{\iota(j+1)}^{\varepsilon(j+1)}.$

For the first constraint $v_{i,j}^{-\varepsilon(j)}\notin
V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+V_{\iota(j)}^{-\varepsilon(j)}$, note
that this subspace is contained in $V_{d,l-2}^{-\varepsilon(j)}+\langle
w_{d}\rangle$ and so its dimension is at most $dl-1$. The case $l=2,j=1$ needs
a special argument, but note that $V_{\iota(1)}^{-\varepsilon(1)}\leq W$ by
its definition and hence we obtain the same bound also in this case. In order
to study the second constraint, we will now distinghish two cases.

_Case 1._ Assume $j\in J_{0}(w)\cup J_{+}(w)$. We must ensure in addition that

$v_{i,j}^{-\varepsilon(j)}.c_{j}\notin
V_{i,j,\iota(j+1)}^{\varepsilon(j+1)}+V_{\iota(j+1)}^{\varepsilon(j+1)}.$

Again, this condition excludes a subspace of dimension at most $dl-1$. Indeed,
we have

$V_{i,j,\iota(j+1)}^{\varepsilon(j+1)}+V_{\iota(j+1)}^{\varepsilon(j+1)}\leq\begin{cases}V_{d,l-2}^{\varepsilon(j+1)}+\langle
w_{d}\rangle&1\leq j\leq l-2\\\ V^{\varepsilon(l)}_{d,l-1}&j=l-1,\end{cases},$

since $V_{\iota(l)}^{\varepsilon(l)}\leq U$ by definition in the second case,
which yields the required bound on the dimension in each case.

Dealing with both constraints, we must have

$v_{i,j}^{-\varepsilon(j)}\notin(V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+V_{\iota(j)}^{-\varepsilon(j)})\cup(V_{i,j,\iota(j+1)}^{\varepsilon(j+1)}+V_{\iota(j+1)}^{\varepsilon(j+1)}).c_{j}^{-1}$

and this is a union of two subspaces of dimension at most $dl-1$. Since
$dl-1<dl<\lvert w^{\prime}\rvert_{\rm crit}\leq\dim(V)$, these subspaces have
codimension at least $2$ and we find an appropriate choice for
$v_{i,j}^{-\varepsilon(j)}$ in this case.

_Case 2._ Assume now that $j\in J_{-}(w)$. We have to ensure that

$\displaystyle v_{i,j}^{-\varepsilon(j)}.c_{j}\notin
V_{i,j,\iota(j+1)}^{\varepsilon(j+1)}+V_{\iota(j+1)}^{\varepsilon(j+1)}$
$\displaystyle=V_{i,j-1,\iota(j)}^{\varepsilon(j+1)}+\langle
v_{i,j}^{\varepsilon(j+1)}\rangle+V_{\iota(j+1)}^{\varepsilon(j+1)}$
$\displaystyle=V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+\langle
v_{i,j}^{-\varepsilon(j)}\rangle+V_{\iota(j+1)}^{\varepsilon(j+1)}.$

In order to get this, it is enough to ensure that for each $\lambda\in K$

$v_{i,j}^{-\varepsilon(j)}.(c_{j}-\lambda)\notin
V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+V_{\iota(j+1)}^{\varepsilon(j+1)}.$

Note that this subspace is again of dimension bounded by $dl-1$, where the
case $l=2$ requires an additional argument that we omit. We further
distinguish three more cases:

_Case 2.1:_ Let’s assume that $\dim(V)=n<\infty$ and let’s first deal with the
case $K=\mathbb{F}_{q}$. Let $\lambda_{k}\in\mathbb{F}_{q}^{\times}$ for
$k\in\\{1,\ldots,m\\}$ be the eigenvalues of the linear map $c_{j}$ and
$d_{k}$ their respective geometric multiplicity. Then for
$\lambda\notin\\{\lambda_{1},\dots,\lambda_{m}\\}$, $c_{j}-\lambda$ is
regular, so the condition above excludes a union of $q-m$ subspaces of
dimension at most $dl-1$. Thus, we need to exclude at most $(q-m)q^{dl-1}$
vectors.

For $\lambda=\lambda_{k}$ we have to exclude the inverse image of
$V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+V_{\iota(j+1)}^{\varepsilon(j+1)}$ under
the map $c_{j}-\lambda.$ Thus, to deal with $\lambda_{k}$, we need to exclude
at most $q^{d_{k}}q^{dl-1}$ vectors. Hence, when dealing with all
$\lambda\in\mathbb{F}_{q}$, we need to avoid at most

$(q-m)q^{dl-1}+\sum_{k=1}^{m}{q^{dl-1}q^{d_{k}}}$

vectors. This covers the second constraint and we need to exclude at most
$q^{dl-1}$ vectors for the first constraint. Hence, in order to show that
there is a good choice for $v_{i,j}^{-\varepsilon(j)}$ it remains to show that

$(q-m+1)q^{dl-1}+\sum_{k=1}^{m}{q^{dl-1}q^{d_{k}}}<q^{n}-1.$

By assumption $n-d_{k}\geq|c_{j}|\geq|w|_{\rm crit}>dl.$ This implies that

$\displaystyle(q-m+1)q^{dl-1}+\sum_{k=1}^{m}{q^{dl-1}q^{d_{k}}}$
$\displaystyle\leq$ $\displaystyle(q-m+1)q^{n-2}+mq^{n-2}<q^{n}-1.$

Hence, we have shown that there is a good choice for
$v_{i,j}^{-\varepsilon(j)}$.

_Case 2.2:_ The proof for $K=\mathbb{C}$ is similar. Consider

$V_{\lambda}\coloneqq\left\\{v\in V\mid v.(c_{j}-\lambda)\in
V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+V_{\iota(j+1)}^{\varepsilon(j+1)}\right\\}$

and note that the estimate $n-d_{k}\geq|c_{j}|\geq|w|_{\rm crit}>dl$ is still
valid and provides a bound of $n-2$ on the dimension of each $V_{\lambda_{k}}$
for $\lambda_{k}$ an eigenvalue of $c_{j}$. In order to deal with the other
choices of $\lambda$, i.e. $\lambda\coloneqq
D=\mathbb{C}\setminus\\{\lambda_{1},\dots,\lambda_{m}\\},$ note that
$V_{\lambda}$ is a $1$-dimensional algebraic family of subspaces of
codimension at least $2$. Thus, $\bigcup_{\lambda\in D}V_{\lambda}$ is of
codimension at least $1$. Thus, we have to avoid finitely many algebraic
subsets of lower dimension, which is of course possible.

_Case 2.3:_ If $K=\mathbb{F}_{q}$ and $\dim(V)=\infty$, there are still only
finitely elements $\lambda\in\mathbb{F}_{q}$ to consider and the inverse image
of $V_{i,j-1,\iota(j)}^{-\varepsilon(j)}+V_{\iota(j+1)}^{\varepsilon(j+1)}$
under the map $c_{j}-\lambda$ is still proper. In fact, together with the
first constraint there are at most $q+1$ of these subspaces to consider and
they are all of codimension at least $2$. This implies that their union has a
non-trivial complement from which we can choose $v_{i,j}^{-\varepsilon(j)}$.
Thus we are done for this case as well.

Thus, this procedure can be carried out until we arrive at $j=l$ and then we
put $v_{i,l}^{-\varepsilon(l)}\coloneqq w_{i}.$ This finishes the proof of the
main technical result. ∎

######  Remark 1.

For finite permutation groups $H\leq S_{n}$, it is a well-known consequence of
the CFSG that the only examples (other than $A_{n}$ and $S_{n}$) of
$k$-transitive groups for $k\geq 4$ are the four Matthieu groups
$M_{11},M_{12},M_{23}$, and $M_{24}$ in degrees $n=11,12,23$ and $24$. There
is a similar classification of multiply linearly transitive subgroups of
$\operatorname{GL}_{n}(q)$, following from the work of Hering
[hering1985transitive]. Indeed the linearly $1$-transitive finite groups (that
is, those acting transitively on the set of non-zero vectors in
$V=\mathbb{F}_{q}^{n}$) are already rather restricted: apart from a finite
list of sporadic examples, and overgroups of $\operatorname{SL}_{n}(q)$, they
all preserve either a non-trivial symplectic form on $V$ or a proper field
extension (meaning a structure on $V$ of an $(n/m)$-dimensional
$\mathbb{F}_{q^{m}}$-vector space, for some proper divisor $m$ of $n$), both
of which are incompatible with linear $2$-transitivity. The sporadic cases
include one linearly $3$-transitive example (an action of $A_{7}$ on
$\mathbb{F}_{2}^{4}$); none of the others are linearly $2$-transitive
[lee2023email].

The rest of the paper consists of consequences that we will explain in some
detail.

##  3\. Diameter and dimension of the word image

###  3.1. Diameter estimates

From Theorem 1 we get the following corollary on the diameter bound of the
word image.

######  Corollary 1.

Let $n\in\mathbb{N}$. Let $V$ be an $n$-dimensional vector space over
$K=\mathbb{F}_{q}$ or $\mathbb{C}$ and
$w\in\operatorname{GL}(V)\ast\mathbf{F}_{r}$ be a reduced word of length
$l\geq 2$. Then the word image $w(\operatorname{SL}(V)^{r})$ has diameter at
least

$\operatorname{diam}(w(\operatorname{SL}(V)^{r}))\geq\frac{\lvert w\rvert_{\rm
crit}}{l}-1$

in the projective rank seminorm.

Equivalently,

(1) $\lvert
w\rvert(\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+1)\geq\lvert
w\rvert_{\rm crit}.$

Let’s now assume that $w$ is non-singular and try to get a lower bound on the
diameter of the word image only depending on the word length and not the
critical length. The argument follows closely the corresponding argument for
permutation groups in [schneiderthom2022word]. If we remove the smallest
critical constant $c_{j}$ in $w$ and perform the corresponding cancellation,
we get a word $w^{\prime}$ such that

$\operatorname{dist}(w(h_{1},\ldots,h_{r}),w^{\prime}(h_{1},\ldots,h_{r}))\leq\operatorname{dist}(c_{j},1)=\lvert
w\rvert_{\rm crit}$

for all $h_{1},\ldots,h_{r}\in H$ by bi-invariance of the rank metric. Hence

$\lvert\operatorname{diam}(w(\operatorname{SL}(V)^{r}))-\operatorname{diam}(w^{\prime}(\operatorname{SL}(V)^{r}))\rvert\leq
2\lvert w\rvert_{\rm crit}.$

so by (1) we get

$\displaystyle\operatorname{diam}(w^{\prime}(\operatorname{SL}(V)^{r}))$
$\displaystyle\leq\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+2\lvert
w\rvert_{\rm crit}$
$\displaystyle\leq\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+2\lvert
w\rvert(\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+1).$

Therefore, we obtain

$\operatorname{diam}(w^{\prime}(\operatorname{SL}(V)^{r}))+1\leq(1+2\lvert
w\rvert)(\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+1)$

Hence, if we iterate these reductions and $w=w_{0},\ldots,w_{m}$ is a chain of
words such that $w_{j}$ is an reduction of $w_{j-1}$ ($1\leq j\leq m$) and
$w_{m}$ is _strong_ , then we get

$\operatorname{diam}(w_{j}(\operatorname{SL}(V)^{r}))+1\leq(1+2\lvert
w_{j-1}\rvert)(\operatorname{diam}(w_{j-1}(\operatorname{SL}(V)^{r}))+1)$

for all $j=1,\ldots,m$, and by iterating we obtain

$\displaystyle\operatorname{diam}(w_{m}(\operatorname{SL}(V)^{r}))+1$
$\displaystyle\leq(\operatorname{diam}(w_{0}(\operatorname{SL}(V)^{r}))+1)\prod_{j=0}^{m-1}(1+2\lvert
w_{j}\rvert)$
$\displaystyle\leq(\operatorname{diam}(w_{0}(\operatorname{SL}(V)^{r}))+1)(1+2\lvert
w_{0}\rvert)^{m}$
$\displaystyle\leq(\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+1)(1+2\lvert
w\rvert)^{\left\lfloor\lvert w\rvert/2\right\rfloor}.$

This is since $\lvert w_{j-1}\rvert\geq\lvert w_{j}\rvert-2$ so that
$m\leq\left\lfloor\lvert w\rvert/2\right\rfloor$. Hence by Inequality (1) we
get

$\displaystyle n=\lvert w_{m}\rvert_{\rm crit}$
$\displaystyle\leq(\operatorname{diam}(w_{m}(\operatorname{SL}(V)^{r}))+1)\lvert
w_{m}\rvert$
$\displaystyle\leq(\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+1)(1+2\lvert
w\rvert)^{\left\lfloor\lvert w\rvert/2\right\rfloor}\lvert w\rvert.$

This gives

$\frac{1}{(1+2\lvert w\rvert)^{\left\lfloor\lvert
w\rvert/2\right\rfloor}\lvert
w\rvert}\leq\frac{\operatorname{diam}(w(\operatorname{SL}(V)^{r}))+1}{n}.$

In particular, this implies that when $\lvert w\rvert$ is bounded and
$n\to\infty$, then also
$\operatorname{diam}(w(\operatorname{SL}(V)^{r}))\to\infty$.

######  Corollary 2.

Let $q$ be a prime power. Let $w$ be a non-singular mixed identity for ${\rm
PSL}_{n}(q)$. Then, we have

$|w|=\Omega\left(\frac{\log(n)}{\log\log(n)}\right).$

######  Corollary 3.

For every $l\geq 1$, there are only finitely many pairs $(n,q)$ such that the
length of the shortest non-singular mixed identity $\operatorname{PSL}_{n}(q)$
is bounded by $l$.

###### Proof.

It is known from [bradfordschneiderthom2023length] that any mixed identity for
$\operatorname{PSL}_{n}(q)$ is of length $\Omega(q)$, where the implied
constant is independent of $n$. This and the previous corollary imply the
claim. ∎

######  Remark 2.

By the proof of [schneiderthom2022word]*Corollary 3, there exists a singular
word $v\in\operatorname{PSL}_{n}(q)\ast\langle x\rangle$ of length
$O(q^{3n^{2}})$ with $v(g)=g^{-1}$ for all $g\in\operatorname{PSL}_{n}(q)$.
Clearly, $w=vx$ is a non-singular word with the property $w(g)=1$ for all
$g\in\operatorname{PSL}_{n}(q)$ and content $\varepsilon(w)=x$. This shows
that there is no reason to assume that non-singularity by itself or just the
length of the content of a word would be enough to give a lower bound on the
diameter of the word image.

###  3.2. Esimates of the dimension

In this section, we follow the notation from the work of
Gordeev–Kunyavskii–Plotkin [gordeevkunyavskiiplotkin2018wordmaps] where we
stick to the special case $G=\operatorname{SL}_{n}$. Their main result is
[gordeevkunyavskiiplotkin2018wordmaps]*Theorem 1.1, which states that for a
fixed non-singular equation

$w=c_{0}x_{\iota(1)}^{\varepsilon(1)}c_{1}\cdots
c_{l-1}x_{\iota(l)}^{\varepsilon(l)}c_{l}\in G\ast\mathbf{F}_{r}=G\ast\langle
x_{1},\ldots,x_{r}\rangle$

with $(c_{0},\dots,c_{l})\in{\rm SL}_{n}(\mathbb{C})^{l+1}$, there is a
Zariski dense set of choices for the constants, where image is full-
dimensional. This is essentially a consequence of the fact that the maximal
dimension is generically attained and the image is full-dimensional by an old
result of Borel [borel1983on] for the trivial choice of constants. However,
their result says nothing about the worst case scenario.

######  Theorem 2.

Let $w=c_{0}x_{\iota(1)}^{\varepsilon(1)}c_{1}\cdots
c_{l-1}x_{\iota(l)}^{\varepsilon(l)}c_{l}\in\operatorname{GL}_{n}(\mathbb{C})\ast\mathbf{F}_{r}$
be reduced of length $l\geq 2$. Then, we have

$\dim w(\operatorname{SL}_{n}(\mathbb{C}))\geq\left\lfloor\frac{\lvert
w\rvert_{\rm crit}-1}{l}\right\rfloor\cdot n.$

###### Proof.

Note that $d=\left\lfloor l^{-1}\lvert w\rvert_{\rm crit}-1\right\rfloor$ is a
valid choice for $d$ in Theorem 1. Hence, by linear $d$-transitivity of the
image of $w$, we get that the space of $(n\times d)$-matrices of rank $d$ is a
subquotient of the image of $w$. This finishes the proof since its dimension
is $dn$. ∎

##  4\. A Hull–Osin type result for linear groups

The following dichotomy theorem is due to Hull and Osin
[hullosin2016transitivity]*Theorem 1.6. Let’s recall it quickly before we
explain our generalization to the linear case.

######  Theorem 3 (Hull–Osin).

Let $\Omega$ be a countably infinite set and let
$H\leq\operatorname{Sym}(\Omega)$ be a highly transitive permutation group.
Then either

$\operatorname{Alt}(\Omega)\trianglelefteq H\quad\mbox{or}\quad\mbox{$H$ has
no mixed identities.}$

There is also a quantitative version of this theorem due to Le Boudec–Matte
Bon [leboudecmattebon2022triple]*Proposition A.1:

######  Theorem 4 (Le Boudec–Matte Bon).

Let $k\in\mathbb{Z}_{\geq 1}$ and $H\leq\operatorname{Sym}(\Omega)$ be a
$k$-transitive permutation group which satisfies a mixed identity of length
$l$. Then either

$\operatorname{Alt}(\Omega)\trianglelefteq H\quad\mbox{or}\quad k<l.$

######  Remark 3.

Now we come back to the discussion of mixed identities for oligomorphic groups
at the end of the introduction. In view of
[leboudecmattebon2022triple]*Proposition A.1, an interesting case to study is
the automorphism group of the random $k$-uniform hypergraph for $k\geq 3$,
which acts $(k-1)$-transitively but not $k$-transitively in its natural
action. It is an open problem if a dense subgroup could have a highly
transitive action and maybe this is obstructed by the existence of a mixed
identity. This provides an interesting test case for
[hullosin2016transitivity]*Question 6.1.

For the statement of the Hull–Osin-type result in the linear case, we need the
following definition: For a countably infinite dimensional vector space $V$
over $\mathbb{F}_{q}$, we define the normal subgroup $\operatorname{GL}_{\rm
fin}(V)\coloneqq\\{g\in\operatorname{GL}(V)\,|\,\lvert g\rvert<\infty\\}$ of
$\operatorname{GL}(V)$. This group consists of all the elements
$g\in\operatorname{GL}(V)$ such that there is a subspace $W\leq V$ of finite
codimension with $w.g=\lambda w$ for all $w\in W$ and a fixed scalar
$\lambda\in\mathbb{F}_{q}^{\times}$.

Before we state the linear version of Theorem 3, we need the following
preparatory lemma.

######  Lemma 1.

Let $g\in\operatorname{GL}_{\rm fin}(V)$ act as the scalar $\lambda$ on the
subspace $W\leq V$ of finite codimension $n\in\mathbb{N}$ and let $U\leq V$ be
a complement of $W$, i.e. $U\oplus W=V$. Then $U^{\prime}\coloneqq U+U.g$ is a
$g$-invariant subspace of dimension $\leq 2n$. There is a complement
$W^{\prime}\leq W$ of $U^{\prime}$ in $V$.

###### Proof.

Take $u^{\prime}=u+w\in U^{\prime}$ arbitrarily for $u\in U$, $w\in W$. Then
$w=u^{\prime}-u\in U^{\prime}+U=U^{\prime}$, hence
$u^{\prime}.g=u.g+w.g=u.g+\lambda w\in U.g+U^{\prime}=U^{\prime}$, showing
that $U^{\prime}$ is $g$-invariant. Define $W^{\prime}$ to be a complement to
$U^{\prime}\cap W$ in $W$. Then $U^{\prime}\cap W^{\prime}=U^{\prime}\cap
W\cap W^{\prime}=\mathbf{0}$ by definition and
$U^{\prime}+W^{\prime}=U^{\prime}+(U^{\prime}\cap
W)+W^{\prime}=U^{\prime}+W=V$, so $U^{\prime}\oplus W^{\prime}=V$. This
completes the proof. ∎

Here comes the analogue of Theorem 3 for linear groups over finite field.

######  Theorem 5.

Let $V$ be a countably infinite dimensional vector space over $\mathbb{F}_{q}$
and let $H\leq\operatorname{GL}(V)$ be a highly linearly transitive linear
group. Then either

  1. (i)

$H\cap\operatorname{GL}_{\rm fin}(V)$ is a highly transitive, locally finite,
normal subgroup, or

  2. (ii)

$H$ has no mixed identities.

We will first prove a quantitative form, generalizing the result of Le
Boudec–Matte Bon to linear groups.

######  Theorem 6.

Let $V$ be a countably infinite dimensional vector space over a finite field.
Let $k\in\mathbb{Z}_{\geq 1}$ and let $H\leq\operatorname{GL}(V)$ be a
subgroup which is linearly $k$-transitive and admits a mixed identity of
length $l$. Then, at least one of the following cases hold:

  1. (i)

$H\cap\operatorname{GL}_{\rm fin}(V)$ is a locally finite, linearly
$(\left\lceil k/2\right\rceil-1)$-transitive, normal subgroup of $H$, or

  2. (ii)

$k\leq 2l.$

###### Proof.

The inequality in Case (ii) can be proven using Theorem 1 provided all non-
trivial elements of $H$, and hence critical constants in a potential mixed
identity $w\in H\ast\mathbf{F}_{r}$, are of length at least $k/2$. Indeed,
$k>2l$ implies $\lvert w\rvert_{\rm crit}-1\geq\left\lceil
k/2\right\rceil-1\geq l$ and we can apply Theorem 1 with $d\coloneqq 1$ and
$D\coloneqq k$. Hence, any such $w$ of length at most $l$ is not a mixed
identity, contradicting an assumption in the statement of the theorem. We
conclude that $k\leq 2l$ must hold.

Hence, in order to show Case (i), we may suppose that there exists a critical
constant $c\in H$ of non-zero length less than $k/2$. We may assume $k\geq 3.$
Using Lemma 1, we get for some $\lambda\in\mathbb{F}_{q}^{\times}$ that there
exists a vector space $V_{1}$ of codimension $k$ with complement
$V_{0}\cong\mathbb{F}_{q}^{k}$, such that $c$ has a matrix of the form

$c=\begin{pmatrix}c^{\prime}&0\\\ 0&\lambda 1_{V_{1}}\end{pmatrix},$

with respect to the decomposition $V=V_{0}\oplus V_{1}$, where
$c^{\prime}\in\operatorname{GL}_{k}(q)$ is a non-scalar matrix. Now, since $H$
is linearly $k$-transitive, it contains elements of the form

$g=\begin{pmatrix}g^{\prime}&*\\\ 0&g^{\prime\prime}\end{pmatrix}.$

for arbitrary $g^{\prime}\in\operatorname{GL}_{k}(q).$ Thus, since
$\operatorname{SL}_{k}(q)$ is quasi-simple and $c^{\prime}$ is non-scalar, the
normal subgroup of $H$ generated by $c$ will contain elements of the form as
above with $g^{\prime}\in\operatorname{SL}_{k}(q)$ arbitrary. Together with
the $k$-linear transitivity of $H$, this shows that the normal subgroup
generated by $c$ is $k^{\prime}$-linearly transitive for
$k^{\prime}\coloneqq\left\lceil k/2\right\rceil-1$. Indeed, any two linearly
independent sets $\\{\alpha_{1},\dots,\alpha_{k^{\prime}}\\}$ and
$\\{\beta_{1},\dots,\beta_{k^{\prime}}\\}$ can be moved into $V_{0}$ by some
element $h\in H$. Now, since $2k^{\prime}<k$, there exists some element
$g^{\prime}\in\operatorname{SL}(V_{0})\cong\operatorname{SL}_{k}(q)$
exchanging the two sets $\\{\alpha_{1}.h,\dots,\alpha_{k^{\prime}}.h\\}$ and
$\\{\beta_{1}.h,\dots,\beta_{k^{\prime}}.h\\}$. Let $g\in\langle\\!\langle
c\rangle\\!\rangle$, which contains $g^{\prime}$ in its upper left corner.
Then, this shows that $\alpha_{i}.hgh^{-1}=\beta_{i}$ for all
$i\in\\{1,\dots,k^{\prime}\\}.$ It is clear that the normal subgroup generated
by $c$ is contained in $H\cap\operatorname{GL}_{\rm fin}(V)$ and that this
group is locally finite. This finishes the proof. ∎

###### Proof of Theorem 5.

Since $H$ is highly linearly transitive, we can apply Theorem 6 for any
$k\in\mathbb{Z}_{\geq 1}$. In particular, the existence of a mixed identity
implies Case (i) in Theorem 6 for arbitrary $k$. Hence,
$H\cap\operatorname{GL}_{\rm fin}(V)$ is a highly transitive, locally finite,
normal subgroup of $H$.

However, the statement of Theorem 5 is a bit stronger since it says that Cases
(i) and (ii) are mutually exclusive. We will now show that Case (i) arises
whenever $H$ contains an element $h$ such that $\lvert h\rvert<\infty$, i.e.
whenever $H\cap\operatorname{GL}_{\rm fin}(V)$ is non-trivial. Let $W\leq V$
be a subspace of codimension $n$ such that $h$ acts as multiplication by the
fixed scalar $\lambda\in\mathbb{F}_{q}^{\times}$ on $W$. For any
$g\in\operatorname{GL}(V)$, the element $h^{g}\in\operatorname{GL}_{\rm
fin}(V)$ acts as $\lambda\operatorname{id}_{W.g}$ on $W.g$. Hence $[h,g]$ acts
as the identity on $W_{0}\coloneqq W\cap W.g$. Then
$\operatorname{codim}(W_{0})\leq 2n$. By Lemma 1, we find $W^{\prime}\leq
W_{0}$ and $U^{\prime}$ a $[h,g]$-invariant complement of $W^{\prime}$, such
that $\dim(U^{\prime})=\operatorname{codim}(W^{\prime})\leq 4n$. Then $[h,g]$
can be seen as an element of
$\operatorname{GL}(U^{\prime})\times\\{1_{W^{\prime}}\\}\leq\operatorname{GL}_{4n}(q)\times\mathbf{1}$,
so $[h,x]^{e}\in H\ast\langle x\rangle=H\ast\mathbf{F}_{1}$ is a mixed
identity for $H$, where $e$ denotes the exponent of
$\operatorname{GL}_{4n}(q)$. ∎

##  5\. No mixed identities for $A(q)$

It was proved in [bradfordschneiderthom2023nonsol]*Proof of Theorem 1.10 that
the topological full group of a faithful and minimal, continuous action of a
group on a Cantor set is mixed identity free. This is essentially a direct
application of the Hull–Osin theorem to the action of the full group on an
infinite orbit. A similar reasoning can be applied to the study of full groups
of an ergodic, probability measure preserving equivalence relation on a non-
atomic standard probability space. In this section, we want to prove a linear
version of this, where a certain group $A(q)$, constructed in
[carderithom2018exotic], replaces the full group.

Let’s briefly recall the construction of $A(q)$, which resembles a
construction of the full group of the hyperfinite equivalence relation. Let
$q$ be a prime power as usual. For $n\in\mathbb{N}$, consider
$\operatorname{GL}_{2^{n}}(q)$ and the diagonal embeddings

$\iota_{n}\colon\operatorname{GL}_{2^{n}}(q)\to\operatorname{GL}_{2^{n+1}}(q);h\mapsto\begin{pmatrix}h&0\\\
0&h\end{pmatrix}.$

Write
$\operatorname{GL}_{2^{\omega}}(q)=\bigcup_{n\in\mathbb{N}}\operatorname{GL}_{2^{n}}(q)$
for the inductive limit of the countable chain of inclusions
$(\operatorname{GL}_{2^{n}}(q),\iota_{n})_{n\in\mathbb{N}}$. Clearly, this is
a countable, locally finite group. Note that the normalized rank metric is
well-defined on $\operatorname{GL}_{2^{\omega}}(q)$. We denote its Cauchy
completion by $A(q)$. Note that we could have started with
$\operatorname{SL}_{2^{n}}(q)$ and arrived at the same completion.
Furthermore, $A(q)/Z$, where $Z$ denotes its center
$Z\cong\mathbb{F}_{q}^{\times}$, is the completion of the inductive limit of
the groups $\operatorname{PGL}_{2^{n}}(q)$ with respect to the normalized
projective rank metric. See [carderithom2018exotic] for more details on the
group $A(q).$ Throughout this section, we work with the normalized projective
rank metric and a suitable version of Theorem 1 and its corollaries.

Now we prove the main result of this section:

######  Theorem 7.

The group $A(q)/Z$ does not admit a mixed identity.

###### Proof.

Let

$w=c_{0}x^{\varepsilon(1)}_{\iota(1)}c_{1}\cdots
c_{l-1}x^{\varepsilon(l)}_{\iota(l)}c_{l}\in A(q)\ast\mathbf{F}_{r}$

be some reduced word with constants, where $l\in\mathbb{Z}_{\geq 2}$. We prove
that $w$ has a non-trivial word image $w(A(q)^{r})$. At first, note that, by
construction, we can find a large number $n=n(\varepsilon)\in\mathbb{N}$, such
that we can choose constants
$c^{\prime}_{0},\dots,c^{\prime}_{l}\in\operatorname{SL}_{2^{n}}(q)$ so that
$d(c^{\prime}_{j},c_{j})\leq\varepsilon$ for all $i\in\\{0,\ldots,l\\}$.
Define the word

$w^{\prime}\coloneqq
c^{\prime}_{0}x^{\varepsilon(1)}_{\iota(1)}c^{\prime}_{1}\cdots
c^{\prime}_{l-1}x^{\varepsilon(l)}_{\iota(l)}c^{\prime}_{l}\in\operatorname{SL}_{2^{n}}(q)\ast\mathbf{F}_{r}<A(q)\ast\mathbf{F}_{r}.$

Note that then for all $h_{1},\ldots,h_{r}\in A(q)$, when have

$\displaystyle d(w^{\prime}(h_{1},\dots,h_{r}),$ $\displaystyle
w(h_{1},\dots,h_{r}))$
$\displaystyle=d(c^{\prime}_{0}h^{\varepsilon(1)}_{\iota(1)}c^{\prime}_{1}\cdots
c^{\prime}_{l-1}h^{\varepsilon(l)}_{\iota(l)}c^{\prime}_{l},c_{0}h^{\varepsilon(1)}_{\iota(1)}c_{1}\cdots
c_{l-1}h^{\varepsilon(l)}_{\iota(l)}c_{l})$ $\displaystyle\leq
d(c^{\prime}_{0}h^{\varepsilon(1)}_{\iota(1)}c^{\prime}_{1}\cdots
c^{\prime}_{l-1}h^{\varepsilon(l)}_{\iota(l)}c^{\prime}_{l},c^{\prime}_{0}h^{\varepsilon(1)}_{\iota(1)}c^{\prime}_{1}\cdots
c^{\prime}_{l-1}h^{\varepsilon(l)}_{\iota(l)}c_{l})$
$\displaystyle+\cdots+d(c^{\prime}_{0}h^{\varepsilon(1)}_{\iota(1)}c_{1}\cdots
c_{l-1}h^{\varepsilon(l)}_{\iota(l)}c_{l},c_{0}h^{\varepsilon(1)}_{\iota(1)}c_{1}\cdots
c_{l-1}h^{\varepsilon(l)}_{\iota(l)}c_{l})$
$\displaystyle\leq\sum_{j=0}^{l}{d(c^{\prime}_{j},c_{j})}\leq(l+1)\varepsilon,$

Hence, we obtain

$\lvert\operatorname{diam}(w^{\prime}(A(q)^{r}))-\operatorname{diam}(w(A(q)^{r}))\rvert\leq
2(l+1)\varepsilon.$

From the definition $\lvert w^{\prime}\rvert_{\rm crit}\geq\lvert w\rvert_{\rm
crit}-\varepsilon$ and for $m\geq n$, Corollary 1 yields

$\operatorname{diam}(w^{\prime}(\operatorname{SL}_{2^{m}}(q)^{r}))\geq\frac{\lvert
w^{\prime}\rvert_{\rm crit}}{l}-\frac{1}{2^{m}}\geq\frac{\lvert w\rvert_{\rm
crit}-\varepsilon}{l}-\frac{1}{2^{m}}$

and hence

$\operatorname{diam}(w(\operatorname{SL}_{2^{m}}(q)^{r}))\geq\frac{\lvert
w\rvert_{\rm crit}-\varepsilon}{l}-\frac{1}{2^{m}}-2(l+1)\varepsilon.$

Thus, taking $\varepsilon>0$ small enough and $m\geq n(\varepsilon)$ large
enough, we see that the word image of $w$ is non-trivial. This finishes the
proof. ∎

## Acknowledgments

We are grateful to Melissa Lee for an enlightening discussion concerning
Remark 1. The third author wants to thank the Mathematisches
Forschungsinstitut in Oberwolfach for its hospitality.

## References