11institutetext: Daria Fokina 22institutetext: Fraunhofer ITWM,
Kaiserslautern, Germany, Technical University Kaiserslautern, Kaiserslautern,
Germany, 22email<EMAIL_ADDRESS>33institutetext: Pavel
Toktaliev 44institutetext: Fraunhofer ITWM, Kaiserslautern, Germany, Technical
University Kaiserslautern, Kaiserslautern, Germany, 44email:
<EMAIL_ADDRESS>55institutetext: Oleg Iliev
66institutetext: Fraunhofer ITWM, Kaiserslautern, Germany, Bulgarian Academy
of Sciences, Sofia, Bulgaria, 66email<EMAIL_ADDRESS>77institutetext: Ivan Oseledets 88institutetext: Skolkovo Institute of Science
and Technology, Moscow, Russia, Artificial Intelligence Research Institute,
Moscow, Russia, 88email<EMAIL_ADDRESS>
# Machine learning methods for prediction of breakthrough curves in reactive
porous media

Daria Fokina  Pavel Toktaliev  Oleg Iliev and Ivan Oseledets

###### Abstract

Reactive flows in porous media play an important role in our life and are
crucial for many industrial, environmental and biomedical applications. Very
often the concentration of the species at the inlet is known, and the so-
called breakthrough curves, measured at the outlet, are the quantities which
could be measured or computed numerically. The measurements and the
simulations could be time-consuming and expensive, and machine learning and
Big Data approaches can help to predict breakthrough curves at lower costs.
Machine learning (ML) methods, such as Gaussian processes and fully-connected
neural networks, and a tensor method, cross approximation, are well suited for
predicting breakthrough curves. In this paper, we demonstrate their
performance in the case of pore scale reactive flow in catalytic filters.

##  1 Introduction

Reactive flows in porous media play an important role in our life and are
crucial for many industrial, environmental and biomedical applications.
Therefore, it is important to understand how the reactive transport works and
how one can control or improve it. Note that chemical reactions occur at pore
scale, however the measurements are usually done at macroscopic scale. Very
often the concentration of the species at the inlet is known, and the so-
called breakthrough curves are the quantities which could be measured or
computed numerically. Breakthrough curves represent the concentration of the
species at the flow outlet. In problems like parameter identification,
optimizations, etc. breakthrough curves for many of the input parameters. The
measurements and simulations could be time-consuming and expensive, and
machine learning and Big Data approaches can help to predict breakthrough
curves at lower costs. Machine learning (ML) methods, such as Gaussian
processes and fully-connected neural networks, and a tensor method, cross
approximation, are well suited for predicting breakthrough curves. In this
paper, we demonstrate their performance in the case of pore scale reactive
flow in catalytic filters. The computational geometry for a piece of the
filtering membrane could come from 3D CT image, or can be generated by proper
software. The governing equation is the convection-diffusion-reaction
equation, and linear and nonlinear reactions are considered. The computational
experiments are followed by discussions about the required size of the
training set, applicability and the performance of each of the methods for
different flow regimes, namely diffusion dominated and reaction dominated. We
expect that the considered methods will show similar performance for
predicting breakthrough curves for other industrial or environmental problems,
supposing similar reactions and flow regimes are considered.

##  2 Problem formulation

We consider at microscale a reactive flow through a piece of a catalytic
filter membrane. A sample of the geometry can be seen on Fig. 1(a). The red
color is reserved for the solid inert particles, the green color denotes the
active washcoat particles, the microscale pores are transparent. There is no
transport through the inert particles. Reactions occur within the washcoat ,
they occupy the domain indicated by $\Omega_{w}$. The resolved microscale
pores occupy a domain denoted as $\Omega_{f}$. Thus, the computational domain
is $\Omega=\Omega_{w}\cup\Omega_{f}$. Mathematically, the process is modeled
by the the following transport equation:

$c_{t}-D\;\Delta c+(\vec{v},\nabla c)+r(c)=0,\quad x\in\Omega,t>0,$ (1)

$D$ is the diffusion coefficient which is different in the pores and in the
washcoat , $v$ is the stationary velocity field which is computed in advance,
see below, $r(c)=0,\,x\in\Omega_{f}$. For the reaction in the washcoat, we
consider linear and nonlinear reactions. Further in the text we refer only to
the processes in the washcoat. If the reaction is linear: $r(c)=kc$, then in a
dimensionless form the equation can be written as follows:

$c_{t}-\mu_{\mathrm{Da}}\;\Delta c+(\vec{v},\nabla c)+\mu_{\mathrm{Pe}}\,c=0,$
(2)

where $\mu_{\mathrm{Da}}=\frac{kL}{D}$ is Damköhler number,
$\mu_{\mathrm{Pe}}=\frac{u_{in}L}{D}$ is Peclet number, $L$ is length of the
computational domain, $u_{in}$ is the velocity at the inlet.

The initial and boundary conditions are specified as follows:

$\begin{gathered}c(x,0)=0,\quad x\in\Omega,\\\ c(x,t)=1,\quad
x\in\Omega_{\mathrm{inlet}},t>0,\\\ \frac{\partial
c}{\partial\vec{n}}(x,t)=0,\quad
x\in\Omega_{\mathrm{outlet}}\cup\Omega_{\mathrm{walls}}.\\\ \end{gathered}$
(3)

Further we use notation: $\vec{\mu}=(\mu_{\mathrm{Da}},\mu_{\mathrm{Pe}})$.

The quantity of interest is the time dependent concentration at the outlet,
i.e. the breakthrough curve:

$s(t;\mu)=s_{\mu}(t)=\int\limits_{\Omega_{\mathrm{outlet}}}c(x,t;\mu)dx.$ (4)

(a)

(b)

Figure 1: Considered geometry (a) and velocity field inside the media (b).

We use _supervised_ machine learning methods to approximate the dependency
$s_{\mu}(t)$ as a function of the input parameters $\mu$. For that we need the
values $s_{\mu}$ for a set of values of $\mu$. This can be done experimentally
or computed numerically. In both cases we do not get a function of time
$s_{\mu}(t)$, but vectors $\vec{s}_{\mu}$ for fixed moments of time. To obtain
$\vec{s}_{\mu}$ we use PoreChem software porechem , which provides us with the
results of numerical simulations. Based on these numerical values we train a
machine learning model and then choose the best configuration based on the
results on new data, validation data, yet not seen by the model. The third
set, the test set, is used to present the final results of the best models.

##  3 Considered methods

In the scope of this paper, we compare two machine learning methods: Gaussian
processes and fully connected neural networks and one tensor method - cross
approximation. In the following, we discuss these methods in detail. Please
note that for our problem other machine learning methods could be considered
as well.

###  3.1 Machine learning methods

The considered machine learning methods considered are the so-called
_regression_ methods. This means that we have a function
$f:\Omega\subset\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, which we would like
to approximate from several data points: $\mathcal{D}=(X,Y)$ which consists of
pairs $\\{(\vec{x}_{i},\vec{y}_{i}=f(x_{i}))\\},i=1,\ldots,N$, $x_{i}$ is the
input, and $y_{i}$ is the target value. The approximation is represented as a
parameterized mapping $\hat{f}_{\theta}$, $\theta$ are the parameters.
Normally optimal parameters $\theta^{*}$ that give the best fit to the data
points are chosen by minimizing mean squared error (MSE):

$l(\theta;\mathcal{D})=\sum\limits_{(\vec{x}_{i},\vec{y}_{i})\in\mathcal{D}}\|\hat{f}_{\theta}(\vec{x}_{i})-\vec{y}_{i}\|^{2}\rightarrow\min\limits_{\theta}.$
(5)

This functional is often used in neural networks.

#### Gaussian processes

The first method to consider is a Gaussian process. For simplicity of the
description, consider a function $f:\Omega\rightarrow\mathbb{R}$. Then set of
given target values $Y=\vec{y}$ is a vector. This is a Bayesian method, which
gives the conditional distribution $p(y_{*}|x_{*},X,Y)$ of $y_{*}=f(x_{*})$ at
a new point $x_{*}$, given a training dataset $(X,Y)$. Consider a Gaussian
process $f(x)\sim GP(0,K)$, where $K$ is a kernel function, or simply a
kernel. This means that for any vector $\vec{y}$,
$y_{i}=f(x_{i}),\,i=1,\ldots,N$, has a normal distribution
$\mathcal{N}(\vec{0},K(X,X)),\,K(X,X)_{ij}=K(x_{i},x_{j}),i,j=1,\ldots,N$.
When several values of $f$ are known for the set of points $X\subset\Omega$,
the posterior distribution of $y_{*}=f(x_{*})$ for the new point
$x_{*}\in\Omega$ is the following:

$\begin{split}{y}_{*}|\vec{y},X,x_{*}\sim\mathcal{N}(&K(x_{*},X)K(X,X)^{-1}\vec{y},\\\
&K(x_{*},x_{*})-K(x_{*},X)K(X,X)^{-1}K(X,x_{*})).\end{split}$ (6)

The approximation is the mean of this distribution, i.e.
$\hat{f}(x_{*})=K(x_{*},X)K(X,X)\vec{y}$.

We consider three kernel functions:

  * •

Radial Basis Function (RBF) kernel:
$K(x,z)=\exp\left(-\frac{d(x,z)^{2}}{2l^{2}}\right)$,

  * •

Matern kernel:
$K(x,z)=\frac{1}{\Gamma(\nu)2^{\nu-1}}\left(\frac{\sqrt{2\nu}d(x,z)}{l}\right)^{l}K_{\nu}\left(\frac{\sqrt{2\nu}d(x,z)}{l}\right)$,

  * •

Rational quadratic (RQ) kernel: $K(x,z)=\left(\frac{1+d^{2}(x,z)}{2\alpha
l^{2}}\right)^{-\alpha}$,

where $K_{\nu}$ is the modified Bessel function, $\Gamma(\nu)$ is the Gamma
function, $d(x,z)$ is the Euclidean distance function.

#### Neural networks

Another class of ML models which has demonstrated a good efficiency are
artificial neural networks (NNs). In our considerations we use a classical
fully-connected neural network. For this kind of neural network, the
approximation is the following:

$f_{\theta}=L_{d}\circ g\circ\ldots\circ L_{2}\circ g\circ L_{1},$ (7)

where
$L_{k}(\vec{x})=\vec{x}W_{k}+\vec{b}_{k},\,W_{k}\in\mathbb{R}^{h_{k-1}\times
h_{k}},\vec{b}_{k}\in\mathbb{R}^{h_{k}}$,
$\theta=\\{(W_{k},\vec{b}_{k}\\}_{k=1}^{d}$, $h_{0}=n,\,h_{d}=m$, $g$ is a
nonlinear function, which is also called an activation function. The
parameters are found by solving the minimization problem 5. This is usually
done by stochastic gradient descent or its modifications. In this work we use
Adam optimizer adam . We experiment with different activation functions:

  * •

Rectified linear unit (ReLU): $g(x)=\max(0,x)$,

  * •

Leaky ReLU: $g(x)=\begin{cases}x,\;\text{if }x>0\\\ \alpha x,\;\text{if }x\leq
0\end{cases},\;\alpha=0.01$,

  * •

Exponential linear unit (ELU): $g(x)=\begin{cases}x,\;\text{if}\,x>0\\\
\alpha\exp(x-1),\;ifx\leq 0\end{cases},\;\alpha=0.01$.

###  3.2 Cross approximation

#### Discrete case

Unlike the previous methods, cross approximation is a tensor method. As it
follows from the name, it operates on tensors, i.e. on $n$-dimensional
matrices. We describe this tensor method, and then explain extension to
continuous case.

The input of the model is a vector
$\mathbf{\mu}=(\mu_{\mathrm{Da}},\mu_{\mathrm{Pe}})$, each component of which
is located within a closed interval of interest:
$\mu_{\mathrm{Da}}\in\gamma_{1}=[p_{1}^{1},p_{2}^{1}],\,\mu_{\mathrm{Pe}}\in\gamma_{1}=[p_{1}^{2},p_{2}^{2}]$.
At the first step, the space of input parameters $\gamma_{1}\times\gamma_{2}$
is restricted to a $k\times k$ grid $A\times B$. This is done by uniformly
discretizing each of intervals $\gamma_{1}$ and $\gamma_{2}$. The notation is
the following:

$\mu_{\mathrm{Da}}\in\\{A_{i}\\}$, |  $A_{i}\in\gamma_{1}$, | $A_{1}=p_{1}^{1},A_{k}=p_{2}^{1}$
---|---|---
$\mu_{\mathrm{Pe}}\in\\{B_{i}\\}$, |  $B_{i}\in\gamma_{2}$, | $B_{1}=p_{1}^{2},A_{k}=p_{2}^{2}$

If we compute the curves $s_{\mu}$ for each element on the grid, we get a
three-dimensional tensor $\mathsf{S}$$\in\mathbb{R}^{k\times k\times N_{T}}$,
where $N_{T}$ is the number of simulated time steps. We suppose that
$\mathsf{S}$ can be represented in a low-rank form using Tucker decomposition:

$\mathsf{S}=\mathsf{G}\times_{1}U\times_{2}V,$ (8)

$G\in\mathbb{R}^{r_{1}\times r_{2}\times N_{T}},\,U\in\mathbb{R}^{r_{1}\times
k},\,V\in\mathbb{R}^{r_{2}\times k}$. The multiplication $\times_{m}$ is
performed along the axis $m$. For $\mathsf{C}=\mathsf{A}\times_{m}B$, where
$\mathsf{A}\in\mathbb{R}^{d_{1},\ldots,d_{n}},B\in\mathbb{R}^{d_{m},l}$, its
elements are expressed as:

$C_{i_{1},\ldots,i_{m-1},j,i_{m+1},\ldots,i_{n}}=\sum\limits_{i_{m}}A_{i_{1},\ldots,i_{n}}B_{i_{m},j}.$
(9)

Cross-approximation provides us with an algorithm to choose $U,\,V$ and $G$.
Denote slices of $\mathsf{S}$ for some $k_{1}$ and $k_{2}$ as
$P_{1}:=S[k_{1},:.:]$, $P_{2}:=S[:,k_{2},:]$. For this matrices we compute a
singular value decomposition:

$P_{1}=\hat{U}S_{1}\Phi_{1};\,P_{2}=\hat{V}S_{2}\Phi_{2},$ (10)

with $U^{T}=\hat{U}[:,:r_{1}]\in\mathbb{R}^{k\times r_{1}}$ and
$V^{T}=\hat{V}[:,:r_{2}]\in\mathbb{R}^{k\times r_{2}}$ are the leftmost
singular vectors $r_{1}$ and $r_{2}$ respectively. To find $\mathsf{G}$ a
MaxVolmaxvol  algorithm is used. MaxVol is a greedy algorithm looking in a
matrix for a submatrix of maximum volume (for a square matrix
$\text{vol}(A)=|\text{det}(A)|$). This algorithm gives us a set of indices,
for which we need to compute values of $\mathsf{S}$, to get the tensor
$\mathsf{G}$, and, respectively, to reconstruct the whole tensor $\mathsf{S}$.
The elements necessary for computing $S$ are indicated in red and blue on Fig.
2

[width=.6]ca_tensor

Figure 2: Elements of tensor, required for cross approximation.

#### Continuous case

We consider the tensor $\mathsf{S}$ in the Tucker decomposition; this means
that its elements can be represented as:

$S_{mnt}=\sum\limits_{i,j=1}^{r_{1},r_{2}}G_{ijt}U_{im}V_{jn},\quad
m=1,\ldots,k,\,n=1,\ldots,k,\,t=1,\ldots,N_{T}.$ (11)

The value of $u_{im}$ depends only on the value $A_{i}$ of
$\mu_{\mathrm{Da}}$, the value of $v_{jn}$ depends only on the value $B_{j}$
of $\mu_{\mathrm{Pe}}$. In the continuous case, it would be written as:

$S(\mu_{\mathrm{Da}},\mu_{\mathrm{Pe}},t)=\sum\limits_{ij}G_{ij}(t)u_{i}(\mu_{\mathrm{Da}})v_{j}(\mu_{\mathrm{Pe}}).$
(12)

We call $u_{i}(\mu_{\mathrm{Da}})$ and $v_{j}(\mu_{\mathrm{Pe}})$ as basis
functions and find them by interpolating the vectors $U_{i}$ and $V_{j}$
correspondingly.

##  4 Computational experiments

We conduct the experiments on a 3D geometry - a random packing of spheres. The
size of it is $50\times 50\times 50$ voxels. The length of the sample is
$L=0.0001$ m, the diameter of the spheres is $l=6$ mkm, the velocity on the
inlet is $u_{\mathrm{in}}=0.047$, porosity is $\phi=0.741$. The geometry is
shown on Fig. 1(a) with the precomputed velocity field on Fig. 1(b). To
calculate the velocities, Geodict software tool geodict  was used. For each
value of input parameters we compute concentration for 200 time steps
($N_{T}$) for $t\in[0.00001,0.002]$.

We start with different regimes for the case of linear reaction and then
consider the case of a strong nonlinear reaction. For the training of Gaussian
processes and neural networks, we uniformly randomly sample the input
parameters from the domain of interest. In each of the cases, the size of the
training set for Gaussian processes is equal to 100 and for neural networks it
is 600. For the cross approximation method we need parameter values for the
middle cross-section of the parameter domain and also at positions found by
MaxVol algorithm. For the evaluation of these three methods, we randomly
sample 100 points for validation and separately 100 points for the test set.
For all of these generated parameters, we numerically compute the breakthrough
curves (target values for the models). In order to compare the results, we
compare two errors:

$\displaystyle e_{l_{2}}$
$\displaystyle=\sqrt{\frac{\sum\limits_{x\in\hat{\Omega}}\sum\limits_{i=1}^{m}(\hat{s}_{i}(x)-s_{i}(x))^{2}}{\sum\limits_{x\in\hat{\Omega}}\sum\limits_{i=1}^{m}s_{i}(x)^{2}}}$
(13) $\displaystyle e_{\max}$
$\displaystyle=\max\limits_{x\in\hat{\Omega}}\sqrt{\frac{\sum\limits_{i=1}^{m}(\hat{s}_{i}(x)-s_{i}(x))^{2}}{\sum\limits_{i=1}^{m}s_{i}(x)^{2}}}$
(14)

###  4.1 Linear reaction

As it was mentioned above, we begin with a linear reaction $r(c)=k\,c$. There
are two regimes considered:

  * •

Diffusion dominated with
$\mu_{\textrm{Da}}\in[0.01,1.1],\,\mu_{\textrm{Pe}}\in[0.01,1.1]$,

  * •

Strong reaction with
$\mu_{\textrm{Da}}\in[800.,1200.],\,\mu_{\textrm{Pe}}\in[0.05,0.1]$.

#### Diffusion dominated case

The results for Gaussian processes for different kernels are provided in Table
4.1, for neural networks for different activation functions are in Table 4.1,
for cross approximation — in Table 4.1. Table 4.1 provides the test errors for
each method. Fig. 3 compares breakthrough curves predicted by the models
(solid line) with numerically computed target values (dotted line) for several
values of input parameters $(\mu_{\textrm{Da}},\mu_{\textrm{Pe}})$. The
predicted values are provided for the best configuration of each method. The
difference between predicted and target values is not seen by a human eye. On
a zoomed-in figure a difference can be noticed, though it is also very small.

As one can see, the best result is achieved by a neural network with ELU
activation function. The cross approximation gives a similar result, but with
a smaller amount of training data.

Table 1: Linear reaction, diffusion dominated case. Results for Gaussian
processes for different kernels.

Kernel | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
RBF |  $0.00161$ |  $0.00489$ |  $0.04281$ |  $0.22586$
Matern |  $1.0\times 10^{-7}$ |  $5.96\times 10^{-7}$ |  $\mathbf{0.04777}$ |  $\mathbf{0.20501}$
RQ |  $1.67\times 10^{-6}$ |  $8.71\times 10^{-6}$ |  $0.06546$ |  $0.33635$

Table 2: Linear reaction, diffusion dominated case. Results for neural
networks for activation functions.

Activation function | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
ReLU |  $0.00011$ |  $0.00040$ |  $0.00166$ |  $0.00905$
Leaky ReLU |  $0.00038$ |  $0.00136$ |  $0.00248$ |  $0.01771$
ELU |  $0.00084$ |  $0.00475$ |  $\mathbf{0.00113}$ |  $\mathbf{0.00428}$

Table 3: Linear reaction, diffusion dominated case. Results for cross-
approximation and different interpolation methods.

Interpolation | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
Polynomial |  $0.01026$ |  $0.03127$ |  $0.01041$ |  $0.03135$
Piecewise linear |  $0.00210$ |  $0.01820$ |  $\mathbf{0.00380}$ |  $\mathbf{0.01974}$

(a)

(b)

(c)

Figure 3: Linear reaction, diffusion dominated case. Breakthrough curve
samples for target values and values predicted by: (a) Gaussian process, (b)
neural network, (c) cross approximation.

Table 4: Linear reaction, diffusion dominated case. Test errors

Method | Test error
---|---
$e_{l_{2}}$ |  $e_{\max}$
Gaussian process |  $0.01716$ |  $0.13296$
Neural network |  $0.00055$ |  $0.00126$
Cross approximation |  $0.00046$ |  $0.00219$

#### Strong reaction

In a case of a strong linear reaction (large Da number) the results are
presented in Table 4.1 for Gaussian processes, Table 4.1 for neural networks
and Table 4.1 for cross approximation. In Table 4.1 test errors are presented.
We also plot the breakthrough curve samples, numerically computed and
predicted by the considered methods.

In this case, the best performance shows cross-approximation + polynomial
interpolation. For Gaussian processes with rational quadratic kernel the error
is bigger, however the difference to cross approximation is minor. Neural
networks have the worst result, although they use the most of the data. On the
Fig. 4 we can notice the difference between the values predicted by the model
and and numerically simulated data for neural networks, but this difference is
noticed only at a closer look at the figure.

Table 5: Linear strong reaction. Results for Gaussian processes for different
kernels.

Kernel | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
RBF |  $6.79\times 10^{-6}$ |  $3.57\times 10^{-5}$ |  $9.32\times 10^{-6}$ |  $4.04\times 10^{-5}$
Matern |  $1.83\times 10^{-5}$ |  $0.00017$ |  $3.97\times 10^{-5}$ |  $0.00028$
RQ |  $2.05\times 10^{-6}$ |  $9.41\times 10^{-6}$ |  $\mathbf{8.37\times 10^{-6}}$ |  $\mathbf{7.35\times 10^{-5}}$

Table 6: Linear strong reaction. Results for neural networks for activation
functions.

Activation function | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
ReLU |  $0.00297$ |  $0.00661$ |  $0.00322$ |  $0.00573$
Leaky ReLU |  $6.14\times 10^{-5}$ |  $0.00016$ |  $\mathbf{0.00055}$ |  $\mathbf{0.00011}$
ELU |  $0.00109$ |  $0.00234$ |  $0.001163$ |  $0.00186$

Table 7: Linear strong reaction. Results for cross-approximation and different
interpolation methods.

Interpolation | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
Polynomial |  $2.24\times 10^{-6}$ |  $4.24\times 10^{-6}$ |  $\mathbf{2.20\times 10^{-6}}$ |  $\mathbf{3.48\times 10^{-6}}$
Piecewise linear |  $4.48\times 10^{-6}$ |  $9.92\times 10^{-6}$ |  $4.32\times 10^{-6}$ |  $8.86\times 10^{-6}$

Table 8: Linear strong reaction. Test errors

Kernel | Test error
---|---
$e_{l_{2}}$ |  $e_{\max}$
Gaussian process |  $8.3\times 10^{-6}$ |  $4.13\times 10^{-5}$
Neural network |  $8.83\times 10^{-5}$ |  $0.00027$
Cross approximation |  $2.11\times 10^{-6}$ |  $4.03\times 10^{-6}$

(a)

(b)

(c)

Figure 4: Linear strong reaction. Breakthrough curve samples for target values
and values predicted by: (a) Gaussian process, (b) neural network, (c) cross
approximation.

###  4.2 Nonlinear reaction

Here we consider a nonlinear reaction

$r(c)=-\frac{k\cdot c}{(1+k_{\text{inhib}}c)^{2}},$ (16)

$k_{\text{inhib}}=4m^{3}/mol$.

Like in the case of the the linear reaction, we introduce here in the same way
$Da$ and $Pe$ numbers. The parameter intervals for this case are the same as
for the strong linear reaction: $Da\in[800.,1200.],\,Pe\in[0.05,0.1]$.

We present errors of the models in Tables 4.2, 4.2, 4.2 for Gaussian
processes, neural networks and cross approximation correspondingly. Test
errors can be found in Table 4.2. In this case we plot also samples of
predicted and precomputed breakthrough curves as well (Fig. 5). The results of
the methods are similar to the results in the linear case. Here it should be
noted that breakthrough curves in both linear and nonlinear cases are quite
similar and can sometimes coincide. This depends on the parameter values, i.e.
ratio $Da/Pe$.

(a)

(b)

(c)

Figure 5: Nonlinear reaction. Breakthrough curve samples for target values and
values predicted by: (a) Gaussian process, (b) neural network, (c) cross
approximation. Table 9: Nonlinear reaction. Results for Gaussian processes for
different kernels.

Kernel | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
RBF |  $9.43\times 10^{-6}$ |  $7.24\times 10^{-5}$ |  $\mathbf{7.20\times 10^{-6}}$ |  $\mathbf{2.76\times 10^{-5}}$
Matern |  $6\times 10^{-5}$ |  $0.00075$ |  $2.15\times 10^{-5}$ |  $0.00013$
RQ |  $1.12\times 10^{-5}$ |  $0.00013$ |  $7.91\times 10^{-6}$ |  $6.42\times 10^{-5}$

Table 10: Nonlinear reaction. Results for neural networks for activation
functions.

Activation function | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
ReLU |  $6.91\times 10^{-5}$ |  $0.00031$ |  $8.48\times 10^{-5}$ |  $0.00031$
Leaky ReLU |  $3.77\times 10^{-5}$ |  $8.87\times 10^{-5}$ |  $\mathbf{5.35\times 10^{-5}}$ |  $\mathbf{0.00021}$
ELU |  $0.00045$ |  $0.00055$ |  $0.00045$ |  $0.00054$

Table 11: Nonlinear reaction. Results for cross-approximation and different
interpolation methods.

Interpolation | Train error | Validation error
---|---|---
$e_{l_{2}}$ |  $e_{\max}$ |  $e_{l_{2}}$ |  $e_{\max}$
Polynomial |  $1.89\times 10^{-6}$ |  $3.64\times 10^{-6}$ |  $\mathbf{1.82\times 10^{-6}}$ |  $\mathbf{3.59\times 10^{-6}}$
Piecewise linear |  $4.06\times 10^{-6}$ |  $9.78\times 10^{-6}$ |  $4.05\times 10^{-6}$ |  $9.92\times 10^{-6}$

Table 12: Nonlinear reaction. Test errors

Kernel | Test error
---|---
$e_{l_{2}}$ |  $e_{\max}$
Gaussian process |  $6.55\times 10^{-6}$ |  $2.15\times 10^{-5}$
Neural network |  $4.68\times 10^{-5}$ |  $0.00016$
Cross approximation |  $1.83\times 10^{-6}$ |  $3.44\times 10^{-6}$

##  5 Conclusion

We have performed a comparison study of two machine learning methods and one
tensor method extended by interpolation when applied to predict breakthrough
curves for a catalytic filter. The training sets were obtained by simulating
at pore scale 3D reactive flow using the software tool PoreChem. All the
methods show good performance, though for different regimes the order of error
differs. In the the case of a strong reaction the errors obtained are much
smaller than in the diffusion-dominated case. This may be due to the fact,
that the curves changed less for the chosen intervals of parameters. The
methods requiring small amount of data (Gaussian processes and cross
approximation) are more beneficial in this case, as numerical simulations for
larger geometries require hours of running to compute just one point. In our
case, it took several minutes to compute one sample.

In the case of cross approximation we have empirically shown, that the
selected basis vectors for the low-rank approximation together with the data
selected by maximum volume algorithm, provide a good approximation to the
considered tensor. Further interpolation of the basis vectors gives us a good
extension to the continuous case.

Future work includes applying the considered methods for cases of more complex
reactions and more complex geometries. The investigation how the intervals
influences the models performances is a matter of further studies.
Additionally, these approaches can be considered for similar tasks in other
areas, e.g. decline curves in oil industry.

## References

  * [1] J. Becker, F. Biebl, E. Glatt, L. Cheng, A. Grießer, M. Groß, S. Linden, D. Mosbach, C. Wagner, A. Weber, and R. Westerteiger.  Geodict (release 2022) [simulation software].  doi.org/10.30423/release.geodict2022, 2021.
  * [2] S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov, and N. L. Zamarashkin.  How to find a good submatrix.  In Vadim Olshevsky and Eugene Tyrtyshnikov, editors, Matrix Methods: Theory, Algorithms and Applications, pages 247–256. WORLD SCIENTIFIC, 2010.
  * [3] O. Iliev and P. Toktaliev.  On pore scale numerical simulation of complex homogeneous reactions with application to filtration processes.  FILTECH, 2022.
  * [4] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.


* Eriksson & Rizzi (1985) Eriksson, L. E. & Rizzi, A. 1985 Computer-aided analysis of the convergence to steady state of discrete approximations to the Euler equations. Journal of Computational Physics 57 (1), 90–128.
  * Falgout & Yang (2002) Falgout, R. D. & Yang, U. M. 2002 hypre: A library of high performance preconditioners. In International Conference on Computational Science, pp. 632–641.
  * Gómez et al. (2016) Gómez, F., Sharma, A. S. & Blackburn, H. M. 2016 Estimation of unsteady aerodynamic forces using pointwise velocity data. Journal of Fluid Mechanics 804, R4.
  * Hairer et al. (1993) Hairer, E., Nørsett, S.P. & Wanner, G. 1993 Solving ordinary differential equations I: nonstiff problems. Springer Berlin Heidelberg.
  * Halko et al. (2011) Halko, N., Martinsson, P. & Tropp, J. A. 2011 Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review 53 (2), 217–288.
  * Hernandez et al. (2005) Hernandez, V., Roman, J. E. & Vidal, V. 2005 Slepc: A scalable and flexible toolkit for the solution of eigenvalue problems. ACM Transactions on Mathematical Software (TOMS) 31 (3), 351–362.
  * Herrmann et al. (2021) Herrmann, B., Baddoo, P. J., Semaan, R., Brunton, S. L. & McKeon, B. J. 2021 Data-driven resolvent analysis. Journal of Fluid Mechanics 918, A10.
  * House et al. (2022) House, D., Skene, C., Ribeiro, J. H. M., Yeh, C.-A. & Taira, K. 2022 Sketch-based resolvent analysis. AIAA Paper $\\#$2022-3335 .
  * Hunt & Crighton (1991) Hunt, R. E. & Crighton, D. G. 1991 Instability of flows in spatially developing media. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 435 (1893), 109–128.
  * Jeun et al. (2016) Jeun, J., Nichols, J. W. & Jovanović, M. R. 2016 Input-output analysis of high-speed axisymmetric isothermal jet noise. Physics of Fluids 28 (4), 047101.
  * Jiménez (2018) Jiménez, J. 2018 Coherent structures in wall-bounded turbulence. Journal of Fluid Mechanics 842, P1.
  * Jimenez-Gonzalez & Brancher (2017) Jimenez-Gonzalez, J. I. & Brancher, P. 2017 Transient energy growth of optimal streaks in parallel round jets. Physics of Fluids 29 (11), 114101.
  * Jovanovic (2004) Jovanovic, M. R. 2004 Modeling, analysis, and control of spatially distributed systems. University of California, Santa Barbara.
  * Jovanović (2021) Jovanović, M. R. 2021 From bypass transition to flow control and data-driven turbulence modeling: An input–output viewpoint. Annual Review of Fluid Mechanics 53 (1), 010719–060244.
  * Kamal et al. (2023) Kamal, O., Lakebrink, M. T. & Colonius, T. 2023 Global receptivity analysis: physically realizable input–output analysis. Journal of Fluid Mechanics 956, R5.
  * Karban et al. (2020) Karban, U., Bugeat, B., Martini, E., Towne, A., Cavalieri, A. V. G., Lesshafft, L., Agarwal, A., Jordan, P. & Colonius, T. 2020 Ambiguity in mean-flow-based linear analysis. Journal of Fluid Mechanics 900, R5.
  * Kato (2013) Kato, T. 2013 Perturbation theory for linear operators. Springer Science & Business Media.
  * Lesshafft et al. (2019) Lesshafft, L., Semeraro, O., Jaunet, V., Cavalieri, A. V. G. & Jordan, P. 2019 Resolvent-based modeling of coherent wave packets in a turbulent jet. Physical Review Fluids 4 (6), 063901.
  * Li & Malik (1996) Li, F. & Malik, M. R. 1996 On the nature of PSE approximation. Theoretical and Computational Fluid Dynamics 8 (4), 253–273.
  * Lumley (1967) Lumley, J. L. 1967 The structure of inhomogeneous turbulent flows. Atmospheric Turbulence and Radio Wave Propagation pp. 166–178.
  * Marant & Cossu (2018) Marant, M. & Cossu, C. 2018 Influence of optimally amplified streamwise streaks on the Kelvin–Helmholtz instability. Journal of Fluid Mechanics 838, 478–500.
  * Marquet & Larsson (2015) Marquet, O. & Larsson, M. 2015 Global wake instabilities of low aspect-ratio flat-plates. European Journal of Mechanics-B/Fluids 49, 400–412.
  * Martini et al. (2020) Martini, E., Cavalieri, A. V. G., Jordan, P., Towne, A. & Lesshafft, L. 2020 Resolvent-based optimal estimation of transitional and turbulent flows. Journal of Fluid Mechanics 900, A2.
  * Martini et al. (2022) Martini, E., Jung, J., Cavalieri, A. V. G., Jordan, P. & Towne, A. 2022 Resolvent-based tools for optimal estimation and control via the Wiener–Hopf formalism. Journal of Fluid Mechanics 938, E2.
  * Martini et al. (2021) Martini, E., Rodríguez, D., Towne, A. & Cavalieri, A. V. G. 2021 Efficient computation of global resolvent modes. Journal of Fluid Mechanics 919, A3.
  * Mattsson & Nordström (2004) Mattsson, K. & Nordström, J. 2004 Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics 199 (2), 503–540.
  * McKeon (2017) McKeon, B. J. 2017 The engine behind (wall) turbulence: perspectives on scale interactions. Journal of Fluid Mechanics 817, P1.
  * McKeon & Sharma (2010) McKeon, B. J. & Sharma, A. S. 2010 A critical-layer framework for turbulent pipe flow. Journal of Fluid Mechanics 658, 336–382.
  * Moarref et al. (2013) Moarref, R., Sharma, A. S., Tropp, J. A. & McKeon, B. J. 2013 Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels. Journal of Fluid Mechanics 734, 275–316.
  * Monokrousos et al. (2010) Monokrousos, A., Åkervik, E., Brandt, L. & Henningson, D. S. 2010 Global three-dimensional optimal disturbances in the Blasius boundary-layer flow using time-steppers. Journal of Fluid Mechanics 650, 181–214.
  * Morra et al. (2019) Morra, P., Semeraro, O., Henningson, D. S. & Cossu, C. 2019 On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows. Journal of Fluid Mechanics 867, 969–984.
  * Nogueira et al. (2019) Nogueira, P. A. S., Cavalieri, A. V. G., Jordan, P. & Jaunet, V. 2019 Large-scale streaky structures in turbulent jets. Journal of Fluid Mechanics 873, 211–237.
  * Nyquist (1928) Nyquist, H. 1928 Certain topics in telegraph transmission theory. Transactions of the American Institute of Electrical Engineers 47 (2), 617–644.
  * de Pando et al. (2012) de Pando, M. F., Sipp, D. & Schmid, P. J. 2012 Efficient evaluation of the direct and adjoint linearized dynamics from compressible flow solvers. Journal of Computational Physics 231 (23), 7739–7755.
  * Pearson & Pestana (2020) Pearson, J. W. & Pestana, J. 2020 Preconditioners for krylov subspace methods: An overview. GAMM-Mitteilungen 43 (4), e202000015.
  * Pickering et al. (2020) Pickering, E., Rigas, G., Nogueira, P. A. S., Cavalieri, A. V. G., Schmidt, O. T. & Colonius, T. 2020 Lift-up, Kelvin–Helmholtz and Orr mechanisms in turbulent jets. Journal of Fluid Mechanics 896, A2.
  * Pickering et al. (2021) Pickering, E., Rigas, G., Schmidt, O. T., Sipp, D. & Colonius, T. 2021 Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets. Journal of Fluid Mechanics 917, A29.
  * Ran et al. (2021) Ran, W., Zare, A. & Jovanović, M. R. 2021 Model-based design of riblets for turbulent drag reduction. Journal of Fluid Mechanics 906, A7.
  * Ribeiro et al. (2020) Ribeiro, J. H. M., Yeh, C.-A. & Taira, K. 2020 Randomized resolvent analysis. Physical Review Fluids 5 (3), 033902.
  * Saad (2003a) Saad, Y. 2003a Finding exact and approximate block structures for ilu preconditioning. SIAM Journal on Scientific Computing 24 (4), 1107–1123.
  * Saad (2003b) Saad, Y. 2003b Iterative methods for sparse linear systems. SIAM.
  * Sasaki et al. (2022) Sasaki, K., Cavalieri, A. V. G., Hanifi, A. & Henningson, D. S. 2022 Parabolic resolvent modes for streaky structures in transitional and turbulent boundary layers. Physical Review Fluids 7 (10), 104611.
  * Schenk et al. (2001) Schenk, O., Gärtner, K., Fichtner, W. & Stricker, A. 2001 Pardiso: a high-performance serial and parallel sparse linear solver in semiconductor device simulation. Future Generation Computer Systems 18 (1), 69–78.
  * Schmid (2007) Schmid, P. J. 2007 Nonmodal stability theory. Annual Review of Fluid Mechanics 39 (1), 129–162.
  * Schmid (2010) Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28.
  * Schmid (2022) Schmid, P. J. 2022 Dynamic mode decomposition and its variants. Annual Review of Fluid Mechanics 54, 225–254.
  * Schmid & Henningson (2001) Schmid, P. J. & Henningson, D. S. 2001 Stability and transition in shear flows. Springer, New York.
  * Schmidt & Towne (2019) Schmidt, O. T. & Towne, A. 2019 An efficient streaming algorithm for spectral proper orthogonal decomposition. Computer Physics Communications 237, 98–109.
  * Schmidt et al. (2017) Schmidt, O. T., Towne, A., Colonius, T., Cavalieri, A. V. G., Jordan, P. & Brès, G. A. 2017 Wavepackets and trapped acoustic modes in a turbulent jet: coherent structure eduction and global stability. Journal of Fluid Mechanics 825, 1153–1181.
  * Schmidt et al. (2018) Schmidt, O. T., Towne, A., Rigas, G., Colonius, T. & Brès, G. A. 2018 Spectral analysis of jet turbulence. Journal of Fluid Mechanics 855, 953–982.
  * Sinha et al. (2016) Sinha, A., Gudmundsson, K., Xia, H. & Colonius, T. 2016 Parabolized stability analysis of jets from serrated nozzles. Journal of Fluid Mechanics 789, 36–63.
  * Sipp & Marquet (2013) Sipp, D. & Marquet, O. 2013 Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer. Theoretical and Computational Fluid Dynamics 27 (5), 617–635.
  * Sirovich (1987a) Sirovich, L. 1987a Turbulence and the dynamics of coherent structures. i. coherent structures. Quarterly of Applied Mathematics 45 (3), 561–571.
  * Sirovich (1987b) Sirovich, L. 1987b Turbulence and the dynamics of coherent structures. ii. symmetries and transformations. Quarterly of Applied Mathematics 45 (3), 573–582.
  * Skeel (1979) Skeel, R. D. 1979 Scaling for numerical stability in Gaussian elimination. Journal of the ACM (JACM) 26 (3), 494–526.
  * Stewart (1993) Stewart, G. W. 1993 On the early history of the singular value decomposition. SIAM review 35 (4), 551–566.
  * Stewartson & Stuart (1971) Stewartson, K. & Stuart, J. T. 1971 A non-linear instability theory for a wave system in plane Poiseuille flow. Journal of Fluid Mechanics 48 (3), 529–545.
  * Süli & Mayers (2003) Süli, E. & Mayers, D. F. 2003 An introduction to numerical analysis. Cambridge university press.
  * Symon et al. (2019) Symon, S., Sipp, D. & McKeon, B. J. 2019 A tale of two airfoils: resolvent-based modelling of an oscillator versus an amplifier from an experimental mean. Journal of Fluid Mechanics 881, 51–83.
  * Taira et al. (2017) Taira, K., Brunton, S. L., Dawson, S. T. M., Rowley, C. W., Colonius, T., McKeon, B. J., Schmidt, O. T., Gordeyev, S., Theofilis, V. & Ukeiley, L. S. 2017 Modal analysis of fluid flows: An overview. AIAA Journal 55 (12), 4013–4041.
  * Theofilis (2011) Theofilis, V. 2011 Global linear instability. Annual Review of Fluid Mechanics 43, 319–352.
  * Thomareis & Papadakis (2018) Thomareis, N. & Papadakis, G. 2018 Resolvent analysis of separated and attached flows around an airfoil at transitional Reynolds number. Physical Review Fluids 3 (7), 073901.
  * Towne (2016) Towne, A. 2016 Advancements in jet turbulence and noise modeling: accurate one-way solutions and empirical evaluation of the nonlinear forcing of wavepackets. PhD thesis, California Institute of Technology.
  * Towne & Colonius (2015) Towne, A. & Colonius, T. 2015 One-way spatial integration of hyperbolic equations. Journal of Computational Physics 300, 844–861.
  * Towne et al. (2015) Towne, A., Colonius, T., Jordan, P., Cavalieri, A. V. & Bres, G. A. 2015 Stochastic and nonlinear forcing of wavepackets in a Mach 0.9 jet. AIAA Paper $\\#$2015-2217 .
  * Towne et al. (2020) Towne, A., Lozano-Durán, A. & Yang, X. 2020 Resolvent-based estimation of space–time flow statistics. Journal of Fluid Mechanics 883, A17.
  * Towne et al. (2019) Towne, A., Rigas, G. & Colonius, T. 2019 A critical assessment of the parabolized stability equations. Theoretical and Computational Fluid Dynamics 33, 359–382.
  * Towne et al. (2022) Towne, A., Rigas, G., Kamal, O., Pickering, E. & Colonius, T. 2022 Efficient global resolvent analysis via the one-way Navier-Stokes equations. Journal of Fluid Mechanics 948, A9.
  * Towne et al. (2018) Towne, A., Schmidt, O. T. & Colonius, T. 2018 Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. Journal of Fluid Mechanics 847, 821–867.
  * Trefethen & Bau III (1997) Trefethen, L. N. & Bau III, D. 1997 Numerical linear algebra. Siam.
  * Vershynin (2018) Vershynin, R. 2018 High-dimensional probability: An introduction with applications in data science. Cambridge University Press.
  * Wang et al. (2021) Wang, C., Lesshafft, L., Cavalieri, A. V. G. & Jordan, P. 2021 The effect of streaks on the instability of jets. Journal of Fluid Mechanics 910, A14.
  * Wanner & Hairer (1996) Wanner, G. & Hairer, E. 1996 Solving ordinary differential equations II: stiff and differential-algebraic problems. Springer Berlin Heidelberg.
  * Yeh et al. (2020) Yeh, C.-A., Benton, S. I., Taira, K. & Garmann, D. J. 2020 Resolvent analysis of an airfoil laminar separation bubble at Re = 500 000. Physical Review Fluids 5 (8), 083906\.
  * Yeh & Taira (2019) Yeh, C.-A. & Taira, K. 2019 Resolvent-analysis-based design of airfoil separation control. Journal of Fluid Mechanics 867, 572–610.
  * Zhu & Towne (2023) Zhu, M. & Towne, A. 2023 Recursive one-way Navier-Stokes equations with PSE-like cost. Journal of Computational Physics 473, 111744.


# Parallelograms and the VC-dimension of the distance sets

Thang Pham  University of Science, Vietnam National University, Hanoi.
<EMAIL_ADDRESS>
###### Abstract.

In this paper, we study the distribution of parallelograms and rhombi in a
given set in the plane over arbitrary finite fields $\mathbb{F}_{q}^{2}$. As
an application, we improve a recent result due to Fitzpatrick, Iosevich,
McDonald, and Wyman (2021) on the Vapnik–Chervonenkis dimension of the induced
distance graph. Our proofs are based on the discrete Fourier analysis.

###### 2020 Mathematics Subject Classification:

52C10, 42B05, 11T23

##  1\. Introduction

Let $\mathbb{F}_{q}$ be a finite field of order $q$. Given
$E_{1},E_{2},E_{3},E_{4}\subset\mathbb{F}_{q}^{2}$, the set of parallelograms
with one side-length $t$ in $E_{1}\times E_{2}\times E_{3}\times E_{4}$ is
denoted by $Par_{t}(E_{1},E_{2},E_{3},E_{4})$, namely,

$Par_{t}(E_{1},E_{2},E_{3},E_{4}):=\\#\\{(x,y,z,w)\in E_{1}\times E_{2}\times
E_{3}\times E_{4}\colon x-y=z-w,~{}||x-y||=t\\}.$

Here $||x-y||=(x_{1}-y_{1})^{2}+(x_{2}-y_{2})^{2}$. Define

$Par(E_{1},E_{2},E_{3},E_{4}):=\sum_{t\in\mathbb{F}_{q}}Par_{t}(E_{1},E_{2},E_{3},E_{4}),$

as the set of all parallelograms. We also denote the set of rhombi with
vertices in $E_{1}\times E_{2}\times E_{3}\times E_{4}$ of side-length $t$ by
$Rhom_{t}(E_{1},E_{2},E_{3},E_{4})$.

The main purpose of this paper is to study the distribution of rhombi and
parallelograms in $\mathbb{F}_{q}^{2}$ and an application on the problem of
determining the Vapnik–Chervonenkis (VC) dimension of the induced distance
graph in $\mathbb{F}_{q}^{2}$. For other similar configurations, say,
rectangles and simplices, we refer the interested reader to [2, 3, 6, 11, 13,
15] for more discussions.

Our first result is the following.

######  Theorem 1.1.

Given $t\neq 0$. Let $E_{1},E_{2},E_{3},E_{4}\subset\mathbb{F}_{q}^{2}$.
Assume that $|E_{1}||E_{2}|,|E_{3}||E_{4}|\gg q^{3}$, then

$\left|Rhom_{t}(E_{1},E_{2},E_{3},E_{4})-\frac{Par_{t}(E_{1},E_{2},E_{3},E_{4})}{q}\right|\leq
q^{1/2}\left(\frac{|E_{1}||E_{2}|}{q}\right)^{1/2}\left(\frac{|E_{3}||E_{4}|}{q}\right)^{1/2}.$

In addition, if $|E_{1}||E_{2}||E_{3}||E_{4}|\gg q^{7}$, then

$Rhom_{t}(E_{1},E_{2},E_{3},E_{4})\gg\frac{|E_{1}||E_{2}||E_{3}||E_{4}|}{q^{4}}.$

This theorem offers a relation between $Rhom_{t}(E_{1},E_{2},E_{3},E_{4})$ and
$Par_{t}(E_{1},E_{2},E_{3},E_{4})$. In this next theorem, we show that when
$Par(E_{1},E_{2},E_{3},E_{4})$ is large enough, then

$Par_{t}(E_{1},E_{2},E_{3},E_{4})\sim\frac{Par(E_{1},E_{2},E_{3},E_{4})}{q}.$

The formal statement reads as follows.

######  Theorem 1.2.

Let $E_{1},E_{2},E_{3},E_{4}\subset\mathbb{F}_{q}^{2}$ and
$t\in\mathbb{F}_{q}^{*}$. We have

$\left|Par_{t}(E_{1},E_{2},E_{3},E_{4})-\frac{Par(E_{1},E_{2},E_{3},E_{4})}{q}\right|\leq
q^{1/2}(|E_{1}||E_{2}||E_{3}||E_{4}|)^{1/2}.$

In addition, for any $t\neq t^{\prime}$, one has

$\left|Par_{t}(E_{1},E_{2},E_{3},E_{4})-Par_{t^{\prime}}(E_{1},E_{2},E_{3},E_{4})\right|\leq
q^{1/2}(|E_{1}||E_{2}||E_{3}||E_{4}|)^{1/2}.$

We note that given sets $E_{1},E_{2},E_{3},E_{4}$ in $\mathbb{F}_{q}^{2}$, the
number parallelograms with vertices in $E_{1}\times E_{2}\times E_{3}\times
E_{4}$ can be very large. A simple example is that

(1.1) $E_{1}=E_{2}=E_{3}=E_{4}=\mathbb{F}_{q}\cdot\begin{pmatrix}1\\\
0\end{pmatrix}+X\cdot\begin{pmatrix}0\\\ 1\end{pmatrix},$

where $X$ is an arithmetic progression in $\mathbb{F}_{q}$. By a direct
computation, the number of parallelograms with vertices in $E_{1}\times
E_{2}\times E_{3}\times E_{4}$ is about $|X|^{3}q^{3}$. So as long as $|X|\gg
q^{1/2}$, we have $Par_{t}(E_{1},E_{2},E_{3},E_{4})\sim|X|^{3}q^{2}$ and
$Rhom_{t}(E_{1},E_{2},E_{3},E_{4})\sim|X|^{3}q$.

We now discuss an application of the two above theorems. Given
$E\subset\mathbb{F}_{q}^{2}$ and $t\in\mathbb{F}_{q}\setminus\\{0\\}$. Define

$\mathcal{F}(E):=\\{N(x)\colon x\in E\\},$

where $N(x)=\\{y\in E\colon||x-y||=t\\}$.

We say that the system $\mathcal{F}(E)$ shatters a finite set $X\subset E$ if
$\mathcal{F}(E)\cap X$ spans all subsets of $X$. We say the VC-dimension of
$E$ is $d$ if there is a set $X\subset E$ of size $d$ such that it is
shattered by $\mathcal{F}(E)$ and no subset of size $d+1$ is shattered by
$\mathcal{F}(E)$.

The VC-dimension of the system $\mathcal{F}(E)$ has been studied recently by
Fitzpatrick, Iosevich, McDonald, and Wyman [4]. More precisely, they proved
the following theorem.

######  Theorem 1.3.

Given $E\subset\mathbb{F}_{q}^{2}$ with $|E|\gg q^{\frac{15}{8}}$, then the
VC-dimension of $\mathcal{F}(E)$ is equal to three.

In the next theorem, we present a threshold-improvement of this result,
namely, we decrease the exponent $\frac{15}{8}$ to $\frac{13}{7}$.

######  Theorem 1.4.

Let $E\subset\mathbb{F}_{q}^{2}$. Assume that $|E|\gg q^{\frac{13}{7}}$, then
the VC-dimension of $\mathcal{F}$ is equal to three.

We do not know whether the exponent $13/7$ is sharp or not, and it is also
difficult to suggest the next improvement. However, we do know that when the
set $E$ has some nice structures, for instance, as in Example 1.1, a better
threshold can be obtained. Some remarks about proofs of Theorems 1.3 and 1.4
will be mentioned in Section 5.

##  2\. Preliminaries

Let $f\colon\mathbb{F}_{q}^{2}\to\mathbb{C}$ be a complex-valued function. The
Fourier transform of $f$ is defined by

$\widehat{f}(m)=\frac{1}{q^{2}}\sum_{x\in\mathbb{F}_{q}^{2}}f(x)\chi(-x\cdot
m).$

Here $\chi$ is a non-trivial additive character of $\mathbb{F}_{q}$. With this
definition, the Fourier inversion formula and the Plancherel theorem are as
follows:

$f(x)=\sum_{m\in\mathbb{F}_{q}^{2}}\chi(m\cdot
x)\widehat{f}(m),~{}\sum_{m\in\mathbb{F}_{q}^{2}}|\widehat{f}(m)|^{2}=q^{-2}\sum_{x\in\mathbb{F}_{q}^{2}}|f(x)|^{2}.$

Given a set $A\subset\mathbb{F}_{q}^{2}$, by abuse of notation, we denote its
characteristic function by $A(x)$, namely, $A(x)=1$ if $x\in A$ and $0$
otherwise.

For any $j\neq 0$, let $S_{j}$ be the circle centered at the origin of radius
$j$ defined as follows:

$S_{j}:=\\{x\in\mathbb{F}_{q}^{2}\colon x_{1}^{2}+x_{2}^{2}=j\\}.$

The next lemma provides the precise form of the Fourier decay of $S_{j}$ for
any $j\in\mathbb{F}_{q}$. A proof can be found in [7] or [10].

######  Lemma 2.1.

For any $j\in\mathbb{F}_{q}$, we have

$\widehat{S_{j}}(m)=q^{-1}\delta_{0}(m)+\frac{1}{q^{3}}G_{1}^{2}(\eta,\chi)\sum_{r\in{\mathbb{F}}_{q}^{*}}\chi\Big{(}jr+\frac{\|m\|}{4r}\Big{)},$

where $\eta$ is the quadratic character, $\delta_{0}(m)=1$ if $m=(0,\ldots,0)$
and $\delta_{0}(m)=0$ otherwise, and $G_{1}(\eta,\chi)$ is the Gauss sum.

In particular, $|\widehat{S_{j}}(m)|\leq q^{-3/2}$ if $m\neq(0,0)$.

The following lemma is on the number of “unit" distances in a pair of given
sets in $\mathbb{F}_{q}^{2}$. A proof can be found in [8, 14].

######  Lemma 2.2.

Given $t\neq 0$. Let $U$ and $V$ be two sets in $\mathbb{F}_{q}^{2}$. We
denote the number of pairs $(u,v)\in U\times V$ with $||u-v||=t$ by
$N_{t}(U,V)$. Then we have

$\left|N_{t}(U,V)-\frac{|U||V|}{q}\right|\leq 2q^{1/2}\sqrt{|U||V|}.$

##  3\. Proof of Theorem 1.1

For $x$ with $||x||=t$, we denote $E_{i}\cap(E_{j}-x)$ by $E_{ij}^{x}$. The
strategy is as follows: if $u\in E_{1}\cap(E_{2}-x)$ and $v\in
E_{3}\cap(E_{4}-x)$, then $u\in E_{1},u+x\in E_{2},v\in E_{3},v+x\in E_{4}$.
If $||u-v||=t$, then we have a rhombus of side-length $t$ with vertices
$(u,u+x,v+x,v)$.

Therefore, applying Lemma 2.2, one has

$\left|Rhom_{t}(E_{1},E_{2},E_{3},E_{4})-\sum_{||x||=t}\frac{|E_{12}^{x}||E_{34}^{x}|}{q}\right|\leq
q^{1/2}\sum_{||x||=t}\sqrt{|E_{12}^{x}||E_{34}^{x}|}.$

Assume that $|E_{1}||E_{2}|\gg q^{3}$ and $|E_{3}||E_{4}|\gg q^{3}$, then the
second term is bounded from above by

$q^{1/2}\left(\frac{|E_{1}||E_{2}|}{q}\right)^{1/2}\left(\frac{|E_{3}||E_{4}|}{q}\right)^{1/2}.$

The first term is exactly equal to $Par_{t}(E_{1},E_{2},E_{3},E_{4})/q$. Thus,
the first estimate is proved. In the next step, we prove the lower bound for
$Rhom_{t}(E_{1},E_{2},E_{3},E_{4})$. In particular, we want to show that

$Par_{t}(E_{1},E_{2},E_{3},E_{4})\gg\frac{|E_{1}||E_{2}||E_{3}||E_{4}|}{q^{3}}.$

Since $|E_{1}||E_{2}|\gg q^{3}$ and $|E_{3}||E_{4}|\gg q^{3}$, the number of
quadruples $(x,y,z,w)\in E_{1}\times E_{2}\times E_{3}\times E_{4}$ such that
$||x-y||=||z-w||=t$ is $\gg\frac{|E_{1}||E_{2}||E_{3}||E_{4}|}{q^{2}}$. Thus,
there is a direction with at least $|E_{1}||E_{2}||E_{3}||E_{4}|/q^{3}$
parallelograms with one side-length $t$. In other words,

$Par_{t}(E_{1},E_{2},E_{3},E_{4})\gg\frac{|E_{1}||E_{2}||E_{3}||E_{4}|}{q^{3}}.$

This gives

$Rhom_{t}(E_{1},E_{2},E_{3},E_{4})\gg\frac{|E_{1}||E_{2}||E_{3}||E_{4}|}{q^{4}}-\frac{(|E_{1}||E_{2}||E_{3}||E_{4}|)^{1/2}}{q^{1/2}}\gg\frac{|E_{1}||E_{2}||E_{3}||E_{4}|}{q^{4}},$

whenever $|E_{1}||E_{2}||E_{3}||E_{4}|\gg q^{7}$. This completes the proof of
the theorem.

######  Remark 3.1.

In the above proof, when counting the number of quadruples $(u,u+x,v,v+x)\in
E_{1}\times E_{2}\times E_{3}\times E_{4}$, it might happen that $v=u+x$. We
refer to this tuple as a degenerate rhombus. The number of such rhombi is at
most $\sum_{||x||=t}|E_{12}^{x}\ll|E_{1}||E_{2}|/q$. In practise, this
scenario can be avoided completely if the sets $E_{i}$ are disjoint.

##  4\. Proof of Theorem 1.2

Given $E_{1},E_{2},E_{3},E_{4}\subset\mathbb{F}_{q}^{2}$. For any $t\neq 0$,
we first show that $Par_{t}(E_{1},E_{2},E_{3},E_{4})$ is equal to

$q^{6}\sum_{m=-m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(m)\widehat{E_{3}}(m)\widehat{E_{4}}(-m)\widehat{S_{t}}(0)+q^{6}\sum_{m\neq-m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(-m^{\prime})\widehat{E_{3}}(m)\widehat{E_{4}}(m^{\prime})\widehat{S_{t}}(-m-m^{\prime}).$

Indeed,

$\displaystyle
Par_{t}(E_{1},E_{2},E_{3},E_{4})=\sum_{v\in\mathbb{F}_{q}^{2}}\sum_{x,y\in\mathbb{F}_{q}^{2}}E_{1}(x)E_{2}(x+v)E_{3}(y)E_{4}(y+v)S_{t}(v)$
$\displaystyle=\sum_{x,y,v}\sum_{m,m^{\prime}}E_{1}(x)E_{2}(y)\widehat{E_{3}}(m)\widehat{E_{4}}(m^{\prime})\chi(m(x+v))\chi(m^{\prime}(y+v))S_{t}(v)$
$\displaystyle=q^{6}\sum_{m,m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(-m^{\prime})\widehat{E_{3}}(m)\widehat{E_{4}}(m^{\prime})\widehat{S_{t}}(-m-m^{\prime})$
$\displaystyle=q^{6}\sum_{m=-m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(m)\widehat{E_{3}}(m)\widehat{E_{4}}(-m)\widehat{S_{t}}(0)+q^{6}\sum_{m\neq-m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(-m^{\prime})\widehat{E_{3}}(m)\widehat{E_{4}}(m^{\prime})\widehat{S_{t}}(-m-m^{\prime})$
$\displaystyle=I+II.$

To simplify the notations, we write $Par_{t}$ instead of
$Par_{t}(E_{1},E_{2},E_{3},E_{4})$.

For all non-zero $t$ and $t^{\prime}$ with $t\neq t^{\prime}$, we observe that
$I(Par_{t})=I(Par_{t^{\prime}})$. On the other hand, using Lemma 2.1, we have

$\displaystyle|II(Par_{t})|=|q^{6}\sum_{m\neq-m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(-m^{\prime})\widehat{E_{3}}(m)\widehat{E_{4}}(m^{\prime})\widehat{S_{t}}(-m-m^{\prime})|$
$\displaystyle\leq
q^{\frac{9}{2}}\sum_{m\neq-m^{\prime}}|\widehat{E_{1}}(-m)||\widehat{E_{2}}(-m^{\prime})||\widehat{E_{3}}(m)||\widehat{E_{4}}(m^{\prime})|$
$\displaystyle\leq
q^{9/2}\left(\sum_{m}|\widehat{E_{1}}(m)|^{2}\right)^{1/2}\left(\sum_{m}|\widehat{E_{2}}(m)|^{2}\right)^{1/2}\left(\sum_{m}|\widehat{E_{3}}(m)|^{2}\right)^{1/2}\left(\sum_{m}|\widehat{E_{4}}(m)|^{2}\right)^{1/2}$
$\displaystyle\leq q^{1/2}(|E_{1}||E_{2}||E_{3}||E_{4}|)^{1/2}.$

This implies that

$|Par_{t}-Par_{t^{\prime}}|\leq q^{1/2}(|E_{1}||E_{2}||E_{3}||E_{4}|)^{1/2}.$

As proved above,

$Par_{t}=q^{6}\left(\sum_{m=-m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(m)\widehat{E_{3}}(m)\widehat{E_{4}}(-m)\widehat{S_{t}}(0)+\sum_{m\neq-m^{\prime}}\widehat{E_{1}}(-m)\widehat{E_{2}}(-m^{\prime})\widehat{E_{3}}(m)\widehat{E_{4}}(m^{\prime})\widehat{S_{t}}(-m-m^{\prime})\right).$

The first term can be computed explicitly, namely,

$\displaystyle I$
$\displaystyle=\frac{1}{q^{3}}\sum_{m}\sum_{x,y,z,w}E_{1}(x)E_{2}(y)E_{3}(z)E_{4}(w)\chi(m\cdot(x+y-z-w))=\frac{}{}\frac{1}{q}\sum_{x,y,z,w}1_{x+y-z-w}=\frac{Par(E_{1},E_{2},E_{3},E_{4})}{q}.$

Since the absolute value of the second term is at most
$q^{1/2}(|E_{1}||E_{2}||E_{3}||E_{4}|)^{1/2}$, we conclude that

$\left|Par_{t}(E_{1},E_{2},E_{3},E_{4})-\frac{Par(E_{1},E_{2},E_{3},E_{4})}{q}\right|\leq
q^{1/2}(|E_{1}||E_{2}||E_{3}||E_{4}|)^{1/2}.$

This completes the proof.

##  5\. Proof of Theorem 1.4

To prove Theorem 1.4, we need to find three distinct vertices
$x^{1},x^{2},x^{3}\in E$ and vertices
$y^{12},y^{13},y^{123},y^{1},y^{2},y^{3},y^{0}$ in $E$ such that the following
hold:

$||y^{123}-x^{1}||=||y^{123}-x^{2}||=||y^{123}-x^{3}||=t.$
$||y^{12}-x^{1}||=||y^{12}-x^{2}||=t,~{}||y^{12}-x^{3}||\neq t.$
$||y^{13}-x^{1}||=||y^{13}-x^{3}||=t,~{}||y^{13}-x^{2}||\neq t.$
$||y^{23}-x^{2}||=||y^{23}-x^{3}||=t,~{}||y^{23}-x^{1}||\neq t.$
$||y^{1}-x^{1}||=t,~{}||y^{1}-x^{2}||\neq t,~{}||y^{1}-x^{3}||\neq t.$
$||y^{2}-x^{1}||\neq t,~{}||y^{2}-x^{2}||=t,~{}||y^{2}-x^{3}||\neq t.$
$||y^{3}-x^{1}||\neq t,~{}||y^{3}-x^{2}||\neq t,~{}||y^{3}-x^{3}||=t.$
$||y^{0}-x^{1}||\neq t,~{}||y^{0}-x^{2}||\neq t,~{}||y^{0}-x^{3}||\neq t.$

This can be described clearly in the following picture.

In this picture, if two vertices are of length $t$ then they are connected by
a solid segment, otherwise by a dotted segment.

Let $E^{\prime}\subset E$ be the set of vertices $u$ such that the number of
$v\in E$ and $||u-v||=t$ is at least $\frac{|E|}{2q}$. By Lemma 2.2, we have

$\frac{|E\setminus E^{\prime}||E|}{q}-2q^{1/2}\sqrt{|E\setminus
E^{\prime}||E|}\leq N(E\setminus E^{\prime},E)<\frac{|E\setminus
E^{\prime}||E|}{2q}.$

This implies

$|E\setminus E^{\prime}|\leq\frac{4q^{3}}{|E|^{2}}=o(|E|),$

when $|E|\gg q^{13/7}$. A similar argument also shows that the set of vertices
with at least $2|E|/q$ neighbors is $o(|E|)$. Thus, by abuse of notations, we
might assume that each vertex in $E$ has at least $|E|/2q$ and at most
$2|E|/q$ neighbors in $E$.

We now use Lemma 2.2 again to find $u\in\mathbb{F}_{q}^{2}$ with $||u||=t$
such that $|E\cap(E-u)|\gg|E|^{2}/q^{2}$. Indeed, we have the number of pairs
$(x,y)\in E^{2}$ such that $||x-y||=t$ is $(1-o(1))|E|^{2}/q$ and there are at
most $q+1$ points on any circle of radius $t$, this implies the existence of
such $u$. In the rest of the proof, we denote $E_{u}=E\cap(E-u)$, so
$|E_{u}|\gg|E|^{2}/q^{2}$. Notice that $E_{u}\subset E$.

The next step is to find a rhombus with vertices
$(x^{1},y^{123},x^{3},y^{13})\in E_{u}\times E_{u}\times E_{u}\times E$ such
that $x^{1}-y^{123}\not\in\\{u,-u\\}$, $y^{123}-y^{13}\not\in\\{u,-u\\}$, and
$y^{123}-x^{3}\not\in\\{u,-u\\}$.

Let $A,B,C\subset E_{u}$ and $D\subset E$ be disjoint sets with
$|A|=|B|=|C|\sim|E_{u}|$ and $|D|\sim|E|$. We want to find such a rhombus with
vertices in $A\times B\times C\times D$. This helps us to avoid the degeneracy
that might happen as noticed in Remark 3.1.

Applying Theorem 1.1, the number of rhombi with one side-length $t$ is at
least

$\frac{|A||B||C||D|}{q^{4}},$

under $|E|\gg q^{13/7}$. Using the fact that the two circles of radius $t$
intersect in at most $2$ points, the number of rhombi
$(x^{1},y^{123},x^{3},y^{13})$ with $y^{123}-y^{13}\in\\{u,-u\\}$ is at most
$2|E_{u}|$. The number of rhombi $(x^{1},y^{123},x^{3},y^{13})$ with
$y^{123}-x^{1}\in\\{u,-u\\}$ is at most $2N_{t}(E_{u},E_{u})$. With the same
argument, the number of rhombi $(x^{1},y^{123},x^{3},y^{13})$ with
$y^{123}-x^{3}\in\\{u,-u\\}$ is at most $2N_{t}(E_{u},E_{u})$. Since $|E|\gg
q^{13/7}$, $|E_{u}|\gg q^{12/7}$. Thus
$N_{t}(E_{u},E_{u})=(1+o(1))|E_{u}|^{2}/q$. In total, the number of these
rhombi is much smaller than $q^{-4}|A||B||C||D|$ whenever $|E|\gg q^{5/3}$.

Hence, there are many rhombi $(x^{1},y^{123},x^{3},y^{13})\in E_{u}\times
E_{u}\times E_{u}\times E$ such that $x^{1}-y^{123}\not\in\\{u,-u\\}$,
$y^{123}-y^{13}\not\in\\{u,-u\\}$, and $y^{123}-x^{3}\not\in\\{u,-u\\}$. We
fix one of them. Set

$y^{12}=x^{1}+u,~{}x^{2}=y^{123+u},~{}y^{23}=x^{3}+u.$

Since $E_{u}=E\cap(E-u)$, we have $y^{12},x^{2},y^{23}\in E$. Since
$y^{123}-x^{1},y^{123}-x^{3}\not\in\\{u,-u\\}$, the three vertices
$x^{1},x^{2},x^{3}$ are distinct. By our constraints, we can see that all
vertices $x^{1},x^{2},x^{3},y^{12},y^{13},y^{23},y^{123}$ are distinct.

The next step is to check $||y^{12}-x^{3}||\neq t$ and $||y^{23}-x^{1}||\neq
t$. This is clear otherwise two distinct circles of radius $t$ intersect in
more than two points.

The final step is to choose $y^{1},y^{2},y^{3}$ with the desired property.
This is also clear since we have assumed that each vertex of $E$ has at least
$|E|/2q$ neighbors in $E$. Notice that $y^{i}\neq y^{j}$ for $i\neq j$,
otherwise, we again have three points in the intersection of two circles of
radius $t$.

Since each of vertices $x^{1},x^{2},x^{3}$ has at most $2|E|/q$ neighbors, the
vertex $y_{0}$ can be chosen arbitrary in $E\setminus(N(x^{1})\cup
N(x^{2})\cup N(x^{3}))$.

This completes the proof of the theorem.

###  5.1. Discussions

In this paragraph, we want to discuss briefly about the proof of Theorem 1.4
and that of Theorem 1.3 in [4]. The main difference comes from the following
step: in stead of counting the number of rhombi with vertices in $E_{u}\times
E_{u}\times E_{u}\times E_{u}$ as in [4], we replace one set $E_{u}$ by a
bigger set $E$, this gives us more room to find desired configurations. Thus,
the exponent $\frac{15}{8}$ can be decreased to $\frac{13}{7}$. In a recent
work of Pham, Senger, Tait, and Thu-Huyen [12], this VC-dimension notation has
been studied in a more general setting for pseudo-random graphs. As mentioned
in [12], when applying their result to the distance graph, Theorem 1.3 can not
be improved. This can be explained in detail as follows. In [12], to determine
the VC-dimension, one has to bound the number of $C_{4}$ from below, but for
pseudo-random graphs with the underlying vertex set $\mathbb{F}_{q}^{2}$,
there is a possibility that it does not contain any non-trivial $C_{4}$. For
example, the graph with the vertex set $\mathbb{F}_{q}^{2}$ and the two
vertices $(x_{1},x_{2})$ and $(y_{1},y_{2})$ are connected by an edge if
$x_{1}y_{1}+x_{2}y_{2}=1$. This graph is a pseudo-random graph. The adjacent
relation can be viewed as the point $(x_{1},x_{2})$ belongs to the line
defined by $y_{1}X+y_{2}Y=1$. Since two lines are either parallel or intersect
in one point, we do not have any non-trivial $C_{4}$ in this graph, i.e. a
$C_{4}$ with four distinct vertices. This explains why in
$\mathbb{F}_{q}^{2}$, for the distance graph, the method in [12] is not
applicable.

If the reader is interested in the VC-dimension of other geometric
configurations, we refer to some recent papers [1, 5, 9] for more details.

##  6\. Acknowledgements

T. Pham would like to thank to the Vietnam Institute for Advanced Study in
Mathematics (VIASM) for the hospitality and for the excellent working
conditions.

## References

  * [1] R. Ascoli, L. Betti, J. Cheigh, A. Iosevich, R. Jeong, X. Liu, B. McDonald, W. Milgrim, S. J. Miller, F. R. Acosta, and S.V. Iannuzzelli, VC-Dimension and Distance Chains in $\mathbb{F}_{q}^{d}$, arXiv:2210.03058, 2023.
  * [2] M. Bennett, A. Iosevich, and J. Pakianathan, Three-point configurations determined by subsets of $\mathbb{F}_{q}^{2}$ via the Elekes-Sharir Paradigm, Combinatorica 34.6: 689–706, 2014.
  * [3] M. Bennett, D. Hart, A. Iosevich, J. Pakianathan, and M. Rudnev, Group actions and geometric combinatorics in $\mathbb{F}_{q}^{d}$, Forum Mathematicum, vol. 29, no. 1, pp. 91-110. De Gruyter, 2017.
  * [4] D. Fitzpatrick, A. Iosevich, B. McDonald, and E. Wyman, The VC-dimension and point configurations in $\mathbb{F}_{q}^{2}$, arXiv:2108.13231, 2023.
  * [5] B. McDonald, A. Sahay, and E. Wyman, The VC-dimension of quadratic residues in finite fields, arXiv:2210.03789, 2022.
  * [6] N. Lyall and A. Magyar, Weak hypergraph regularity and applications to geometric Ramsey theory, Transactions of the American Mathematical Society, Series B, 9(5) (2022): 160–207.
  * [7] D. Hart, A. Iosevich, D. Koh, and M. Rudnev, Averages over hyperplanes, sum-product theory in vector spaces over finite fields and the Erdős-Falconer distance conjecture, Transactions of the American Mathematical Society 363, no. 6 (2011): 3255-3275.
  * [8] A. Iosevich and M. Rudnev, Erdős distance problem in vector spaces over finite fields, Transactions of the American Mathematical Society 359(12) (2007): 6127–6142.
  * [9] A. Iosevich, B. McDonald, and M. Sun, Dot products in $\mathbb{F}_{q}^{3}$ and the Vapnik-Chervonenkis dimension, Discrete Mathematics 346(2023), no. 1, 113096.
  * [10] D. Koh, and H.S. Sun, Distance sets of two subsets of vector spaces over finite fields, Proceedings of the American Mathematical Society 143.4 (2015): 1679-1692.
  * [11] D. Koh, S. Lee, T. Pham, and C-Y. Shen, Configurations of rectangles in $\mathbb{F}_{q}^{2}$, Finite Fields and Their Applications, 86(2023): 102147.
  * [12] T. Pham, S. Senger, M. Tait, and N. Thu-Huyen, VC-dimension and pseudo-random graphs, arXiv:2303.07878, 2023.
  * [13] I. D. Shkredov, On some problems of Euclidean Ramsey theory, Analysis Mathematica, 41(4) (2015): 299-310.
  * [14] I. G. Shparlinski, On the set of distances between two sets over finite fields, International Journal of Mathematics and Mathematical Sciences, 2006.
  * [15] L. A. Vinh, On kaleidoscopic pseudo-randomness of finite Euclidean graphs, Discus- siones Mathematicae: Graph Theory, 32(2) (2012): 279-287.