High-dimensional limit of one-pass SGD on least squares High-dimensional limit
of one-pass SGD on least squares. Elizabeth Collins–Woodfin111McGill
University, Canada<EMAIL_ADDRESS>https://sites.google.com/view/e-collins-woodfin and Elliot Paquette 222McGill
University, Canada<EMAIL_ADDRESS>http://elliotpaquette.github.iostochastic gradient descent ; random matrix
theory ; optimization ; stochastic differential equations 60H30 TBD TBD 0 0 0
0 We give a description of the high-dimensional limit of one-pass single-batch
stochastic gradient descent (SGD) on a least squares problem. This limit is
taken with non-vanishing step-size, and with proportionally related number of
samples to problem-dimensionality. The limit is described in terms of a
stochastic differential equation in high dimensions, which is shown to
approximate the state evolution of SGD. As a corollary, the statistical risk
is shown to be approximated by the solution of a convolution-type Volterra
equation with vanishing errors as dimensionality tends to infinity. The sense
of convergence is the weakest that shows that statistical risks of the two
processes coincide. This is distinguished from existing analyses by the type
of high-dimensional limit given as well as generality of the covariance
structure of the samples.

##  1 Introduction

Stochastic optimization methods are the modern day standard for many large-
scale computational tasks, especially those that arise in machine learning.
There is a long history of analyses of these algorithms, beginning with the
seminal work of [RM51], which focused on long-time behavior in a fixed
dimensional space. However, modern applications of stochastic optimization
have motivated a different regime of analysis, where the problem
dimensionality grows proportionally with the run-time of the algorithm.

In this article, we derive the exact scaling behavior of stochastic gradient
descent (SGD) on a least squares problem, in the _one-pass_ setting (see
below) when dimension tends to infinity. We further draw a comparison to the
recent work [Paq+22, Paq+21], in which the _multi-pass_ version of this
problem was considered.

#### Stochastic gradient descent for empirical risk minimization

Most versions of (minibatch) SGD can be formulated in the context of _finite-
sum problems_ :

$\min_{\mathbold{x}\in\mathbb{R}^{d}}\leavevmode\nobreak\
\bigg{\\{}f(\mathbold{x})\coloneqq\frac{1}{n}\sum_{i=1}^{n}f_{i}(\mathbold{x})\bigg{\\}}.$
(1.1)

Empirical risk minimization fits in this context by supposing that there are
$n$ independent samples from some data distribution, and each $f_{i}$
represents the loss of how the parameters $\mathbold{x}$ in some model fit the
$i$-th datapoint. In this article we will exclusively consider the case of
linear regression with $\ell^{2}$–regularizer. So we suppose that there are
$n$ iid samples $\left((\mathbold{a}_{i},\mathbold{b}_{i}):1\leq i\leq
n\right)$ from some distrbution $\mathcal{D}$, with some assumptions to be
specified. We arrange this data into a design matrix $\mathbold{A}$ and label
vector $\mathbold{b}$, whose $i$-th row is given by $\mathbold{a}_{i}$.
Finally, we specify the functions $f_{i}$ in (1.1) by setting

$f_{i}(\mathbold{x})=\tfrac{1}{2}(\mathbold{a}_{i}\cdot\mathbold{x}-\mathbold{b}_{i})^{2}+\tfrac{\delta}{2}\|\mathbold{x}\|^{2}.$

The parameter $\delta\geq 0$ is fixed and is the strength of the regularizer
and throughout $\|\cdot\|$ will be the Euclidean norm.

Minibatch stochastic gradient descent in this context can be described as

$\displaystyle\mathbold{x}_{k+1}$
$\displaystyle=\mathbold{x}_{k}-\gamma_{k}\nabla f_{i_{k}}(\mathbold{x}_{k})$
(1.2)
$\displaystyle=\mathbold{x}_{k}-\gamma_{k}\mathbold{A}^{T}\mathbold{e}_{i_{k}}\mathbold{e}_{i_{k}}^{T}(\mathbold{A}\mathbold{x}_{k}-\mathbold{b})-\gamma_{k}\delta\mathbold{x}_{k}\,$

where $\\{\gamma_{k}\\}$ are stepsize parameters, $\mathbold{e}_{i}$ is the
$i$-th standard basis vector, and $\\{i_{k}\\}$ is a sequence of choices data.

In this article we consider the _one-pass_ case, in which $i_{k}=k$ but the
algorithm is terminated after $n$ steps. In practice, the order of the data
points might be shuffled once before, but in the setting we have posed, with
iid data, there is no point to including this additional randomization. There
are other choices for how to pick $i_{k},$ and we highlight three of them, all
of which are _multi-pass_ variants.

In _random (with replacement) sample_ SGD, each $i_{k}$ is chosen uniformly at
random from $\\{1,2,\dots,n\\}$. This is the setting considered in [Paq+22],
and we shall refer simply to this flavor of SGD simply as _multi-pass_ SGD in
the bulk of the paper. But for context, we also mention _single shuffle_ SGD,
in which one takes $i_{k}=k\text{ mod }n$, and so only differs from one-pass
SGD in that the algorithm performs the same operations every _epoch_ 333We
take epoch to mean $n$ steps of the algorithm in all these contexts (including
the with-replacement case). In _random shuffle_ SGD, one modifies the above
strategy by randomly permuting the data between each epoch. All of these
strategies are extensively studied in the optimization literature: it is
generally thought that the _single shuffle_ and _random shuffle_ strategies
are faster than the random sample strategy [YSJ21] (see also [RR12, GOP21,
SS20, AYS20]).

The one-pass case is the fundamental point of comparison for all of these
methods, being both simpler phenomenologically and also representing an
idealization of SGD, in which the run-time of the algorithm is the amount
data. Appropriately, running for longer (meaning increasing $n$) can only
improve the statistical performance of the SGD estimator $\mathbold{x}_{n}$,
in which context this is usually referred to as _streaming SGD_.

The performance of SGD is measured through the _population risk_ $\mathscr{P}$
and sometimes through an $\ell_{2}$–regularized risk $\mathscr{R}$, which are
given by

$\displaystyle\mathscr{P}(\mathbold{x})\coloneqq\tfrac{1}{2}\mathbb{E}_{(\mathbold{a},\mathbold{b})}(\mathbold{a}\cdot\mathbold{x}-\mathbold{b})^{2},\quad$
$\displaystyle(\mathbold{a},\mathbold{b})\sim\mathcal{D},$ (1.3)
$\displaystyle\mathscr{R}(\mathbold{x})\coloneqq\tfrac{1}{2}\mathbb{E}_{(\mathbold{a},\mathbold{b})}(\mathbold{a}\cdot\mathbold{x}-\mathbold{b})^{2}+\tfrac{\delta}{2}\|\mathbold{x}\|^{2},\quad$
$\displaystyle(\mathbold{a},\mathbold{b})\sim\mathcal{D}.$ (1.4)

This regularized risk appears naturally as the mean behavior of one-pass SGD,
in that

$\mathbold{x}_{k+1}=\mathbold{x}_{k}-\gamma_{k}(\nabla\mathscr{R}(\mathbold{x}_{k})+\xi_{k+1}),$

for martingale increments $(\xi_{k}:1\leq k\leq n)$. Another natural
statistical setting to consider is out-of-distributional regression, in which
case we would measure the performance of SGD trained on $\mathcal{D}$ but
tested, as in (1.3) after replacing $\mathcal{D}$ with another distribution
$\mathcal{D}^{\prime}$. We shall not pursue this case in detail, but we note
that all of the above examples are some quadratic functionals of the SGD state
$\mathbold{x}$.

#### Data and stepsize assumptions

The goal of this analysis is to allow the number of samples $n$ to be large
and proportional to the dimension of the problem, here $d$. This means that
the data must be normalized to be nearly dimension independent. Further, we
shall need good tail properties of some of the random variables involved, and
so we recall the Orlicz norms $\|\cdot\|_{\psi_{p}}$ for $p\geq 1$ which are
given by

$\|X\|_{\psi_{p}}=\inf\\{t:\mathbb{E}e^{|X|^{p}/t^{p}}\leq 2\\}.$

We refer the reader to [Ver18] for further exposition, properties and
equivalent formulations.

We shall suppose throughout that under $\mathcal{D}$, the labels are given by
an underlying linear model with noise. Formally, we suppose that: For
$(\mathbold{a},\mathbold{b})$ sampled from $\mathcal{D}$, conditionally on
$\mathbold{a}$, the distribution of $\mathbold{b}$ is given by
$\mathbold{a}\cdot\tilde{\mathbold{x}}+w$ where $w$ is mean $0$, variance
$\eta^{2}\geq 0$ and is subgaussian with $\|w\|_{\psi_{2}}\leq
d^{\varepsilon}.$ The ground truth $\tilde{\mathbold{x}}$ is assumed to have
norm at most $d^{\varepsilon}$. The constant $\varepsilon$ will be small and
fixed throughout. Anything less than $\tfrac{1}{18}$ will do.

The data covariance is assumed to be normalized in such a way that it is
almost bounded in norm, which is to say: The covariance matrix
$\mathbold{K}\coloneqq\mathbb{E}[\mathbold{a}\mathbold{a}^{T}]$ has operator
norm bounded independent of $d.$ Note that while we do not explicitly assume
that $\mathbold{a}$ is centered, the mean would have to be small in some sense
to achieve the assumption above.

Finally, we suppose that $\mathbold{a}$ has good tail properties, namely that
The data vector $\mathbold{a}$ satisfies that, for any deterministic
$\mathbold{x}$ of norm less than $1$,
$\|\mathbold{a}\cdot\mathbold{x}\|_{\psi_{2}}\leq d^{\varepsilon}$, and the
data vector $\mathbold{a}$ satisfies the Hanson-Wright inequality: for all
$t\geq 0$ and for any deterministic matrix $\mathbold{B}$

$\mathbb{P}\left(\left|\mathbold{a}^{T}\mathbold{B}\mathbold{a}-\mathbb{E}\mathbold{a}^{T}\mathbold{B}\mathbold{a}\right|\geq
t\right)\leq
2\exp\left(-\min\left\\{\frac{t^{2}}{d^{4\varepsilon}||\mathbold{B}||_{HS}^{2}},\;\frac{t}{d^{2\varepsilon}||\mathbold{B}||}\right\\}\right).$

We remark that these assumptions hold for two important settings:

  1. (a)

when $\mathbold{a}=\sqrt{\mathbold{K}}\mathbold{u}$ where $\mathbold{K}$ is
some deterministic matrix of bounded operator norm and $\mathbold{u}$ is a
vector of iid subgaussian random variables or

  2. (b)

when $\mathbold{a}$ is a vector with the _convex concentration property_ , see
[Ada15] for details.

There are natural examples of the second case, such as random features models
[RR08] with Lipschitz activation functions (see also [Paq+22, Proposition 6.2]
for specifics in the case of random features). We note that by truncation, it
is also possible to work in the setting (a) above but solely under a uniform
bound on a sufficiently high but finite moment, but we do not pursue this.

Finally, the step-size parameters $\gamma_{k}$ must be normalized
approporiately: The stepsize $\gamma_{k}=\tfrac{\gamma}{d}$ for all $k$ and
fixed $\gamma>0$. We note that we may also pick $\gamma_{k}=\gamma(k/d)/d$ for
a bounded continuous function $\gamma:[0,\infty)\to[0,\infty)$, and this leads
to no change anywhere in the arguments.

In light of all these assumptions, we note that SGD finally reduces to the
following stochastic recurrence

$\mathbold{x}_{k}-\tilde{\mathbold{x}}=(\mathbold{I}(1-\tfrac{\gamma\delta}{d})-\gamma\mathbold{m}_{k}\mathbold{m}_{k}^{T})(\mathbold{x}_{k-1}-\tilde{\mathbold{x}})-\tfrac{\gamma\delta}{d}\tilde{\mathbold{x}}+\gamma\mathbold{m}_{k}\eta_{k},$
(1.5)

where $\eta_{k}=w_{k}/\sqrt{d}$ and
$\mathbold{m}_{k}=\mathbold{a}_{k}/\sqrt{d}$.

#### Homogenized SGD

Our theorem is most easily formulated as showing that the state of SGD can be
compared to a certain diffusion model in high dimensions. Homogenized SGD is
defined to be a continuous time process with initial condition
$\mathbold{X}_{0}=\mathbold{x}_{0}$ that solves the stochastic differential
equation

$\mathrm{d}\mathbold{X}_{t}=-\gamma\nabla\mathscr{R}(\mathbold{X}_{t})\mathrm{d}t+\gamma\sqrt{\tfrac{2}{d}\mathscr{P}(\mathbold{X}_{t})\mathbold{K}}\mathrm{d}B_{t}$
(1.6)

where $B_{t}$ is standard Brownian motion in dimension $d$, and $\mathscr{R}$
and $\mathscr{P}$ are the regularized and unregularized risks, respectively
(recall (1.3) and (1.4)).

Our main theorem shows that for quadratic statistics (in particular the risks
(1.3) and (1.4)), homogenized SGD and SGD are interchangeable to leading
order. We use the probabilistic modifier _with overwhelming probability_ to
mean a statement holds except on an event of probability at most
$e^{-\omega(\log d)}$ where $\omega(\log d)$ tends to $\infty$ faster than
$\log d$ as $d\to\infty.$ We further introduce a norm $\|\cdot\|_{C^{2}}$ on
quadratic functions $q:\mathbb{R}^{d}\to\mathbb{C}$

$\|q\|_{C^{2}}\coloneqq\|\nabla^{2}q\|+\|\nabla q(0)\|+|q(0)|,$

with the norms on the right hand side being given by the operator and
Euclidean norm respectively.

######  Theorem 1.1.

For any quadratic $q:\mathbb{R}^{d}\to\mathbb{R}$, and for any deterministic
initialization $\mathbold{x}_{0}$ with $\|\mathbold{x}_{0}\|\leq 1$, there is
a constant $C\left(\|\mathbold{K}\|\right)$ so that the processes
$\\{\mathbold{x}_{k}\\}_{k=0}^{n}$ and $\\{\mathbold{X}_{t}\\}_{t=0}^{n/d}$
satisfy for any $n$ satisfying $n\leq d\log d/C(\|\mathbold{K}\|)$

$\sup_{0\leq k\leq
n}\Big{|}q(\mathbold{x}_{k})-q(\mathbold{X}_{k/d})\Big{|}<\|q\|_{C^{2}}\cdot
e^{C\left(\|\mathbold{K}\|\right)\frac{n}{d}}\cdot
d^{-\frac{1}{2}+9\varepsilon}$ (1.7)

with overwhelming probability.

The processes $\mathbold{x}_{k}$ and $\mathbold{X}_{t}$ are independent, and
hence this is also a statement about concentration. In particular, the
statement is also true if we replace $q(\mathbold{X}_{k/d})$ by
$\mathbb{E}q(\mathbold{X}_{k/d})$.

#### Explicit risk curves

Using existing theory (see [Paq+22, Theorem 1.1]),
$\mathscr{P}(\mathbold{X}_{k/d})$ and $\mathscr{R}(\mathbold{X}_{k/d})$ can be
seen to concentrate around their means, which solve a convolution Volterra
equation. Specifically, $\mathbb{E}\mathscr{P}(\mathbold{X}_{t})=\Psi(t)$ and
$\mathbb{E}\mathscr{R}(\mathbold{X}_{t})=\Omega(t)$ where

$\begin{pmatrix}\Psi(t)\\\
\Omega(t)\end{pmatrix}=\begin{pmatrix}\mathscr{P}(\mathscr{X}_{\gamma t})\\\
\mathscr{R}(\mathscr{X}_{\gamma t})\\\
\end{pmatrix}+\int_{0}^{t}\begin{pmatrix}K(t-s;\nabla^{2}\mathscr{R})\Psi(s)\\\
K(t-s;\nabla^{2}\mathscr{P})\Psi(s)\end{pmatrix}\mathrm{d}s,\quad\begin{cases}K(t;\mathbold{M})\coloneqq\tfrac{\gamma^{2}}{d}\operatorname{tr}(\mathbold{K}\mathbold{M}e^{-2\gamma
t(\mathbold{K}+\delta\mathbold{I})}),\\\ \mathrm{d}\mathscr{X}_{\gamma
t}\coloneqq-\gamma\nabla\mathscr{R}(\mathscr{X}_{\gamma
t}),\quad\mathscr{X}_{0}=\mathbold{X}_{0}.\end{cases}$

The process $\mathscr{X}_{t}$ is gradient flow, and for the problems here, is
explicitly solvable. Note the equation for $\Psi$ is autonomous, and the
solution of $\Omega$ is then solvable in terms of it. From this Volterra
equation, it is easy to derive convergence rates, as well as convergence
thresholds, as well as expressions for the limiting risk (in the double
scaling limit, $d\to\infty$ followed by $t\to\infty$). See the discussions in
[Paq+21] and see Figure 1 for an example.

#### Comparison to the multi-pass case

In [Paq+22], the analog of Theorem 1.1 was proven for the (random sample)
multi-pass case. In that case, the diffusion is different. Introduce the
_empirical risk_
$\mathscr{L}(\mathbold{x})=\tfrac{1}{2n}\|\mathbold{A}\mathbold{x}-\mathbold{b}\|^{2}$
and the regularized empirical risk
$f(\mathbold{x})=\mathscr{L}(\mathbold{x})+\delta\|\mathbold{x}\|^{2}$. Then
the homogenized SGD for the multi-pass case becomes

$\mathrm{d}\mathbold{X}_{t}=-\gamma\nabla
f(\mathbold{X}_{t})\mathrm{d}t+\gamma\sqrt{\tfrac{2}{d}\mathscr{L}(\mathbold{X}_{t})(\tfrac{1}{n}\mathbold{A}^{T}\mathbold{A})}\mathrm{d}B_{t}$
(1.8)

Hence the difference between the multi-pass and one-pass cases is that the
population risks are traded for the empirical risk and the data covariance
matrix $\mathbold{K}$ is traded for the empirical data covariance matrix
$(\tfrac{1}{n}\mathbold{A}^{T}\mathbold{A})$. Note that if one conditions on
the data $(\mathbold{A},\mathbold{b})$, then multi-pass SGD in fact _is_
streaming SGD, but with a finitely supported data distribution (specifically
the empirical distribution of data); however the empirical distribution of
data is very far from satisfying Assumption 1.

Figure 1: Risk curves for a simple linear regression problem in $d=2000$.
Multi-pass SGD, its high dimensional equivalent (the expected risk under
homogenized SGD, i.e. “Volterra”), Streaming SGD (i.e. one-pass with varying
dataset size), and the expected risk of homogenized SGD (“Streaming Volterra”)
are all plotted. Risk levels for streaming SGD at various levels $n$ are
plotted for comparison against the corresponding multi-pass version. Note that
at smaller dataset sizes, multi-pass SGD improves greatly over one-pass SGD.
At higher dataset sizes, they are similar and in fact multi-pass SGD always
underperforms.

#### Discussion

We have presented an approach to taking the high-dimensional limit of one-pass
SGD on a least squares problem in which the number of steps is proportional to
the dimension of the problem. The limit object is described in terms of a
Langevin type diffusion, which can be directly compared to the same object in
the multi-pass case.

The literature on scaling limits of one-pass SGD training is large, and so we
mention just some of the closest literature. [AGJ22] is perhaps the closest
high-dimensional diffusion approximation, and it applies in cases where there
is a hidden finite dimensional structure; it covers the case studied here when
$\mathbold{K}=\mathbold{I}$ as well as cases in which $\mathbold{K}$ has
boundedly many eigenvalues. See also [AGJ21].

There are other scaling limits that pursue a different formulation than the
one here. [WL19, WML17] give a PDE limit for the state for a generalized
linear model, with identity covariance. [BGH23] give a scaling limit of the
SGD under smoothness assumptions on the covariance $\mathbold{K},$ when
interpreted geometrically; they further describe fluctuations of SGD in a
certain sense. Note that Theorem 1.1 essentially gives the law of large
numbers for the risks and not the fluctuations.

[Ger+22] uses dynamical field theory to give closely related results for
shallow neural networks with minibatch SGD of large batch-size; dynamical mean
field theory provides an implicit characterization of the autocorrelation of
the minibatch noise and a few other processes. See also the related work of
[CCM21]. We comment that in the case of proportional batch sizes, there is
also a discrete Volterra description in [Lee+22].

#### Organization

In Section 2 we give an overview of the main proof, reducing it to its main
technicalities. In Section 3, we bound the stochastic error terms,
representing the main technical contribution of the paper.

##  2 Main argument of proof

In order to compare the SGD and homogenized SGD, we use a version of the
martingale method in diffusion approximations (see [EK09]). In effect we show
that $q(\mathbold{x}_{k})$ nearly satisfies the conclusion of Itô’s lemma.
Further, we show the martingale terms in both of the Doob decompositions are
small, and hence it suffices to show the predictable parts of
$q(\mathbold{x}_{k})$ and $q(\mathbold{X}_{t})$ are close.

To advance the discussion, we compute this Doob decomposition. To take
advantage of the simpler structure afforded by removing
$\tilde{\mathbold{x}}$, introduce

$\mathbold{v}_{k}\coloneqq\mathbold{x}_{k}-\tilde{\mathbold{x}}\quad\text{and}\quad\mathbold{V}_{t}\coloneqq\mathbold{X}_{t}-\tilde{\mathbold{x}}.$
(2.1)

We shall extend the first integer indexed function to real-valued indices by
setting $\mathbold{v}_{t}=\mathbold{v}_{\lfloor t\rfloor}$. We also let
$(\mathcal{F}_{t}:t\geq 0)$ be the filtration generated by
$(\mathbold{v}_{t}:t\geq 0)$ and $(\mathbold{V}_{t/d}:t\geq 0)$. Hence for all
$k\in\mathbb{N}$, $\mathbold{v}_{k}$ is measurable with respect to
$\mathcal{F}_{k}$. Recalling the recurrence (1.5) for a quadratic $q$

$\begin{split}q(\mathbold{v}_{k})-q(\mathbold{v}_{k-1})&=-\gamma(\nabla
q(\mathbold{v}_{k-1}))^{T}(\mathbold{u}_{k-1}+\Delta_{k})+\tfrac{\gamma^{2}}{2}(\mathbold{u}_{k-1}+\Delta_{k})^{T}(\nabla^{2}q)(\mathbold{u}_{k-1}+\Delta_{k}),\quad\text{where}\quad\\\
\Delta_{k}&=\mathbold{m}_{k}(\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}-\eta_{k})\quad\text{and}\quad\mathbold{u}_{k-1}=\tfrac{\delta}{d}(\mathbold{v}_{k-1}+\tilde{\mathbold{x}}).\end{split}$
(2.2)

The equation above can each be decomposed as a predictable part and two
martingale increments

$\begin{split}q(\mathbold{v}_{k})-q(\mathbold{v}_{k-1})=&-\gamma(\nabla
q(\mathbold{v}_{k-1}))^{T}\bigl{(}(\tfrac{\delta}{d}\mathbold{I}+\tfrac{1}{d}\mathbold{K})\mathbold{v}_{k-1}+\tfrac{\delta\tilde{\mathbold{x}}}{d}\bigr{)}+\Delta\mathcal{M}_{k}^{\operatorname{lin}}\\\
&+\tfrac{\gamma^{2}}{2}\operatorname{tr}(\tfrac{1}{d}\mathbold{K}(\nabla^{2}q))\left(\tfrac{1}{d}\mathbold{v}_{k-1}^{T}\mathbold{K}\mathbold{v}_{k-1}+\mathbb{E}[\eta_{k}^{2}]\right)+\Delta\mathcal{E}_{k}^{\operatorname{quad}}+\Delta\mathcal{M}_{k}^{\operatorname{quad}},\\\
\text{where}\quad\Delta\mathcal{M}_{k}^{\operatorname{quad}}\coloneqq&\Delta_{k}^{T}\nabla^{2}q\Delta_{k}-\mathbb{E}[\Delta_{k}^{T}\nabla^{2}q\Delta_{k}\mid\mathcal{F}_{k-1}].\end{split}$
(2.3)

The remainder of the martingale increments are given by
$\Delta\mathcal{M}_{k}^{\operatorname{lin}}$ and are all linear in
$\Delta_{k}$. The predictable parts have been further decomposed into the
leading order terms and an error term
$\Delta\mathcal{E}_{k}^{\operatorname{quad}}$.

These predictable parts, in turn, depend on different statistics
$q_{1}(\mathbold{v}_{k-1})$. In finite dimensional settings, we would be able
to relate this (or some suitably large finite set of summary statistics
$q,q_{1},\dots,q_{r}$ to itself through a closed system of recurrences. In
this setting, this is not possible. On the other hand, for the problem at
hand, we show there is a manifold of functions which approximately closes.
Specifically, we let

$\begin{split}Q_{n}(q)\coloneqq
Q_{n}(q,\mathbold{K})=&\Big{\\{}q(\mathbold{x}),\quad(\nabla
q(\mathbold{x}))^{T}R(z;\mathbold{K})\mathbold{x},\quad\mathbold{x}^{T}R(y;\mathbold{K})(\nabla^{2}q)R(z;\mathbold{K})\mathbold{x},\\\
&\qquad(\nabla
q(\mathbold{x}))^{T}R(z;\mathbold{K})\tilde{\mathbold{x}},\quad\mathbold{x}^{T}R(y;\mathbold{K})(\nabla^{2}q)R(z;\mathbold{K})\tilde{\mathbold{x}},\quad\forall\,z,y\in\Gamma\Big{\\}}.\end{split}$
(2.4)

Here $R(z;\mathbold{K})=(\mathbold{K}-z\mathbold{I})^{-1}$ is the resolvent
matrix, and $\Gamma$ is a circle of radius $\max\\{1,3\|\mathbold{K}\|\\}$. In
order to control the martingales, it is convenient to impose a stopping time

$\tau\coloneqq\inf\left\\{k:||\mathbold{v}_{k}||>d^{\varepsilon}\\}\cup\\{td:||\mathbold{V}_{t}||>d^{\varepsilon}\right\\},$
(2.5)

and we introduce the corresponding stopped processes

$\mathbold{v}_{k}^{\tau}=\mathbold{v}_{k\wedge\tau},\quad\mathbold{V}_{t}^{\tau}=\mathbold{V}_{t\wedge(\tau/d)}.$
(2.6)

We prove a version of our theorem for the stopped processes and then show that
the stopping time is greater than $n$ with overwhelming probability.

Our key tool for comparing $\mathbold{v}_{td}$ and $\mathbold{V}_{t}$ is the
following lemma.

######  Lemma 2.1.

Given a quadratic $q$ with $\|q\|_{C^{2}}\leq 1$, with $Q=Q_{n}(q)\cup
Q_{n}(\mathscr{P})\cup Q_{n}(\|\cdot\|^{2})$ as above,

$\begin{split}\max_{0\leq
t\leq\tfrac{n}{d}}|q(\mathbold{v}_{td}^{\tau})-q(\mathbold{V}_{t}^{\tau})|\leq&\sup_{0\leq
t\leq{n}/{d}}\left(|\mathcal{M}_{\lfloor
td\rfloor}^{\operatorname{lin},\tau}|+|\mathcal{M}_{\lfloor
td\rfloor}^{\operatorname{quad},\tau}|+|\mathcal{E}_{\lfloor
td\rfloor}^{\operatorname{quad},\tau}|+|\mathcal{M}_{t}^{HSGD,\tau}|\right)\\\
&+C(\|\mathbold{K}\|)\cdot\sup_{g\in
Q}\int_{0}^{n/d}|g(\mathbold{v}_{sd}^{\tau})-g(\mathbold{V}_{s}^{\tau})|ds.\end{split}$
(2.7)

Here $\mathcal{M}_{t}^{HSGD,\tau}$ is the martingale part in the
semimartingale decomposition of $q(\mathbold{V}_{t}^{\tau})$.

######  Proof 2.2 (Sketch of Proof).

Owing to the similarities of this claim with the proof in [Paq+22, Proposition
4.1], we just illustrate the main idea. The idea is that if we take a $g\in
Q$, and we apply (2.3), then in the predictable part of $g(\mathbold{v}_{t})$
we have

$I_{1}\coloneqq\int_{0}^{t}\nabla
g(\mathbold{v}_{sd})^{T}(\delta\mathbold{I}+\mathbold{K})\mathbold{v}_{sd}\mathrm{d}s,\quad
I_{2}\coloneqq\int_{0}^{t}\nabla
g(\mathbold{v}_{sd})^{T}\tilde{\mathbold{x}}\mathrm{d}s,\quad
I_{3}\coloneqq\int_{0}^{t}\mathbold{v}_{sd}^{T}\mathbold{K}\mathbold{v}_{sd}\mathrm{d}s.$

These also appear with coefficients that can be bounded solely using
$\|g\|_{C^{2}}$ and $\|\mathbold{K}\|$. We get the same, applying Itô’s lemma
to $g(\mathbold{V}_{t})$, albeit with the replacement
$\mathbold{v}_{t}\to\mathbold{V}_{t}$. We wish to bound for example
$I_{1}(\mathbold{v}_{t})-I_{1}(\mathbold{V}_{t})$. We do this by expressing
its integrand as $p(\mathbold{v}_{t})-p(\mathbold{V}_{t})$ for polynomial $p$.
If $g$ is linear (the final row of (2.4)), then $p$ is again linear. For
example, if it is $g(\mathbold{x})=\nabla
q(\mathbold{x})^{T}R(z;\mathbold{K})\tilde{\mathbold{x}},$ then $p$ is again
linear and is given by

$p(\mathbold{x})=\mathbold{x}^{T}(\delta\mathbold{I}+\mathbold{K})R(z;\mathbold{K})\tilde{\mathbold{x}}=\delta\mathbold{x}^{T}R(z;\mathbold{K})\tilde{\mathbold{x}}+\mathbold{x}^{T}R(z;\mathbold{K})\tilde{\mathbold{x}}-z\mathbold{x}^{T}\tilde{\mathbold{x}},$

where we have used the resolvent identity
$(\mathbold{K}-z)R(z;\mathbold{K})=\mathbold{I}$. Note the function
$\mathbold{x}^{T}R(z;\mathbold{K})\tilde{\mathbold{x}}$ is contained in $Q$ by
virtue of being in $Q_{n}(\|\cdot\|^{2})$. Moreover, by Cauchy’s integral
formula, we can represent $\mathbold{x}^{T}\tilde{\mathbold{x}}$ by averaging
$\tfrac{-1}{2\pi i}\mathbold{x}^{T}R(y;\mathbold{K})\tilde{\mathbold{x}}$ over
$y\in\Gamma$. Hence

$|p(\mathbold{v}_{td})-p(\mathbold{V}_{t})|\leq(\delta+1+3\|\mathbold{K}\|)\max_{g\in
Q}|g(\mathbold{v}_{td})-g(\mathbold{V}_{t})|.$

The same manipulations lead finally to showing every term included in $Q$ can
be controlled in a similar manner, using the other elements of the class $Q$.

The second important idea is to discretize the set $Q$.

######  Lemma 2.3.

There exists $\bar{Q}\subseteq Q$ with $|\bar{Q}|\leq
C(\|\mathbold{K}\|)d^{4m}$ such that, for every $q\in Q$, there is some
$\bar{q}\in\bar{Q}$ satisfying $||q-\bar{q}||_{C^{2}}\leq d^{-2m}$.

######  Proof 2.4.

On the spectral curve $\Gamma$, we can bound the norm of the resolvent. Since

$\frac{\mathrm{d}}{\mathrm{d}z}R(z;\mathbold{K})=(\mathbold{K}-z\mathbold{I})^{-2},$

we have it is norm bounded by an absolute constant. The arc length of the
curve is at most $C(\|\mathbold{K}\|),$ and so by choosing a minimal net
$d^{-2\varepsilon}$ of the manifold $\Gamma\times\Gamma$, the lemma follows.

Now the main technical part of the argument is to control the martingales and
errors. As we work with the stopped process $\mathbold{v}_{k}^{\tau}$ we
introduce the stopped proccesses
$\mathcal{M}_{k}^{\operatorname{lin},\tau},\mathcal{M}_{k}^{\operatorname{quad},\tau},\mathcal{E}_{k}^{\operatorname{quad},\tau}$,
which are defined analogously to (2.6).

######  Lemma 2.5.

For any quadratic $q$ with $\|q\|_{C^{2}}\leq 1$, the terms
$\mathcal{M}_{k}^{\operatorname{lin},\tau},\mathcal{M}_{k}^{\operatorname{quad},\tau},\mathcal{E}_{k}^{\operatorname{quad},\tau}$
satisfy the following bounds with overwhelming probability (with a bound which
is uniform in $q$) for $n\leq d\log d$

  1. (i)

$\sup_{1\leq k\leq n}|\mathcal{M}_{k}^{\operatorname{lin},\tau}|\leq
d^{-\frac{1}{2}+5\varepsilon}$,

  2. (ii)

$\sup_{1\leq k\leq n}|\mathcal{M}_{k}^{\operatorname{quad},\tau}|\leq
d^{-\frac{1}{2}+9\varepsilon}$,

  3. (iii)

$\sup_{1\leq k\leq n}|\mathcal{E}_{k}^{\operatorname{quad},\tau}|\leq
d^{-1+9\varepsilon}.$

Combining Lemmas 2.1 and 2.3, along with the above (see also Lemma 3.4 in
which the homogenized SGD martingales are bounded), we conclude that, for any
$\bar{q}\in\bar{Q}$ with $\|q\|_{C^{2}}=1$,

$|\bar{q}(\mathbold{v}_{td}^{\tau})-\bar{q}(\mathbold{V}_{t}^{\tau})|\leq
4d^{-\frac{1}{2}+9\varepsilon}+C(\|\mathbold{K}\|)\max_{g\in{Q}}\int_{0}^{t}|g(\mathbold{v}_{sd}^{\tau})-g(\mathbold{V}_{s}^{\tau})|ds.$
(2.8)

Hence by Lemma 2.3 and by bounding $\|g\|_{C^{2}}$ over all $Q$,

$\max_{g\in{Q}}|{q}(\mathbold{v}_{td}^{\tau})-{q}(\mathbold{V}_{t}^{\tau})|\leq
C(\|\mathbold{K}\|)\biggl{(}d^{-2}+d^{-\frac{1}{2}+9\varepsilon}+\int_{0}^{t}\max_{g\in{Q}}|g(\mathbold{v}_{sd}^{\tau})-g(\mathbold{V}_{s}^{\tau})|ds\biggr{)}.$
(2.9)

By Gronwall’s inequality, this gives us that with overwhelming probability

$\max_{g\in{Q}}\max_{0\leq t\leq
n/d}|{g}(\mathbold{v}_{td}^{\tau})-{g}(\mathbold{V}_{t}^{\tau})|\leq
C(\|\mathbold{K}\|)(d^{-2}+4d^{-\frac{1}{2}+9\varepsilon})e^{C(\|\mathbold{K}\|)n/d}.$
(2.10)

Now we note that the norm function $\mathbold{x}\mapsto\|\mathbold{x}\|^{2}$
is one of the quadratics included in $Q.$ Hence if we let $\mathcal{G}$ be the
event in the above display, and we let $\mathcal{E}=\\{\max_{0\leq s\leq
n/d}\|\mathbold{V}_{s}\|\leq d^{\varepsilon/2}\\}$, then we have

$\mathcal{G}\cap\mathcal{E}\cap\\{\tau\leq
n/d\\}\subseteq\\{\|\mathbold{v}_{\tau}\|-\|\mathbold{v}_{\tau-1}\|\geq
d^{\varepsilon/2}\\}\cap\\{\tau\leq n/d\\}.$

This is because on the event $\\{\tau\leq n/d\\}\cap\mathcal{E}$ we must have
had $\|\mathbold{v}_{\tau}\|>d^{\varepsilon}$, but in the step before $\tau,$
we had $\mathbold{v}_{\tau-1}$ could be compared to $\mathbold{V}_{\tau-1}$
(due to $\mathcal{G}$, and we had the norm of $\mathbold{V}_{\tau-1}$ was
small. Now it is easily seen that with overwhelming probability, no increment
of SGD between time $0$ and $n/d$ can increase the norm by a power of $d$. So
to complete the proof it suffices to show $\mathcal{E}$ holds with
overwhelming probability.

Thus the proof is completed by the following:

######  Lemma 2.6.

For any $\delta>0$ and any $t>0$ with overwhelming probability

$\max_{0\leq s\leq t}\|\mathbold{X}_{s}\|^{2}\leq
e^{C(\|\mathbold{K}\|)t}d^{\delta}.$

######  Proof 2.7.

We apply Itô’s formula to
$\phi(\mathbold{X}_{t})\coloneqq\log(1+\|\mathbold{X}_{t}\|^{2})$, from which
we have

$\mathrm{d}\phi(\mathbold{X}_{t})=-2\gamma\tfrac{\mathbold{X}_{t}\cdot\nabla\mathscr{R}(\mathbold{X}_{t})}{1+\|\mathbold{X}_{t}\|^{2}}\mathrm{d}t+\tfrac{\mathbold{X}_{t}\cdot\gamma\sqrt{\tfrac{2}{d}\mathscr{P}(\mathbold{X}_{t})\mathbold{K}}\mathrm{d}B_{t}}{1+\|\mathbold{X}_{t}\|^{2}}+\bigl{(}\tfrac{\mathscr{P}(\mathbold{X}_{t})}{1+\|\mathbold{X}_{t}\|^{2}}\tfrac{2\gamma^{2}}{d}\operatorname{tr}(\mathbold{K})-\tfrac{2\gamma^{2}\mathscr{P}(\mathbold{X}_{t})\mathbold{X}_{t}^{T}\mathbold{K}\mathbold{X}_{t}}{d}\bigr{)}\mathrm{d}t$

The drift terms and the quadratic variation terms can be bounded by some
$C(\|\mathbold{K}\|)$. Hence with this constant, for all $r\geq 0,$

$\mathbb{P}(\max_{0\leq s\leq t}\phi(\mathbold{X}_{s})\geq
C(\|\mathbold{K}\|)(t+r\sqrt{t}))\leq 2\exp(-r^{2}/2).$

Taking $r=\sqrt{\log d\log\log d}$, we conclude that with overwhelming
probability

$\max_{0\leq s\leq t}\phi(\mathbold{X}_{s})\leq
C(\|\mathbold{K}\|)(t+\sqrt{t\log d\log\log d}).$

##  3 Controlling the errors

The main goal of this section is to control the martingale terms and error
terms; in particular we prove Lemma 2.5. In order to obtain these bounds, we
will need the following concentration lemma, which is standard (c.f. [Ver18,
Theorem 2.8.1], where the nonmartingale bound is proven. The adaptation to the
martingale case is a small extension):

######  Lemma 3.1 (Martingale Bernstein inequality).

If $(M_{n})_{1}^{N}$ is a martingale on the filtered probability space
$(\Omega,(\mathcal{F}_{n})_{1}^{N},\mathbb{P}))$ and we define

$\sigma_{n,p}:=\left\|\inf\\{t\geq
0:\mathbb{E}\left(e^{|M_{n}-M_{n-1}|^{p}/t^{p}}|\mathcal{F}_{n-1}\right)\leq
2\\}\right\|_{L^{\infty}(\mathbb{P})},$ (3.1)

then there is an absolute constant $C>0$ so that, for all $t>0$,

$\mathbb{P}\left(\sup_{1\leq n\leq N}|M_{n}-\mathbb{E}M_{0}|\geq t\right)\leq
2\exp\left(-\min\left\\{\frac{t}{C\max\sigma_{n,1}},\frac{t^{2}}{C\sum_{1}^{N}\sigma_{n,1}^{2}}\right\\}\right).$
(3.2)

We will also record for future use an estimate on $\nabla q$ that follows from
$\|\cdot\|_{C^{2}}$ control.

$||\nabla q(\mathbold{x})||\leq||\nabla^{2}q||\cdot||\mathbold{x}||+||\nabla
q(0)||\leq\|q\|_{C^{2}}\cdot(||\mathbold{x}||+1).$ (3.3)

###  3.1 Martingale for gradient part of recurrence

######  Proof 3.2 (Proof of Lemma 2.5 part (i)).

Comparing (2.2) and (2.3), we see that for $k\leq\tau$

$\begin{split}\Delta\mathcal{M}_{k}^{\operatorname{lin},\tau}&=\Big{[}\left(\mathbold{w}_{k-1}^{T}\mathbold{m}_{k}\right)\left(\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}-\eta_{k}\right)-\tfrac{1}{d}\mathbold{w}_{k-1}^{T}\mathbold{K}\mathbold{v}_{k-1}^{\tau}\Big{]}\eqqcolon[\Delta\mathcal{M}_{k}^{\operatorname{lin}1,\tau}-\Delta\mathcal{M}_{k}^{\operatorname{lin}2,\tau}],\\\
&\text{where}\quad\mathbold{w}_{k-1}\coloneqq-\gamma\nabla
q(\mathbold{v}_{k-1}^{\tau})+\tfrac{\gamma^{2}\delta}{d}(\mathbold{v}^{\tau}_{k-1}+\tilde{\mathbold{x}}).\end{split}$
(3.4)

Note for $k>\tau$, the stopped martingale increment is $0.$ Using (3.3),
$\|\mathbold{w}_{k-1}\|\leq C(\gamma,\delta)d^{\varepsilon}.$ We will
separately bound the contributions from
$\Delta\mathcal{M}_{k}^{\operatorname{lin}1,\tau}$ and
$\Delta\mathcal{M}_{k}^{\operatorname{lin}2,\tau}$ in terms of their Orlicz
norms. For the first part, for any fixed $k$, we condition on
$\mathcal{F}_{k-1}$ and Assumption 1, we conclude

$\begin{split}||\Delta\mathcal{M}_{k}^{\operatorname{lin}1,\tau}||_{\psi_{1}}\leq\left\|\mathbold{w}_{k-1}^{T}\mathbold{m}_{k}\right\|_{\psi_{2}}\left\|\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}-\eta_{k}\right\|_{\psi_{2}}\leq
Cd^{-\frac{1}{2}+2\varepsilon}\cdot d^{-\frac{1}{2}+2\varepsilon}\end{split}$
(3.5)

where $C$ is some absolute constant. For the second part, we have

$|\Delta\mathcal{M}_{k}^{\operatorname{lin}2,\tau}|=|\tfrac{1}{d}\mathbold{w}_{k-1}^{T}\mathbold{K}\mathbold{v}_{k-1}^{\tau}|\leq
Cd^{-1+2\varepsilon}.$ (3.6)

Combining these, we see that, for every $k$,

$\sigma_{k,1}:=\inf\\{t>0:\mathbb{E}[\exp(|\Delta\mathcal{M}_{k}^{\operatorname{lin}1,\tau}-\Delta\mathcal{M}_{k}^{\operatorname{lin}2,\tau}|/t)|\mathcal{F}_{k-1}]\leq
2\\}\leq Cd^{-1+4\varepsilon}$ (3.7)

and, by the martingale Bernstein inequality,

$\begin{split}\mathbb{P}\left(\sup_{1\leq k\leq
n}|\mathcal{M}_{k}^{\operatorname{lin},\tau}-\mathbb{E}\mathcal{M}_{0}^{\operatorname{lin}}|\geq
t\right)&\leq
2\exp\left(-\min\left\\{\frac{t}{c\max\sigma_{k,1}},\frac{t^{2}}{c\sum_{k=1}^{n}\sigma_{k,1}}\right\\}\right)\\\
&\leq
2\exp\left(-\min\left\\{Ctd^{1-4\varepsilon},Ct^{2}d^{2-8\varepsilon}n^{-1}\right\\}\right).\end{split}$
(3.8)

As we assume that $n\leq d\log d$ then this gives us

$\sup_{1\leq k\leq n}|\mathcal{M}_{k}^{\operatorname{lin},\tau}|\leq
d^{-\frac{1}{2}+5\varepsilon}$ (3.9)

with overwhelming probability.

###  3.2 Martingale for Hessian part of recurrence

######  Proof 3.3 (Proof of Lemma 2.5 parts (ii) and (iii)).

Next we consider the contribution from the Hessian part of the recurrence. We
write

$\begin{split}&\tfrac{\gamma^{2}}{2}(\mathbold{m}_{k}\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}-\mathbold{m}_{k}\eta_{k})^{T}(\nabla^{2}q)(\mathbold{m}_{k}\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}-\mathbold{m}_{k}\eta_{k})\\\
&=\mathbb{E}\left[\tfrac{\gamma^{2}}{2}(\mathbold{m}_{k}\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}-\mathbold{m}_{k}\eta_{k})^{T}(\nabla^{2}q)(\mathbold{m}_{k}\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}-\mathbold{m}_{k}\eta_{k})|\mathcal{F}_{k-1}\right]+\Delta\mathcal{M}_{k}^{\operatorname{quad}}.\end{split}$
(3.10)

Rearranging the terms, we get

$\Delta\mathcal{M}_{k}^{\operatorname{quad}}=A_{k}B_{k}-\mathbb{E}[A_{k}B_{k}|\mathcal{F}_{k-1}]$
(3.11)

where

$A_{k}:=\mathbold{m}_{k}^{T}(\nabla^{2}q)\mathbold{m}_{k},\quad
B_{k}:=(\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}-\eta_{k})^{2}.$ (3.12)

This can be expanded as

$\begin{split}\Delta\mathcal{M}_{k}^{\operatorname{quad}}=&(A_{k}-\mathbb{E}[A_{k}])(B_{k}-\mathbb{E}[B_{k}])+\mathbb{E}[A_{k}]\mathbb{E}[B_{k}]-\mathbb{E}[A_{k}B_{k}]\\\
&+(A_{k}-\mathbb{E}[A_{k}])\mathbb{E}[B_{k}]+(B_{k}-\mathbb{E}[B_{k}])\mathbb{E}[A_{k}],\end{split}$
(3.13)

so we focus first on obtaining subexponential bounds for the quantities
$A_{k}-\mathbb{E}[A_{k}]$ and $B_{k}-\mathbb{E}[B_{k}]$ using the Hanson-
Wright inequality. For $A_{k}$, we have

$\begin{split}\mathbb{P}(|A_{k}-\mathbb{E}A_{k}|\geq t)&\leq
2\exp\left[-c\min\left(\frac{t^{2}}{d^{-2+4\varepsilon}||\nabla^{2}q||_{HS}^{2}},\frac{t}{d^{-1+2\varepsilon}||\nabla^{2}q||}\right)\right]\\\
&\leq 2\exp[-c^{\prime}\min(t^{2}d^{1-4\varepsilon},td^{1-2\varepsilon})]\leq
2\exp[-c^{\prime\prime}td^{\frac{1}{2}-2\varepsilon}]\end{split}$ (3.14)

and thus we have the subexponential bound

$||A_{k}-\mathbb{E}[A_{k}]||_{\psi_{1}}<Cd^{-\frac{1}{2}+2\varepsilon}.$
(3.15)

Next we obtain a subexponential bound for $B_{k}$. For the part of $B_{k}$ not
involving $\eta_{k}$, we use Hanson-Wright to get

$\begin{split}&\mathbb{P}\left(\left|\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}(\mathbold{v}_{k-1}^{\tau})^{T}\mathbold{m}_{k}-\mathbb{E}\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}(\mathbold{v}_{k-1}^{\tau})^{T}\mathbold{m}_{k}\right|\geq
t\right)\\\ &\leq
2\exp\left[-c\min\left(\frac{t^{2}}{d^{-2+4\varepsilon}||\mathbold{v}_{k-1}^{\tau}(\mathbold{v}_{k-1}^{\tau})^{T}||_{HS}^{2}},\frac{t}{d^{-1+2\varepsilon}||\mathbold{v}_{k-1}^{\tau}(\mathbold{v}_{k-1}^{\tau})^{T}||}\right)\right]\\\
&\leq 2\exp[-c\min(t^{2}d^{2-8\varepsilon},td^{1-4\varepsilon})].\end{split}$
(3.16)

For the terms involving $\eta_{k}$, we use the Orlicz bounds from the
assumptions in the set-up to obtain

$\begin{split}||\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}\eta_{k}||_{\psi_{1}}&\leq||\mathbold{m}_{k}^{T}\mathbold{v}_{k-1}^{\tau}||_{\psi_{2}}\cdot||\eta_{k}||_{\psi_{2}}=d^{-\frac{1}{2}+2\varepsilon}d^{-\frac{1}{2}+\varepsilon}\\\
&=d^{-1+3\varepsilon}.\end{split}$ (3.17)

Since also $||\eta_{k}^{2}||_{\psi_{1}}=d^{-1+2\varepsilon}$ combining the
bounds (3.16) and (3.17), we have

$||B_{k}-\mathbb{E}[B_{k}]||_{\psi_{1}}<Cd^{-1+4\varepsilon}.$ (3.18)

Furthermore, we have

$\mathbb{E}[A_{k}]=O(1),\qquad\mathbb{E}[B_{k}]=O(d^{-1}),$ (3.19)

uniformly for all $k$ based on the assumptions on $\eta_{k}$ and
$\mathbold{m}_{k}$. We now use (3.15), (3.18), (3.19) to bound each term of
(3.13) in turn.

To bound the contribution from
$(A_{k}-\mathbb{E}[A_{k}])(B_{k}-\mathbb{E}[B_{k}])$, we observe that, for
each $k$, with overwhelming probability,
$|A_{k}-\mathbb{E}[A_{k}]|<d^{-\frac{1}{2}+3\varepsilon}$ and
$|B_{k}-\mathbb{E}[B_{k}]|<d^{-1+5\varepsilon}$, so we can conclude that, with
overwhelming probability,

$\sum_{k=1}^{n}\Big{|}(A_{k}-\mathbb{E}[A_{k}])(B_{k}-\mathbb{E}[B_{k}])\Big{|}<nd^{-\frac{3}{2}+8\varepsilon}<d^{-\frac{1}{2}+9\varepsilon}.$
(3.20)

For the second term of (3.13) we have

$\Big{|}\mathbb{E}[A_{k}]\mathbb{E}[B_{k}]-\mathbb{E}[A_{k}B_{k}]\Big{|}=\Big{|}\mathbb{E}\big{[}(A_{k}-\mathbb{E}A_{k})(B_{k}-\mathbb{E}B_{k})\big{]}\Big{|}\leq\mathbb{E}\Big{|}(A_{k}-\mathbb{E}A_{k})(B_{k}-\mathbb{E}B_{k})\Big{|}.$
(3.21)

We can bound this quantity using

$\begin{split}\mathbb{P}\left(\big{|}(A_{k}-\mathbb{E}A_{k})(B_{k}-\mathbb{E}B_{k})\big{|}\geq
t\right)&\leq\mathbb{P}(|A_{k}-\mathbb{E}A_{k}|\geq\sqrt{t})+\mathbb{P}(|B_{k}-\mathbb{E}B_{k}|\geq\sqrt{t})\\\
&\leq
4\exp\left[-c\min(td^{1-4\varepsilon},\sqrt{t}d^{1-4\varepsilon})\right]\end{split}$
(3.22)

where the bound in the last line comes from combining (3.14) and (3.18). Using
this bound, we obtain

$\begin{split}\Big{|}\mathbb{E}[A_{k}]\mathbb{E}[B_{k}]-\mathbb{E}[A_{k}B_{k}]\Big{|}&\leq\int_{0}^{\infty}x\mathbb{P}\left(\big{|}(A_{k}-\mathbb{E}A_{k})(B_{k}-\mathbb{E}B_{k})\big{|}\geq
x\right)dx\\\
&\leq\int_{0}^{1}4x\exp(-cxd^{1-4\varepsilon})dx+\int_{1}^{\infty}4x\exp(-c\sqrt{x}d^{1-4\varepsilon})dx\end{split}$
(3.23)

Making the change of variables $y=xd^{1-4\varepsilon}$ in the first integral
and $z=\sqrt{x}d^{1-4\varepsilon}$ in the second integral, this becomes

$4d^{-2+8\varepsilon}\int_{0}^{d^{1-4\varepsilon}}y\exp(-cy)dy+4d^{-4+8\varepsilon}\int_{d^{1-4\varepsilon}}^{\infty}z^{2}\exp(-cz)dz=O(d^{-2+8\varepsilon}).$
(3.24)

Thus,

$\sum_{k=1}^{n}\Big{|}\mathbb{E}[A_{k}]\mathbb{E}[B_{k}]-\mathbb{E}[A_{k}B_{k}]\Big{|}=O(nd^{-2+8\varepsilon}).$
(3.25)

Finally, we note that the remaining terms of (3.13), namely
$(A_{k}-\mathbb{E}[A_{k}])\mathbb{E}[B_{k}]$ and
$(B_{k}-\mathbb{E}[B_{k}])\mathbb{E}[A_{k}]$, are martingale increments with

$||(A_{k}-\mathbb{E}[A_{k}])\mathbb{E}[B_{k}]||_{\psi_{1}}\leq
Cd^{-\frac{3}{2}+2\varepsilon},\qquad||(B_{k}-\mathbb{E}[B_{k}])\mathbb{E}[A_{k}]||_{\psi_{1}}\leq
Cd^{-1+4\varepsilon}.$ (3.26)

Applying the Martingale Bernstein inequality, we conclude

$\begin{split}\mathbb{P}&\left(\sup_{1\leq k\leq
n}\left|\sum_{j=1}^{k}(A_{j}-\mathbb{E}[A_{j}])\mathbb{E}[B_{j}]+(B_{j}-\mathbb{E}[B_{j}])\mathbb{E}[A_{j}]\right|\geq
t\right)\\\ &\leq
2\exp\left(-\min\left\\{\frac{t}{c\max\sigma_{k,1}},\frac{t^{2}}{c\sum_{k=1}^{n}\sigma_{k,1}}\right\\}\right)\\\
&\leq
2\exp\left(-\min\left\\{Ctd^{1-4\varepsilon},Ct^{2}d^{2-8\varepsilon}n^{-1}\right\\}\right).\end{split}$
(3.27)

Thus, for $n\leq d\log d$, we get

$\sup_{1\leq k\leq
n}\left|\sum_{j=1}^{k}(A_{j}-\mathbb{E}[A_{j}])\mathbb{E}[B_{j}]+(B_{j}-\mathbb{E}[B_{j}])\mathbb{E}[A_{j}]\right|\leq
d^{-\frac{1}{2}+5\varepsilon}$ (3.28)

with overwhelming probability. Finally, combining the bounds from (3.20),
(3.25), (3.28), we conclude that, for $n\leq d\log d$,

$\sup_{1\leq k\leq n}|\mathcal{M}_{k}^{\operatorname{quad},\tau}|\leq
d^{-\frac{1}{2}+8\varepsilon}$ (3.29)

with overwhelming probability. This completes the proof of part (ii) of the
lemma.

For part (iii), we observe that
$\Delta\mathcal{E}_{k}^{\operatorname{quad},\tau}=\mathbb{E}[A_{k}B_{k}]-\mathbb{E}[A_{k}]\mathbb{E}[B_{k}]+O(d^{-2+4\epsilon})$,
the error terms arising from $\mathbold{u}_{k}$ cross terms, so that the bound
of $\mathcal{E}_{k}^{\operatorname{quad},\tau}$ follows immediately from
(3.25).

###  3.3 Martingale for HSGD

Recall that the HSGD process $\\{\mathbold{V}_{t}\\}$ satisfies the
differential equation

$d\mathbold{V}_{t}=-\gamma\nabla\mathscr{R}(\mathbold{V}_{t}+\tilde{\mathbold{x}})dt+\gamma\sqrt{\tfrac{2}{d}\mathscr{P}(\mathbold{V}_{t}+\tilde{\mathbold{x}})\mathbold{K}}dB_{t}$
(3.30)

where

$\mathscr{P}(\mathbold{x})=(\mathbold{x}-\tilde{\mathbold{x}})^{T}\mathbold{K}^{2}(\mathbold{x}-\tilde{\mathbold{x}})+\eta^{2}.$
(3.31)

Using Itô’s Lemma, this gives us

$q(\mathbold{V}_{t}^{\tau})=q(\mathbold{V}_{0}^{\tau})-\gamma\int_{0}^{t}(\nabla
q(\mathbold{V}_{s}^{\tau}))^{T}\nabla\mathscr{R}(\mathbold{V}_{s}^{\tau}+\tilde{\mathbold{x}})ds+\tfrac{\gamma^{2}}{2}\int_{0}^{t}\mathscr{P}(\mathbold{V}_{s}^{\tau}+\tilde{\mathbold{x}})\operatorname{tr}(\tfrac{1}{d}\mathbold{K}(\nabla^{2}q))ds+\mathcal{M}_{t}^{HSGD,\tau},$
(3.32)

where

$\mathcal{M}_{t}^{HSGD,\tau}=\gamma\int_{0}^{t}(\nabla
q(\mathbold{V}_{s}^{\tau}))^{T}\sqrt{\tfrac{1}{d}\mathscr{P}(\mathbold{V}_{s}^{\tau}+\tilde{\mathbold{x}})\mathbold{K}}dB_{s}.$
(3.33)

######  Lemma 3.4.

For any quadratic $q$ with $\|q\|_{C^{2}}\leq 1$, after imposing the stopping
time $\tau$, the resulting martingale $\mathcal{M}^{HSGD,\tau}_{t}$ satisfies

$\sup_{0\leq t\leq n/d}\left[\mathcal{M}^{HSGD,\tau}_{t}\right]\leq
Cd^{-\frac{1}{2}+3\varepsilon},$ (3.34)

provided $n\leq d\log d$.

######  Proof 3.5.

This martingale has quadratic variation

$\left[\mathcal{M}^{HSGD,\tau}_{t}\right]=\frac{\gamma^{2}}{d}\int_{0}^{t}\mathscr{P}(\mathbold{V}_{s}^{\tau}+\tilde{\mathbold{x}})(\nabla
q(\mathbold{V}_{s}^{\tau}))^{T}\nabla^{2}\mathscr{P}(\mathbold{V}_{s}^{\tau}+\tilde{\mathbold{x}})(\nabla
q(\mathbold{V}_{s}^{\tau}))ds$ (3.35)

and, for all $s$, we have the bound

$\begin{split}&|\mathscr{P}(\mathbold{V}_{s}^{\tau}+\tilde{\mathbold{x}})|\leq||\mathbold{V}_{s}^{\tau}||^{2}||\mathbold{K}^{2}||+\eta^{2}\leq
C_{1}d^{2\varepsilon},\\\ \end{split}$ (3.36)

Thus, we have

$\sup_{0\leq t\leq n/d}\left[\mathcal{M}^{HSGD,\tau}_{t}\right]\leq C(\log
d)d^{-1+4\varepsilon}$ (3.37)

almost surely. By the Gaussian tail bound for continuous martingales of
bounded quadratic varation,
$\mathbb{P}(\sup_{t}[\mathcal{M}^{HSGD,\tau}_{t}]>a)<\exp(-a^{2}/(2C(\log
d)d^{-1+4\varepsilon}))$. Taking $t=d^{-\frac{1}{2}+3\varepsilon}$, we
conclude that

$\sup_{0\leq t\leq n/d}\left[\mathcal{M}^{HSGD,\tau}_{t}\right]\leq
Cd^{-\frac{1}{2}+3\varepsilon}$ (3.38)

with overwhelming probability.

## References

  * [Ada15] R. Adamczak  “A note on the Hanson-Wright inequality for random vectors with dependencies”  In _Electron. Commun. Probab._ 20, 2015, pp. no. 7213  DOI: 10.1214/ECP.v20-3829
  * [AGJ21] Gerard Ben Arous, Reza Gheissari and Aukosh Jagannath  “Online stochastic gradient descent on non-convex losses from high-dimensional inference”  In _The Journal of Machine Learning Research_ 22.1 JMLRORG, 2021, pp. 4788–4838
  * [AGJ22] Gerard Ben Arous, Reza Gheissari and Aukosh Jagannath  “High-dimensional limit theorems for SGD: Effective dynamics and critical scaling”  In _Advances in Neural Information Processing Systems_ , 2022  URL: https://openreview.net/forum?id=Q38D6xxrKHe
  * [AYS20] Kwangjun Ahn, Chulhee Yun and Suvrit Sra  “SGD with shuffling: optimal rates without component convexity and large epoch requirements”  In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 17526–17535  URL: https://proceedings.neurips.cc/paper_files/paper/2020/file/cb8acb1dc9821bf74e6ca9068032d623-Paper.pdf
  * [BGH23] Krishnakumar Balasubramanian, Promit Ghosal and Ye He  “High-dimensional scaling limits and fluctuations of online least-squares SGD with smooth covariance”  In _arXiv e-prints_ , 2023, pp. arXiv:2304.00707  DOI: 10.48550/arXiv.2304.00707
  * [CCM21] Michael Celentano, Chen Cheng and Andrea Montanari  “The high-dimensional asymptotics of first order methods with random data”  In _arXiv e-prints_ , 2021, pp. arXiv:2112.07572  DOI: 10.48550/arXiv.2112.07572
  * [EK09] Stewart N Ethier and Thomas G Kurtz  “Markov processes: characterization and convergence”  John Wiley & Sons, 2009
  * [Ger+22] Cedric Gerbelot, Emanuele Troiani, Francesca Mignacco, Florent Krzakala and Lenka Zdeborova  “Rigorous dynamical mean field theory for stochastic gradient descent methods”  In _arXiv e-prints_ , 2022, pp. arXiv:2210.06591  DOI: 10.48550/arXiv.2210.06591
  * [GOP21] Mert Gürbüzbalaban, Asu Ozdaglar and Pablo A Parrilo  “Why random reshuffling beats stochastic gradient descent”  In _Mathematical Programming_ 186 Springer, 2021, pp. 49–84
  * [Lee+22] Kiwon Lee, Andrew N. Cheng, Courtney Paquette and Elliot Paquette  “Trajectory of Mini-Batch Momentum: Batch Size Saturation and Convergence in High Dimensions”  In _To Appear in NeurIPS 2022_ , 2022, pp. 38pp  arXiv:2206.01029 [math.OC]
  * [Paq+21] Courtney Paquette, Kiwon Lee, Fabian Pedregosa and Elliot Paquette  “SGD in the Large: Average-case Analysis, Asymptotics, and Stepsize Criticality”  In _Proceedings of Thirty Fourth Conference on Learning Theory_ 134, Proceedings of Machine Learning Research  PMLR, 2021, pp. 3548–3626  arXiv:2102.04396 [math.OC]
  * [Paq+22] Courtney Paquette, Elliot Paquette, Ben Adlam and Jeffrey Pennington  “Homogenization of SGD in high-dimensions: Exact dynamics and generalization properties”  In _arXiv e-prints_ , 2022, pp. 64pp  arXiv:2205.07069 [math.ST]
  * [RM51] H. Robbins and S. Monro  “A Stochastic Approximation Method”  In _Ann. Math. Statist._ The Institute of Mathematical Statistics, 1951  DOI: 10.1214/aoms/1177729586
  * [RR08] A. Rahimi and B. Recht  “Random features for large-scale kernel machines”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 20, 2008, pp. 1177–1184
  * [RR12] Benjamin Recht and Christopher Re  “Toward a Noncommutative Arithmetic-geometric Mean Inequality: Conjectures, Case-studies, and Consequences”  In _Proceedings of the 25th Annual Conference on Learning Theory_ 23, Proceedings of Machine Learning Research  Edinburgh, Scotland: PMLR, 2012, pp. 11.1–11.24  URL: https://proceedings.mlr.press/v23/recht12.html
  * [SS20] Itay Safran and Ohad Shamir  “How Good is SGD with Random Shuffling?”  In _Proceedings of Thirty Third Conference on Learning Theory_ 125, Proceedings of Machine Learning Research  PMLR, 2020, pp. 3250–3284  URL: https://proceedings.mlr.press/v125/safran20a.html
  * [Ver18] Roman Vershynin  “High-dimensional probability: An introduction with applications in data science”  Cambridge university press, 2018
  * [WL19] Chuang Wang and Yue M Lu  “The scaling limit of high-dimensional online independent component analysis*”  In _Journal of Statistical Mechanics: Theory and Experiment_ 2019.12 IOP PublishingSISSA, 2019, pp. 124011  DOI: 10.1088/1742-5468/ab39d6
  * [WML17] Chuang Wang, Jonathan Mattingly and Yue M. Lu  “Scaling Limit: Exact and Tractable Analysis of Online Learning Algorithms with Applications to Regularized Regression and PCA”  In _arXiv e-prints_ , 2017, pp. arXiv:1712.04332  DOI: 10.48550/arXiv.1712.04332
  * [YSJ21] Chulhee Yun, Suvrit Sra and Ali Jadbabaie  “Open problem: Can single-shuffle SGD be better than reshuffling SGD and gd?”  In _Conference on Learning Theory_ , 2021, pp. 4653–4658  PMLR


# Steady periodic irrotational blood flow with time-dependent body force111
This paper was jointly supported from Natural Science Foundation of Zhejiang
Province (No.LZ23A010001), National Natural Science Foundation of China (No.
11671176, 11931016).

Yuchao Hea Yongli Songb Yonghui Xiaa222Corresponding author.
yhxia<EMAIL_ADDRESS>
a School of Mathematical Science, Zhejiang Normal University, 321004, Jinhua,
China
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
b School of Mathematics, Hangzhou Normal University, 311121, Hangzhou, China
Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we study steady 2D periodic blood flow propagating along blood
vessels with free boundary conditions. In particular, we focus on the dynamic
behavior of irrotational flows with time dependent body force. An equivalent
formulation with fix boundary is developed by utilizing flow force functions.
The local bifurcation result is obtained by Crandall-Rabinowitz theorem. The
existence of a local $C^{1}$-curve of small-amplitude solution is strictly
proved.
Keywords: Blood flow; Free boundary; Local bifurcation

MSC: 76D33; 76B15; 35J60; 47J15

##  1 Introduction

The study of human arterial blood flow model has a long history, it tracks
back to the original work by Euler [17] in 1775. He proposed a one-dimension
inviscid flow with the conservation of mass and momentum. Euler’s work was
innovative, but regrettably failed in considering the wave-like nature of the
flow and finding a solution for his equations. Later, in Young [37], the wave
nature of fluids was incorporated into Euler’s model. After 100 years later,
Anliker et al. [27] and Skalak [28] applied Riemann’s characteristics method
(see [26]) to develop the research on the artistic flow.

Up till to 21-th century, the study of the blood flows has once again
attracted widespread attentions. In 2002, Quarteroni and Formaggia [31]
studied the problem of haemodynamics and 2D, 3D blood flow models. For the
convenience of research, the model of blood flow was often seen as an
incompressible non-Newtonian fluid with viscosity. In 2010, Dong [16] studied
the Cauchy problem to the viscous incompressible micropolar non-Newtonian
fluid flows. And Dong concerned with the asymmetric stress tensor in the
system. In 2015, Friedman [19] et al. considered a free boundary problem for
small plaques in the artery and studied the stability of blood flow system
with some conditions. In 2016, Quarteroni et al. [30] use the geometric
multiscale approach to deal with the 3D blood flow problem and emphasize the
need for rigorous mathematical models to fill the gap between theory and
practice. In order to better simulate the dynamic behavior in actual blood
flow, Tabi et al. [32] considered the impact of viscosity, radiations,
magnetic fields on velocity and temperature on the blood flow. And in 2021, Bi
et al. [4] studied one kind of incompressible blood flow under the weak
viscosity condition. A traveling wave solution and some properties of it are
given. The numerical algorithms on blood flow model have also attracted a
large number of scholars to conduct research, one can refer to [1, 20, 25, 29,
36]. In the systematic monograph on the mathematical modelling of the
cardiovascular system [31], the authors mentioned that a few of factors should
be incorporated into the blood vessel artery based on the basic model. These
factors includes axial symmetry, radial displacements, constant pressure,
dominance of axial velocity and body forces. Quarteroni and Formaggia [31]
incorporated the body force into the blood model described by a viscous and
incompressible Navier-Stokes equation as follows:

$\frac{\partial u}{\partial t}+(u\cdot\nabla)u+\nabla p-\mu\Delta u=f,$

where $\mu$ is the viscosity coefficient, $p$ is the pressure function and $f$
is the body force. For more detailed applications of the body force, one can
refer to [18, 23].

However, most of the above mentioned works are on the numerical computation of
the blood flow model. The problem of solving the dynamic equations of blood
flow with boundary conditions is rarely solved by scholars. It is known that
if the surface of flow is considered as a free boundary (see Constantin [6]),
solving the PDEs becomes particularly complex. An effective way to describe
the fluid dynamics is to find the exact nonlinear geophysical Gerstner-like
solutions (see Constantin [9]) and study qualitative behaviors (see [11, 13,
33, 34, 35]) such as bifurcation and stability, develop a flow force
reformulation (see [5, 7, 8, 10, 24, 22]) and find laminar solutions (see [2,
3, 30, 21]). In this paper, we use a generalized Euler control equation with
viscosity terms and free boundary conditions to describe the dynamic behavior
of blood flow in arteries. In particular, we note that the volume force in the
blood flow is not constant, but changes over time (For more details on the
time-dependent body force, readers can refer to [23, 18]). We pay our
particular attentions to the irrotational blood flow. Constantin and Strauss
[6], Constantin et al. [10], Constantin and Escher [con1] made laying a
foundation contribution to the works of the rotational water waves. Basu [2],
Chu [12], Henry [21] made contributions to the irrotational water waves. In
this paper, we focus ourselves on the irrotational blood flow with time-
dependent body force. By introducing the flow force function and height
transformation, we obtain a class of equivalent equations and corresponding
exact laminar solutions. Further, according to the bifurcation theorem in
[15], local bifurcation results have been presented. Specially, we prove the
existence of a local $C^{1}$-curve of small-amplitude solution. In addition,
we find that if the radius of the blood vessel is small enough, the small-
amplitude solutions exist. As far as we know, most of the work on the dynamic
behavior of the blood flow are numerical (e.g. [1, 29, 20, 36, 25]). Different
from their works on the numerical computation, motivated by the works of Basu
[2], Chu [12], Henry [21], we develop a flow force reformulation and find the
exact solutions, which is interesting from the point of mathematics. This
helps us easier to observe the dynamics of the blood flow. Consequently, our
method sheds some new light on the study of the blood flow. Our flow force
formulation of the blood flow would benefit to design numerical algorithms for
the blood flow model.

The outline of the rest of this paper: in Section 2, we state a blood flow
model with time-varying body forces, and the problem is equivalent to a second
order quasi-linear elliptic equation by some transformations. Finally, in
Section 3, the existence of a local $C^{1}$-curve of small-amplitude solution
is strictly proved by Crandall-Rabinowitz bifurcation theorem.

##  2 Formulation of blood flow model and flow force reformulation

###  2.1 Formulation of blood flow model

In the vascular system, choose $(x,y)\in\Phi\subset\mathbb{R}^{2}$ and let the
$x$-axis be the axial direction of blood vessels and the $y$-axis be the
radial direction of blood vessels. In actual vascular fluid models, there are
usually two types of situations: the blood vessels are fully filled with
blood; the blood vessels are not fully filled with blood.

If the blood vessels are not fully filled with blood, the free surface,
denoted by $y=\Omega(x,t)$, on the domain bounded above satisfies

$\int_{-\pi}^{\pi}\Omega(x,t)dx=0.$ (2.1)

And below by the central axis of blood vessels $y=-d$ with $0<d<\infty$, the
period respect to $x$ is chosen as $2\pi$. The situation when the blood flow
fully fills the blood vessels can be seen as a special case of the above
situation and $\Omega(x,t)$ reduces to $0$. Without loss of generality, we
only consider the case when the blood vessel is not fully filled in the
following (see figure 1).

Figure 1: Cross section of blood flow in blood vessels.

Assume that the speed $c$ of traveling waves is positive ($c>0$). That is, the
free surface $\Omega(x,t)$ takes the form of $\Omega(x,t)=\Omega(x-ct)$.

By traveling wave transformation, we transform $(x-ct,y)$ into $(x,y)$. And
$(u(x,y),v(x,y))$ is used to denote the blood velocity field in the vascular
region $\\{(x,y)\in R^{2}:-d<y<\Omega(x)\\}$. Set that the blood flow has no
stagnation points, i.e.

$u-c<0\quad for\quad-d\leq y\leq\Omega(x),$ (2.2)

through blood vessels.
In this paper, we formulate a blood flow model described by the following
modified Navier-Stokes equations:

$\begin{cases}u_{t}+uu_{x}+vu_{y}+\mu\Delta
u=-\frac{1}{\rho}P_{x}+f_{1}(x-ct,y),\\\ v_{t}+uv_{x}+vv_{y}+\mu\Delta
v=-\frac{1}{\rho}P_{y}-f_{2}(x-ct,y),\end{cases}$ (2.3)

with the condition of mass conservation

$u_{x}+v_{y}=0.$ (2.4)

The movement of blood flow satisfies the following boundary conditions:

$v=0\quad\rm{on}\ \ y=-d,$ (2.5)

and

$P=P_{0}\quad on\ \ y=\Omega(x).$ (2.6)

If the blood flow does not fully fill the blood vessels, $v$ satisfies

$v=(u-c)\Omega^{\prime}(x)\quad\rm{on}\ \ y=\Omega(x).$ (2.7)

Otherwise,

$v=0\quad\rm{on}\ \ y=\Omega(x).$ (2.8)

In the above equations, $\mu$ is the viscosity coefficient, $\rho$ is the
constant density, $P$ is the pressure, $f_{1}$ and $f_{2}$ are non-negative
body forces satisfying the irrotational condition

$\frac{\partial f_{1}}{\partial y}=-\frac{\partial f_{2}}{\partial x}.$ (2.9)

And we assume that for any $x$, the body force $f_{1}$ in vertical direction
satisfies

$\int_{-d}^{\Omega(x)}f_{1}(x,y)dy=0,$

$f_{1}(x,y)$ vanishes on the boundary: $f_{1}(x,\Omega(x))=f_{1}(x,-d)=0.$
We defined the stream function $\Phi(x,y)$: $\Phi\to\mathbb{R}$ by as follows

$\begin{cases}\Phi_{x}=-v,\\\ \Phi_{y}=u-c.\end{cases}$ (2.10)

Moreover, $\Phi$ is normalized by choosing $\Phi=0$ on $y=\Omega(x)$ and then
$\Phi=-p_{0}$ on $y=-d$, where

$p_{0}=\int_{-d}^{\Omega(x)}(u(x,y)-c)dy.$

It is easy to see that $p_{0}$ is a constant by taking the derivative with
respect to $x$. $P_{0}$ varies with the solution of (2.3)-(2.8). Vorticity
(see [5]) is an important factor in characterizing fluid dynamics behavior. We
are more concerned with the dynamic behavior of the flow in the irrotational
condition, that is

$u_{y}-v_{x}=0.$ (2.11)

###  2.2 Equivalent formulations

No-rotational condition (2.11) plays a key role in later computations,
especially in the conversion of equivalent formulations. Combining (2.11) with
(2.4), we obtain $u$, $v$ are all harmonic functions. In fact,

$\Delta u=u_{xx}+u_{yy}=u_{xx}+v_{xy}=(u_{x}+v_{y})_{x}=0.$

This allows the effect of viscosity on the model to be ignored. By taking the
traveling wave transform $(x-ct,y)\to(x,y)$, (2.3) finally reduces to

$\begin{cases}(u-c)u_{x}+vu_{y}=-\frac{1}{\rho}P_{x}+f_{1}(x-ct,y),\\\
(u-c)v_{x}+vv_{y}=-\frac{1}{\rho}P_{y}-f_{2}(x-ct,y),\end{cases}$ (2.12)

Distinctly, according to (2.12), the total energy

$E=\frac{\Phi_{x}^{2}+\Phi_{y}^{2}}{2}+\frac{P}{\rho}+F(x,y)$ (2.13)

is a constant throughout the blood fluids, where $F_{x}=-f_{1},\quad
F_{y}=f_{2}$. The existence of $F(x,y)$ follows form (2.9). Based on the
analysis above, we obtain that equations (2.4)-(2.8) and (2.12) are
transformed into the following equations

$\begin{cases}\Delta\Phi=0\quad\ -d<y<\Omega(x),\\\
|\nabla\Phi|^{2}+2F(x,y)=Q\quad\ on\ y=\Omega(x),\\\ \Phi=0\quad\ on\
y=\Omega(x),\\\ \Phi=-p_{0}\quad\ on\ y=0,\end{cases}$ (2.14)

where $Q=2(E-\frac{P_{0}}{\rho})$ is constant i.e. $\frac{\partial Q}{\partial
x}=\frac{\partial Q}{\partial y}=0.$

###  2.3 Flow force function

We now introduce the flow force function

$z(x,y)=\int_{-d}^{\Omega(x)}[\frac{p(x,r)}{\rho}+(u(x,r)-c)^{2}]dr,$ (2.15)

where $p(x,r)=P(x,r)-P_{0}$. Visibly, $z=0$ on the axis of blood vessels
$y=-d$. Denote $z_{0}=z|_{y=\Omega(x)}$. Moreover,

$z_{0}=\int_{-d}^{\Omega(x)}[\frac{p(x,y)}{\rho}+(u(x,y)-c)^{2}]dy.$

Take the derivative of $z$ with respect to $x$, then

$\begin{split}\frac{dz_{0}}{dx}&=\left(\frac{p(x,\Omega(x))}{\rho}+(u(x,\Omega(x)))^{2}\right)\Omega^{\prime}(x)+\int_{-d}^{\Omega(x)}\left(\frac{p_{x}}{\rho}+2(u-c)u_{x}\right)dy\\\
&=(u(x,\Omega(x))-c)((u-c)\Omega^{\prime}(x)-v)+\int_{-d}^{\Omega(x)}f_{1}(x,y)dy\\\
&=0.\end{split}$ (2.16)

And (2.16) means that $z_{0}$ is a constant.

Take $w=x$, $z=z(x,y)$. It is worth mentioning that under the new variables,
the original free boundary blood flow region is transformed into a new fixed
boundary region, greatly simplifying the problem. And it is intuitive in the
graph 2.

Figure 2: Sketch of the blood flow region transformation.

Note that

$\begin{split}\frac{\partial}{\partial x}&=\frac{\partial}{\partial
w}+\frac{\partial z}{\partial x}\cdot\frac{\partial}{\partial z},\\\
\frac{\partial}{\partial y}&=\frac{\partial}{\partial w}\cdot\frac{\partial
w}{\partial y}+\frac{\partial}{\partial z}\cdot\frac{\partial z}{\partial
y},\end{split}$

then

$\begin{split}z_{x}&=v(c-u)+\int_{-d}^{y}f_{1}(x,r)dr,\\\
z_{y}&=\frac{p(x,y)}{\rho}+(u-c)^{2}.\end{split}$

And it is easy to see that $z_{y}>0.$
Moreover,

$\begin{split}z_{xx}&=v_{x}(u-c)+vu_{x}+\int_{0}^{y}\frac{\partial
f_{1}(x,r)}{\partial x}dr,\\\
z_{yy}&=(u-c)u_{y}-vv_{y}-f_{2}(x,y).\end{split}$

According to (2.4) and (2.11), we obtain

$\Delta z=\int_{-d}^{y}\frac{\partial f_{1}(x,r)}{\partial x}dr-
f_{2}(x,y),\quad\ 0<y<\Omega(x).$ (2.17)

According to (2.9), $\int_{-d}^{y}\frac{\partial f_{1}(x,r)}{\partial x}dr-
f_{2}(x,y)$ in (2.17) is reduced to a constant $-c_{1}$. And we assume
$c_{1}>0$ in this paper.
In order to convert a free boundary problem into a fixed boundary problem, the
height function in the new $(w,z)$-variables is defined as follows,

$h(w,z)=y(w,z)+d.$

Even and $2\pi$-periodic in $w$ solutions are sought in the following. The
condition (2.7) turns to be

$\int_{-\pi}^{\pi}h(w,z_{0})dw=2\pi d.$ (2.18)

Motivated by Basu [2], the flow force function is invariant along a ‘flow
force line’, that is

$z_{x}+z_{y}h_{w}=0.$ (2.19)

Moreover,

$\begin{split}h_{w}&=-\frac{z_{x}}{z_{y}},\\\
h_{z}&=\frac{1}{z_{y}}.\end{split}$ (2.20)

Then, we obtain

$z_{x}=v(c-u)+\int_{-d}^{y}f_{1}(x,r)dr=-\frac{h_{w}}{h_{z}}$ (2.21)

and

$z_{y}=\frac{p(x,y)}{\rho}+(u-c)^{2}=\frac{1}{h_{z}}.$ (2.22)

It follows from (2.21) and (2.22),

$\begin{split}z_{xx}=v_{x}(c-u)-vu_{x}+\int_{-d}^{y}\frac{\partial
f_{1}(x,r)}{\partial x}dr&=(\frac{\partial}{\partial
w}-\frac{h_{w}}{h_{z}}\cdot\frac{\partial}{\partial
z})(-\frac{h_{w}}{h_{z}}),\\\
z_{yy}=\frac{p_{y}(x,y)}{\rho}+2(u-c)u_{y}&=(\frac{1}{h_{z}}\cdot\frac{\partial}{\partial
z})\frac{1}{h_{z}}.\end{split}$ (2.23)

Combining with (2.17) and (2.23), we obtain

$-c_{1}=\int_{-d}^{y}\frac{\partial f_{1}(x,r)}{\partial x}dr-
f_{2}(x,y)=-(\frac{\partial}{\partial
w}-\frac{h_{w}}{h_{z}}\frac{\partial}{\partial
z})(\frac{h_{w}}{h_{z}})+\frac{1}{h_{z}}\frac{\partial}{\partial
z}(\frac{1}{h_{z}}),$

i.e.

$(1+h_{w}^{2})h_{zz}-2h_{z}h_{w}h_{zw}+h_{z}^{2}h_{ww}=c_{1}h_{z}^{3}.$ (2.24)

Moreover, the second equation in (2.14) is transformed into

$(1+h_{w}^{2})=\left[Q-2F(w,z_{0})\right]h_{z},\quad\ \rm{on}\quad z=z_{0}$

and

$h=0,\quad\rm{on}\quad z=0.$

Thus, the previous free boundary problem is equivalent to the equations of $h$

$\begin{cases}(1+h_{w}^{2})h_{zz}-2h_{z}h_{w}h_{wz}+h_{z}^{2}h_{ww}=c_{1}h_{z}^{3},\quad\
0<z<z_{0},\\\ 1+h_{w}^{2}=\left[Q-2F(w,z_{0})\right]h_{z},\quad\ \rm{on}\quad
z=z_{0},\\\ h=0,\quad\ on\quad z=0,\\\
\int_{-\pi}^{\pi}h(w,z_{0})dw=d.\end{cases}$ (2.25)

By resemble statements as in [21], (2.25) is equivalent to (2.12).

##  3 Main results

###  3.1 Laminar solution

Note that, if $z=z_{0}$, then $y=\Omega(x)$ and $\frac{\rm
dF(x,\Omega(x))}{\rm dx}=f_{1}(x,\Omega(x))=0$, $F(x,y)$ is independent of $w$
on $y=\Omega(x)$ or $z=z_{0}$. Set $F_{1}(y)=F(w,z_{0})$, moreover
$F_{1}(y)|_{y=\Omega(x)}=F_{1}(H(1)-d)$. Since $F$ is the first integral, it
can be corrected by adding or subtracting a constant. Without loss generality,
we let $F(0,-d)=0$, or else, if $F(0,-d)=constant\neq 0$, take
$\widehat{F}=F-constant.$ Firstly, we seek the laminar solutions $H(z)$ to
(2.25), which is independent of $w$. Consequently, $H(z)$ solves

$\begin{cases}H_{zz}=c_{1}H_{z}^{3},\quad\ 0<z<z_{0},\\\
1=[Q-2F_{1}(H(z_{0})-d)]H_{z}(z_{0}),\\\ H(0)=0,\\\ H(z_{0})=d.\end{cases}$
(3.1)

By direct calculation, the first equation in (3.1) has the solution

$H_{z}=[\xi^{2}+2c_{1}(z_{0}-z)]^{-\frac{1}{2}},$ (3.2)

where $\xi=\frac{1}{H_{z}(z_{0})}>0.$
The second equation of (3.1) is transformed into

$H(z)=-\frac{1}{c_{1}}[\xi^{2}+2c_{1}(z_{0}-z)]^{\frac{1}{2}}+F^{-1}(\frac{Q-\xi}{2})+d+\frac{\xi}{c_{1}},$
(3.3)

where the existence of $F_{1}^{-1}$ follows from the nonnegativity of
$f_{1},f_{2}.$ Combining with $H(0)=0$ and $H(z_{0})=d$, we obtain

$\begin{cases}-\frac{1}{g}(\xi^{2}+2gz_{0})^{\frac{1}{2}}+F_{1}^{-1}(\frac{Q-\xi}{2})+d+\frac{\xi}{g}=0,\\\
F_{1}^{-1}(\frac{Q-\xi}{2})=-d.\end{cases}$ (3.4)

From (3.4), we obtain

$z_{0}=\frac{d(dc_{1}+2\xi)}{2}$ (3.5)

and

$Q=\xi.$ (3.6)

In this part, we seek solutions $h\in C^{3,\alpha}_{per}(\bar{R}_{z})$ even in
the $w$-variable and zero on $z=0.$ In the following part of this paper, we
restrict $\xi>0.$
Then

$H(z,\xi)=-\frac{1}{c_{1}}[\xi^{2}+2c_{1}(\frac{d(dc_{1}+2\xi)}{2}-z)]^{\frac{1}{2}}+d+\frac{\xi}{c_{1}}$

holds for $\xi\in(0,+\infty).$

###  3.2 Linearised problem

To deal with the linearised problem, we decompose $h(w,z)$ into two parts,
that is

$h(w,z)=H(z,\xi)+\epsilon m(w,z).$ (3.7)

where $m\in C^{3,\alpha}_{per}(\bar{R}_{z})$ is even in the $w$-variable. Note
that the Taylor expansion of $F|_{z=z_{0}}$ at $\epsilon=0$, we obtain

$\begin{split}F(w,h-d)&=F(0,-d)+\int_{-d}^{h-d}f_{2}(w,r)dr-\int_{0}^{w}f_{1}(r,-d)dr\\\
&=\int_{-d}^{\epsilon m}f_{2}(w,r)dr\\\
&=\int_{-d}^{0}f_{2}(w,r)dr+f_{2}(w,0)m\epsilon+o(\epsilon^{2}).\end{split}$

Therefore,

$1+(H_{w}+\epsilon
m_{w})^{2}=\left[Q-2\left(\int_{-d}^{0}f_{2}(x,r)dr+\epsilon
mf_{2}(x,0)\right)\right](H_{z}+\epsilon m_{z})\quad\quad on\quad z=z_{0}.$
(3.8)

We only focus on the first order of (2.25) in $\epsilon$. According to (2.25)
and (3.7), for $m$ the boundary problem is on the form:

$\begin{cases}(a^{3}m_{z})_{z}+am_{ww}=0,\quad\ in\ \mathbb{R}_{z}\\\
\xi\left(\frac{Q}{2}-\int_{-d}^{0}f_{2}(x,r)dr\right)m_{z}=f_{2}(x,0)m,\quad\
z=z_{0},\\\ m=0,\quad\ z=0,\end{cases}$ (3.9)

where

$a=\sqrt{\xi^{2}+2c_{1}(z_{0}-z)},\quad\xi\in(0,z_{0}).$

Note that $m$ is an even function, we have the following equation by using
Fourier transform.

$m(w,z)=\sum_{k=0}^{\infty}m_{k}(z)cos(kw),\ in\
C^{2}_{per}(\widetilde{R}_{z})$

with the coefficients

$m_{0}(z)=\frac{1}{2\pi}\int_{-\pi}^{\pi}m(w,z)dz,\quad\
m_{k}(z)=\frac{1}{2\pi}\int_{-\pi}^{\pi}m(w,z)cos(kw)dw,\ k\geq 1.$

And set $f_{2}(x,0)=c_{2},$ $\int_{-d}^{0}f_{2}(x,r)dr=c_{3}.$
Obviously, $m$ is a solution to (3.9) iff for $M$ the Sturm-Liouville problem
is solved by each $m_{k}=0$.

$\begin{cases}(a^{3}M_{z})_{z}=k^{2}aM,\quad z\in(0,z_{0}),\\\
\xi(\frac{Q}{2}-\int_{-d}^{0}f_{2}(x,r)dr)M_{z}=f_{2}(x,0)M,\quad\rm{on}\quad
z=z_{0},\\\ M=0,\quad on\quad z=0.\end{cases}$ (3.10)

As a special case, if $k=0$, the solution $m_{0}(z)$ of following equations is
zero.

$\begin{cases}(a^{3}{m_{0}}_{z})_{z}=k^{2}am_{0},\quad z\in(0,z_{0}),\\\
\xi(\frac{Q}{2}-\int_{-d}^{0}f_{2}(x,r)dr){m_{0}}_{z}=f_{2}(x,0)m_{0},\quad\rm{on}\quad
z=z_{0},\\\ m_{0}=0,\quad\rm{on}\quad z=0.\end{cases}$ (3.11)

The case $k=1$ is considered to find the $2\pi$-periodic solutions of (3.10).
According to the first equation of (3.10), we obtain

$-k^{2}\int_{0}^{z_{0}}aM^{2}dz=\int_{0}^{z_{0}}a^{3}M_{z}^{2}dz-2c_{1}a(z_{0})M^{2}(z_{0}).$

Denote

$\mu(\xi)=\inf\limits_{\zeta\in H^{1}(0,z_{0}),\zeta(0)=0,\zeta\neq
0}\mathbb{G}(\zeta,\xi)=\inf\limits_{\zeta\in
H^{1}(0,z_{0}),\zeta(0)=0,\zeta\neq
0}\frac{-\frac{c_{3}}{\xi(\frac{Q}{2}-c_{2})}a^{3}(z_{0})\zeta^{2}(z_{0})+\int_{0}^{z_{0}}a^{3}\zeta_{z}^{2}dz}{\int_{0}^{z_{0}}a\zeta^{2}dz},$
(3.12)

where $H^{1}$ is a Hilbert space. Obviously, a minimizer $M\in
C^{3,\alpha}[0,z_{0}]$ of (3.10) is a classical solution of following
equation.

$\begin{cases}(a^{3}M_{z})_{z}=-\mu(\xi)aM,\quad z\in(0,z_{0}),\\\
\xi\left(\frac{Q}{2}-\int_{-d}^{0}f_{2}(x,r)dr\right)M_{z}=f_{2}(x,0)M,\quad
z=z_{0},\\\ M=0,\quad z=0.\end{cases}$ (3.13)

We ensure that the solutions to (3.10) exist, if $\mu(\xi)=-1$ holds for some
$\xi>0.$

######  Lemma 3.1.

(i) The existence of $\mu(\xi)\in\mathbb{R}$ always holds for each
$\xi\in[\epsilon_{0},+\infty)$ with $0<\epsilon_{0}<1$ small.
(ii) There is an even function $M\in C^{3,\alpha}[0,z_{0}]$ with $2\pi$
periodic in $w$, such that the infimum in (3.12) is ensured. Furthermore, if
$\xi>2(c_{2}+c_{3})$, then the infimum $\mu(\xi)>-1$.

###### Proof.

Bear in mind that $a\leq\xi+\sqrt{2c_{1}z_{0}}$, we obtain

$\begin{split}\int_{0}^{z_{0}}a^{3}\zeta_{z}^{2}dz+\frac{4c_{3}^{2}(\xi+\sqrt{2c_{1}z_{0}})^{2}}{\xi^{2}(\frac{Q}{2}-c_{2})^{2}}\int_{0}^{z_{0}}a\zeta^{2}dz\geq&\frac{4c_{3}}{\xi(\frac{Q}{2}-c_{2})}\int_{0}^{z_{0}}a^{3}\zeta\zeta_{z}dz\\\
\geq&\frac{2c_{3}}{\xi(\frac{Q}{2}-c_{2})}a^{3}(z_{0})\zeta^{2}(z_{0}).\end{split}$
(3.14)

(3.14) implies
$\mu(\xi)\geq-\frac{4c_{3}^{2}(\xi+\sqrt{2c_{1}z_{0}})^{2}}{\xi^{2}(\frac{Q}{2}-c_{2})^{2}}.$
Now we prove that there is an even function $M\in C^{3,\alpha}[0,z_{0}]$ with
$2\pi$ periodic in $w$, such that the infimum in (3.12) is ensured. Without
loss of generality, we take a minimizing sequence $\zeta_{n}\in
H^{1}(0,z_{0})$ with $\zeta_{n}(0)=0$ satisfying
$G(\zeta_{n},\xi)\to\mu(\xi)$. And we normalize $\\{\zeta_{n}\\}$ such that

$\int_{0}^{z_{0}}a\zeta^{2}dz=1.$

Then

$\begin{split}G(\zeta_{n},\xi)&\geq-\frac{1}{2}\int_{0}^{z_{0}}a^{3}({\zeta_{n}}_{z})^{2}dz-\frac{2c_{3}^{2}(\xi+\sqrt{2c_{1}z_{0}})^{2}}{\xi^{2}(\frac{Q}{2}-c_{2})^{2}}\int_{0}^{z_{0}}a\zeta^{2}dz+\int_{0}^{z_{0}}a^{3}({\zeta_{n}}_{z})^{2}dz\\\
&=\frac{1}{2}\int_{0}^{z_{0}}a^{3}({\zeta_{n}}_{z})^{2}dz-\frac{2c_{3}^{2}(\xi+\sqrt{2c_{1}z_{0}})^{2}}{\xi^{2}(\frac{Q}{2}-c_{2})^{2}}\int_{0}^{z_{0}}a\zeta^{2}dz\\\
&\geq\frac{\xi^{3}}{2}\int_{0}^{z_{0}}(\zeta_{n})_{z}^{2}dz-\frac{2c_{3}^{2}(\xi+\sqrt{2c_{1}z_{0}})^{2}}{\xi^{2}(\frac{Q}{2}-c_{2})^{2}}\end{split}$

Since $G(\zeta_{n},\xi)\to\mu(\xi)$, $n\to\infty$,
$\\{\int_{0}^{z_{0}}({\zeta_{n}}_{z})^{2}dz\\}_{n\geq 1}$ is bounded. From

$\frac{1}{\xi}=\frac{1}{\xi}\int_{0}^{z_{0}}a\zeta_{n}^{2}dz\geq\int_{0}^{z_{0}}\zeta_{n}^{2}dz,$

$\\{\zeta_{n}\\}_{n\geq 1}$ is bounded in $H^{1}(0,z_{0})$, consequently,
there is a subsequence $\\{\zeta_{n_{k}}\\}$, which is weakly convergent to a
limit $M\in H^{1}(0,z_{0})$. It follows from $\frac{\partial}{\partial
z}\zeta_{n_{k}}\rightharpoonup M_{z}$ in $L^{2}(0,z_{0})$ and
$\zeta_{n_{k}}(0)=0$ that

$\zeta_{n_{k}}(z)\int_{0}^{z}\frac{\partial}{\partial
z}\zeta_{n_{k}}(r)dr\to\int_{0}^{z}M_{z}(r)dr=M(z),\quad z\in[0,z_{0}].$
(3.15)

Hence

$\int_{0}^{z_{0}}aM^{2}dz=\lim_{n\to\infty}\int_{0}^{z_{0}}a\zeta_{n}^{2}dz=1.$

The following property holds, because $\\{\zeta_{n_{k}}\\}$ is minimizing.

$\int_{0}^{z_{0}}a^{3}M_{z}^{2}dz\leq\liminf_{n_{k}\to\infty}\int_{0}^{z_{0}}a^{3}(\frac{\partial}{\partial
z}\zeta_{n_{k}})^{2}dz.$ (3.16)

In fact,

$\int_{0}^{z_{0}}a^{3}(\frac{\partial}{\partial
z}\zeta_{n_{k}})^{2}dz-\int_{0}^{z_{0}}a^{3}M_{z}^{2}dz=\int_{0}^{z_{0}}a^{3}(\frac{\partial}{\partial
z}\zeta_{n_{k}}-M_{z})^{2}dz+2\int_{0}^{z_{0}}a^{3}(\frac{\partial}{\partial
z}\zeta_{n_{k}})M_{z}dz-2\int_{0}^{z_{0}}a^{3}M_{z}^{2}dz$

and $a^{3}M_{z}\in L^{2}(0,z_{0})$ together with $\frac{\partial}{\partial
z}\zeta_{n_{k}}\rightharpoonup M_{z}$ weakly in $L^{2}(0,z_{0})$, then

$2\int_{0}^{z_{0}}a^{3}(\frac{\partial}{\partial
z}\zeta_{n_{k}})M_{z}dz-2\int_{0}^{z_{0}}a^{3}M_{z}^{2}dz\to 0,\quad
n_{k}\to\infty.$

According to (3.15) and (3.16), we obtain

$-\frac{c_{3}}{\xi(\frac{Q}{2}-c_{2})}a^{3}(z_{0})M^{2}(z_{0})+\int_{0}^{z_{0}}a^{3}M_{z}^{2}dz\leq\liminf_{n_{k}\to\infty}\left[-\frac{c_{3}}{\xi(\frac{Q}{2}-c_{2})}a^{3}(z_{0})\zeta_{n_{k}}^{2}(z_{0})+\int_{0}^{z_{0}}a^{3}(\frac{\partial}{\partial
z}{\zeta_{n_{k}}})^{2}dz\right]$

the infimum is a minimum when $M\in H^{1}(0,z_{0})$. By analogous arguments in
[12], $M$ is in $C^{3,\alpha}[0,z_{0}]$, because $\\{\zeta_{n_{k}}\\}$ is
minimizing for $G(\cdot,\xi)$. Keep in mind that $Q>2(c_{2}+c_{3})$, we obtain

$\begin{split}\int_{0}^{z_{0}}(a^{3}M_{z}^{2}dz+aM^{2})dz&\geq\xi\int_{0}^{z_{0}}(\xi^{2}M_{z}^{2}+M^{2})dz\\\
&\geq 2\xi M_{z}Mdz\\\ &\geq\xi^{2}M^{2}(z_{0})\\\
&\geq\frac{c_{3}M^{2}}{\xi(\frac{Q}{2}-c_{2})}a^{3}(z_{0})\zeta^{2}(z_{0}),\end{split}$

which implies $\mu(\xi)>-1$.

∎

######  Lemma 3.2.

On any interval such that $\mu(\xi)\leq-1$, if $\xi>2(c_{2}+c_{3})$, then
$\mu(\xi)$ is increasing in $\xi$. If $\mu(\xi)=-1$ exists a solution, then
the solution is unique, denoted by $\xi^{\star}$. Moreover, $\mu(\xi)>-1$, for
$\xi>\xi^{\star}$.

###### Proof.

The eigenfunction to the eigenvalue problem with eigenvalue $\mu(\xi)$ is
denoted by $M(\xi,z)$. Normalize $M(\xi,1)=1$, then

$M(\xi,z)=\zeta(\xi,z,\mu(\xi))$

where $\zeta$ is the unique solution of

$(a^{3}\zeta_{z})_{z}=-\mu a\zeta,\quad z\in(0,z_{0}),$ (3.17)

with initial data

$\zeta(z_{0})=1,\ \zeta^{\prime}(z_{0})=\frac{c_{3}}{\xi(\frac{Q}{2}-c_{2})}.$

The variational method ensures the existence of $\mu(\xi)$ and the follow
equation

$\frac{c_{3}}{\xi(\frac{Q}{2}-c_{2})}\xi^{3}-\int_{0}^{z_{0}}a^{3}\zeta_{z}^{2}dz+\mu\int_{0}^{z_{0}}a\zeta^{2}dz=0.$
(3.18)

By differentiating, (3.17) is transformed into

$-3c_{1}a\zeta_{z\mu}+a^{3}\zeta_{zz\mu}=-a\zeta-a\mu\zeta_{\mu}.$ (3.19)

Multiplying (3.19) by $\zeta=M$ and using integrating, we obtain

$\mu\int_{0}^{z_{0}}aMM_{\mu}dz-\int_{0}^{z_{0}}a^{3}M_{z}M_{z\mu}dz=-\int_{0}^{z_{0}}aM^{2}dz.$
(3.20)

Consequently, the derivative with respect to $\mu$ of the left hand side of
(3.17) is written as

$\int_{0}^{z_{0}}aM^{2}dz+2\mu\int_{0}^{z_{0}}aMM_{\mu}dz-2\int_{0}^{z_{0}}a^{3}M_{z}M_{z\mu}dz=-\int_{0}^{z_{0}}aM^{2}dz<0.$

Moreover, for $\xi\in(0,+\infty]$, the real-analytic dependence of $\mu(\xi)$
is deduced by implicit function theorem.
We denote $\dot{a}=\frac{\partial a}{\partial\xi}$, then

$\dot{a}=\frac{\xi}{a},\ \dot{a}_{z}=-\frac{\xi a_{z}}{a^{2}}.$

It follows from (3.13) that

$\begin{cases}3\xi(aM_{z})_{z}+(a^{3}\dot{M}_{z})_{z}=-\dot{\mu}(\xi)aM-\mu(\xi)\frac{\xi}{a}M-\mu(\xi)a\dot{M},\quad
0<z<z_{0},\\\ (\xi-
c_{2})M_{z}+\xi(\frac{\xi}{2}-c_{2})\dot{M}_{z}=c_{3}\dot{M},\quad z=z_{0},\\\
\dot{M}=0,\quad z=0.\end{cases}$ (3.21)

From the first equation of (3.21) and the first equation of (3.13), we obtain

$\int_{0}^{z_{0}}3\xi(aM_{z})_{z}Mdz+\int_{0}^{z_{0}}(a^{3}\dot{M}_{z})_{z}Mdz=-\int_{0}^{z_{0}}\dot{\mu}(\xi)aM^{2}dz-\int_{0}^{z_{0}}\mu(\xi)\frac{\xi}{a}M^{2}dz-\int_{0}^{z_{0}}\mu(\xi)a\dot{M}Mdz$

and

$\int_{0}^{z_{0}}\dot{M}(a^{3}M_{z})_{z}=-\int_{0}^{z_{0}}\mu(\xi)aM\dot{M}dz.$

Subtracting the two equations above, then

$\dot{\mu}\int_{0}^{z_{0}}aM^{2}dz=-\xi\mu\int_{0}^{z_{0}}a^{-1}M^{2}dz+3\xi\int_{0}^{z_{0}}aM_{z}^{2}dz-\frac{3c_{3}\xi}{\frac{\xi}{2}-c_{2}}+\frac{c_{3}a^{3}(\xi-
c_{2})}{\xi^{2}(\frac{\xi}{2}-c_{2})^{2}}.$

If $\mu\leq-1$, we use the Minkowski inequality and obtain that

$\begin{split}\dot{\mu}\int_{0}^{z_{0}}aM^{2}dz&\geq
2\sqrt{3}\xi|\mu|^{\frac{1}{2}}\left(\int_{0}^{z_{0}}a^{-1}M^{2}dz\right)^{\frac{1}{2}}\left(\int_{0}^{z_{0}}aM_{z}^{2}dz\right)^{\frac{1}{2}}-\frac{3c_{3}\xi}{\frac{\xi}{2}-c_{2}}+\frac{c_{3}a^{3}(\xi-
c_{2})}{\xi^{2}(\frac{\xi}{2}-c_{2})^{2}}\\\ &\geq
2\sqrt{3}\xi|\mu|^{\frac{1}{2}}\int_{0}^{z_{0}}MM_{z}dz-\frac{3c_{3}\xi}{\frac{\xi}{2}-c_{2}}+\frac{c_{3}a^{3}(\xi-
c_{2})}{\xi^{2}(\frac{\xi}{2}-c_{2})^{2}}\\\
&\geq\frac{\xi[\frac{\sqrt{3}}{4}\xi^{2}-(\sqrt{3}c_{2}+\frac{1}{2}c_{3})\xi+\sqrt{3}c_{2}^{2}+2c_{2}c_{3}]}{(\frac{\xi}{2}-c_{2})^{2}}.\end{split}$

Obviously, assume $\Delta=\frac{1}{4}c_{3}^{2}-\sqrt{3}c_{2}c_{3}<0.$ Then
$\frac{\sqrt{3}}{4}\xi^{2}-(\sqrt{3}c_{2}+\frac{1}{2}c_{3})\xi+\sqrt{3}c_{2}^{2}+2c_{2}c_{3}>0$
holds. Or else, if $\Delta=\frac{1}{4}c_{3}^{2}-\sqrt{3}c_{2}c_{3}\geq 0$, it
is easy to see that
$\frac{\sqrt{3}}{4}\xi^{2}-(\sqrt{3}c_{2}+\frac{1}{2}c_{3})\xi+\sqrt{3}c_{2}^{2}+2c_{2}c_{3}=0$
has two roots, namely,
$\xi_{1}=\frac{1}{\sqrt{3}}(2\sqrt{3}c_{2}+c_{3}-\sqrt{c_{3}^{2}-4\sqrt{3}c_{2}c_{3}})$,
$\xi_{2}=\frac{1}{\sqrt{3}}(2\sqrt{3}c_{2}+c_{3}+\sqrt{c_{3}^{2}-4\sqrt{3}c_{2}c_{3}})$.
Clearly, $\xi>2(c_{2}+c_{3})\geq\xi_{2}$, which implies
$\frac{\sqrt{3}}{4}\xi^{2}-(\sqrt{3}c_{2}+\frac{1}{2}c_{3})\xi+\sqrt{3}c_{2}^{2}+2c_{2}c_{3}>0$.
Therefore, $\dot{\mu}>0$ follows immediately. ∎

######  Lemma 3.3.

Assuming that the depth of vascular flow satisfies the following conditions

$\begin{split}\frac{\sqrt{2}}{24}c_{1}^{\frac{5}{2}}\left[\frac{2(c_{2}+c_{3})^{2}}{c_{1}}+\frac{d(dc_{1}+4c_{2}+4c_{3})}{2}\right]^{\frac{3}{2}}+\frac{\sqrt{2}}{20}c_{1}^{\frac{3}{2}}\left[\frac{2(c_{2}+c_{3})^{2}}{c_{1}}+\frac{d(dc_{1}+4c_{2}+4c_{3})}{2}\right]^{\frac{5}{2}}<(c_{2}+c_{3})^{4}.\end{split}$
(3.22)

Then $\mu(\xi)=-1$ exists a unique solution $\widetilde{\xi}$. Moreover,
$\widetilde{\xi}>2(c_{2}+c_{3}).$

###### Proof.

Denote

$\widetilde{\xi}=2(c_{2}+c_{3}),\quad\widetilde{z}=\frac{\widetilde{\xi}}{2c_{1}}+z_{0}.$
(3.23)

Then,

$a(z,\widetilde{\xi})=\sqrt{2c_{1}}\sqrt{\widetilde{z}-z},$

where $\widetilde{z}=\frac{\xi^{2}}{2c_{1}}+z_{0}.$
For $n>0$, set

$\Phi_{n}(z)=\begin{cases}0,\quad\quad\quad\quad 0\leq z\leq z_{n},\\\
(\widetilde{z}-z)^{\frac{1}{2}},\quad z_{n}\leq z\leq z_{0}.\end{cases}$
(3.24)

Obviously, $\zeta_{n}\in H^{1}(0,z_{0})$ and $\zeta_{n}(0)=0$, $\zeta_{n}\neq
0$.
By direct calculation, we obtain

$\begin{split}&\int_{0}^{z_{0}}a^{3}(z,\widetilde{\xi})(\frac{\partial\zeta_{n}(z)}{\partial
z})^{2}dz+\int_{0}^{z_{0}}a(z,\widetilde{\xi})\zeta_{n}^{2}(z)dz\\\
=&\int_{z_{n}}^{z_{0}}a^{3}(z,\widetilde{\xi})(\frac{\partial\zeta_{n}(z)}{\partial
z})^{2}dz+\int_{z_{n}}^{z_{0}}a(z,\widetilde{\xi})\zeta_{n}^{2}(z)dz\\\
=&\frac{1}{6}(2c)^{\frac{3}{2}}[(\frac{\widetilde{\xi}^{2}}{2c_{1}}+z_{0}-\frac{z_{0}}{n})^{\frac{3}{2}}-(\frac{\widetilde{\xi}^{2}}{2c_{1}})^{\frac{3}{2}}]+\frac{2}{5}\sqrt{2c_{1}}[(\frac{\widetilde{\xi}^{2}}{2c_{1}}+z_{0}-\frac{z_{0}}{n})^{\frac{5}{2}}-(\frac{\widetilde{\xi}^{2}}{2c_{1}})^{\frac{5}{2}}].\end{split}$

On the other hand,

$\begin{split}\frac{c_{3}}{\widetilde{\xi}(\frac{Q}{2}-c_{2})}a^{3}(z_{0})\zeta_{n}^{2}(z_{0})&=\frac{c_{3}}{\widetilde{\xi}(\frac{Q}{2}-c_{2})}(\sqrt{2c_{1}}\sqrt{\widetilde{z}-z_{0}})^{3}(\widetilde{z}-z_{0})\\\
&=\frac{c_{3}}{\widetilde{\xi}(\frac{\widetilde{\xi}}{2}-c_{2})}(2c_{1})^{\frac{3}{2}}(\frac{\widetilde{\xi}^{2}}{2c_{1}})^{\frac{5}{2}}.\end{split}$

Achieving our goal is equivalent to proof

$\int_{0}^{z_{0}}a^{3}(z,\widetilde{\xi})(\frac{\partial\zeta_{n}(z)}{\partial
z})^{2}dz+\int_{0}^{z_{0}}a(z,\widetilde{\xi})\zeta_{n}^{2}(z)dz<\frac{c_{3}}{\widetilde{\xi}(\frac{Q}{2}-c_{2})}a^{3}(z_{0})\zeta_{n}^{2}(z_{0}).$
(3.25)

The above inequality holds, if

$\frac{1}{6}(2c)^{\frac{3}{2}}(\frac{\widetilde{\xi}^{2}}{2c_{1}}+z_{0}-\frac{z_{0}}{n})^{\frac{3}{2}}+\frac{2}{5}\sqrt{2c_{1}}(\frac{\widetilde{\xi}^{2}}{2c_{1}}+z_{0}-\frac{z_{0}}{n})^{\frac{5}{2}}<\frac{c_{3}}{\widetilde{\xi}(\frac{\widetilde{\xi}}{2}-c_{2})}(2c_{1})^{\frac{3}{2}}(\frac{\widetilde{\xi}^{2}}{2c_{1}})^{\frac{5}{2}}.$

According to conditions (3.5) and (3.23), the above inequality can be proven
by (3.22). This proves the existence of $\zeta_{n}\in H^{1}(0,z_{0})$
satisfying (3.25). And for the specific $\widetilde{\xi}$ in (3.23),
$\mu(\xi)<-1$ holds. Note that $\mu(\xi)>-1$, if $\xi>2(c_{2}+c_{3})$. By the
continuity of function $\xi\to\mu(\xi)$, $\mu(\xi)=-1$ exists one solution
$\xi^{\star}$. ∎

###  3.3 Local bifurcation

In this section, we will apply Crandall-Rabinowitz bifurcation theorem ([15],
pp325, Theorem 1.7) to study the local bifurcation of the non-trivial
solutions of the blood flow model. In order to apply Crandall-Rabinowitz
bifurcation theorem, let $B_{T}=\\{(w,z):w\in[-\pi,\pi],z=z_{0}\\}\\}$ be the
top boundary and $B_{B}=\\{(w,z):w\in[-\pi,\pi],z=0\\}$ be the bottom
boundary. And we define the Banach spaces

$\mathbb{X}=\\{\gamma\in
C_{per}^{3,\alpha}(\bar{R}_{z}):\int_{-\pi}^{\pi}\gamma(w,z_{0})=d,\ \gamma=0\
on\ B_{B}\\},$ $\mathbb{Y}=C_{per}^{1,\alpha}(\bar{R}_{z})\times
C^{2,\alpha}_{per}(B_{T}).$

Moreover, system (2.25) can be rewritten as

$G(\xi,\gamma)=0\quad with\quad\gamma\in\mathbb{X},$

where $G:(\epsilon_{0},z_{0})\times\mathbb{X}\to\mathbb{Y}$ is given by

$G_{1}(\xi,\gamma)=(1+\gamma_{w}^{2})(H_{zz}+\gamma_{zz})-2\gamma_{w}\gamma_{zw}(H_{z}+\gamma_{z})+(H_{z}+\gamma_{z})^{2}\gamma_{ww}-c_{1}(H_{z}+\gamma_{z})^{3}$

and

$G_{2}(\xi,\gamma)=1+\gamma_{w}^{2}-\left[Q-2F(w,H+\gamma-d)\right](H_{z}+\gamma_{z})\bigg{|}_{B_{B}}.$

The linearized operator
$G_{\gamma}(\xi,0)=(G_{1\gamma},G_{2\gamma})\Bigg{|}_{\gamma=0}$ is given by

$\begin{cases}G_{1\gamma}(\xi,0)=\frac{\partial^{2}}{\partial
z^{2}}+H^{2}_{z}\frac{\partial^{2}}{\partial
w^{2}}-3c_{1}H^{2}_{z}\frac{\partial}{\partial z},\\\
G_{2\gamma}(\xi,0)=\frac{2f_{2}(w,0)}{\xi}-[\xi-2\int_{-d}^{0}f_{2}(x,r)dr]\frac{\partial}{\partial
z}\bigg{|}_{B_{T}}.\end{cases}$ (3.26)

According to Crandall-Rabinowitz bifurcation theorem ([15], pp325, Theorem
1.7), in what follows, we divide the proof of the local bifurcation into two
steps. Firstly, we prove that the null space
$ker\\{G_{\gamma}(\xi^{\star},0)\\}$ is one-dimensional. Secondly, we prove
the transversality condition $[G_{\xi\gamma}(\xi^{*},0)](1,\zeta^{*})\notin
R(G_{\gamma}(\xi^{*},0))$ holds.

######  Lemma 3.4.

The null space $ker\\{G_{\gamma}(\xi^{\star},0)\\}$ is one-dimensional.

###### Proof.

Based on the limitations on depth in (3.22), there exists a unique solution
$\xi^{\star}$ of the equation $\mu(\xi)=-1$ with $\xi>\frac{4}{\sqrt{3}d+4}.$
Moreover, (3.13) has unique normalized ($M(z_{0})=1$) eigenfunction $M\in
C^{3,\alpha}[0,z_{0}]$ such that $\mu(\xi)=-1$. Then there exists at least one
element $\gamma(w,z)=M(z)\cos w$ in the space
$ker\\{G_{\gamma}(\xi^{\star},0)\\}$.
By assuming that $m\in X$ is in the null space
$ker{G_{\gamma}(\xi^{\star},0)}$, we obtain that $m_{k}$ solve the equation
(3.10). Therefore, $m_{1}(z)$ is a constant multiple of $M(z)$. Moreover, if
$m_{k}\neq 0$ for some $k>2$, then

$\frac{-\frac{c_{3}}{\xi(\frac{Q}{2}-c_{2})}+\int_{0}^{z_{0}}a^{3}(\frac{\partial}{\partial
z}m_{k})^{2}dz}{\int_{0}^{z_{0}}d^{2}am_{k}^{2}dz}=-k^{2}<-1.$

This is a contradiction, because a minimum is attained corresponding to
$\mu(\xi^{\star})=-1$. And from the differential equation and from the
boundary condition at $z=0$ in (3.11) with $k=0$ we get that

$m_{0}(z)=A_{0}\int_{0}^{z}a^{-3}(r,\xi^{\star})dr,\quad z\in[0,z_{0}]$

holds for some $A_{0}\in\mathbb{R}$. The condition
$\int_{\pi}^{\pi}m_{0}(z_{0})dw=2\pi m_{0}(z_{0})=0$ and
$A_{0}\int_{0}^{z_{0}}a^{-3}(z,\xi^{\star})dz=0.$ However, it is not difficult
to show that $\int_{0}^{z_{0}}a^{-3}(z,\xi^{\star})dz\neq 0$ by direct
computation. Hence $A_{0}=0$ and $m_{0}=0$ ∎

Based on the method in [12, 14], for $\xi\in(0,z_{0})$ the derivative
$G_{\gamma}:\mathbb{X}\to\mathbb{Y}$ is a Fredholm operator of index zero. And
we claim that $G_{\gamma}(\xi^{\star},0)$ has a closed range of codimension
one.

######  Lemma 3.5.

The following transversality condition holds.

$[G_{\xi\gamma}(\xi^{*},0)](1,\zeta^{*})\notin R(F_{\gamma}(\xi^{*},0)).$

###### Proof.

Keep in mind that $a_{z}=-2(1-\xi^{*})a^{-1}$ where $a=a(\cdot,\xi^{*}).$
It follows from (3.26) that

$G_{\xi\gamma}(\xi^{*},0)=-\left(2\xi^{*}a^{-4}\frac{\partial^{2}}{\partial
w^{2}}+6\xi^{*}a^{-3}a_{z}\frac{\partial}{\partial
z},\left(\frac{2f_{2}(x,0)}{{\xi^{*}}^{2}}+\frac{\partial}{\partial
z}\right)\bigg{|}_{B_{T}}\right).$

Lemma is equivalent to proof

$\iint_{R_{z}}a^{3}\zeta^{*}2\xi^{*}(a^{-4}\frac{\partial^{2}\zeta^{*}}{\partial
w^{2}}+3a^{-3}a_{z}\frac{\partial\zeta^{*}}{\partial
z})dwdz+\int_{B_{T}}\\{2f_{2}(x,0){\zeta^{*}}^{2}+{\xi^{*}}^{2}\zeta^{*}\zeta^{*}_{z}\\}dw\neq
0.$ (3.27)

Note that $\zeta^{*}(w,z)=M(z)\cos(w)$. Since $M\neq 0$ solves (3.13) with
$\mu=-1$, and $a(z_{0})=\xi^{*}$, then we obtain

$\begin{cases}a^{3}\zeta_{zz}^{*}+3a^{2}a_{z}\zeta_{z}^{*}-a\zeta^{*}=0,\\\
\xi^{*}(\frac{Q}{2}-c_{3})\zeta^{*}_{z}=c_{2}\zeta^{*},\\\
\zeta^{*}=0.\end{cases}$ (3.28)

Hence, we obtain

$\begin{split}\iint_{R_{z}}a_{z}\zeta^{*}\zeta_{z}^{*}dzdw&=\int_{B_{T}}a\zeta^{*}\zeta_{z}^{*}dw-\iint_{R_{z}}a\zeta^{*}\zeta_{zz}^{*}dwdz-\iint_{R_{z}}a(\zeta^{*}_{z})^{2}dwdz\\\
&=\int_{B_{T}}a\zeta^{*}\zeta^{*}_{z}dw+\iint_{R_{z}}\zeta^{*}(3a_{z}\zeta^{*}_{z}-a^{-1}\zeta^{*})dwdz-\iint_{R_{z}}a(\zeta^{*}_{z})^{2}dwdz.\end{split}$

Furthermore,

$\iint_{R_{z}}a_{z}\zeta^{*}\zeta^{*}_{z}dwdz=-\frac{1}{2}\int_{B_{T}}a\zeta\zeta_{z}dw+\frac{1}{2}\iint_{R_{z}}\left({a^{-1}(\zeta^{*})^{2}+a(\zeta_{z}^{*})^{2}}\right)dwdz.$

Keep in mind $\zeta_{ww}=-\zeta$ and $a(z_{0})=\xi$. And form (3.27), we need
to show that

$-2\xi^{*}\iint_{R_{z}}a^{-1}(\zeta^{*})^{2}dwdz+3\xi^{*}\iint_{R_{z}}\left({a^{-1}(\zeta^{*})^{2}+a(\zeta_{z}^{*})^{2}}\right)dwdz+\int_{B_{T}}2c_{2}\zeta^{2}dw-2\xi\int_{B_{T}}a\zeta^{*}\zeta^{*}_{z}dw\neq
0,$

which is equal to

$\iint_{R_{z}}\xi^{*}a^{-1}{\zeta^{*}}^{2}dwdz+3\xi^{*}\iint_{R_{z}}a(\zeta_{z}^{2})dwdz+\int_{B_{T}}2c_{2}[1-\xi^{*}(\frac{Q}{2}-c_{3})^{-1}](\zeta^{*})^{2}dw\neq
0.$

With the average inequality and Cauchy-Schwartz inequality, we obtain

$\begin{split}&\iint_{R_{z}}\xi^{*}a^{-1}{\zeta^{*}}^{2}dwdz+3\xi^{*}\iint_{R_{z}}a(\zeta_{z}^{2})dwdz+\int_{B_{T}}2c_{2}[1-\xi^{*}(\frac{Q}{2}-c_{3})^{-1}](\zeta^{*})^{2}dw\\\
\geq&2\sqrt{3}\xi^{*}\left\\{\iint_{R_{z}}a^{-1}{\zeta^{*}}^{2}dwdz\right\\}^{\frac{1}{2}}\left\\{\iint_{R_{z}}a(\zeta_{z}^{2})dwdz\right\\}^{\frac{1}{2}}+\int_{B_{T}}2c_{2}[1-\xi^{*}(\frac{Q}{2}-c_{3})^{-1}](\zeta^{*})^{2}dw\\\
\geq&2\sqrt{3}\xi^{*}\iint_{R_{z}}a^{-\frac{1}{2}}{\zeta^{*}}a^{\frac{1}{2}}\zeta^{*}_{z}dwdz+\int_{B_{T}}2c_{2}[1-\xi^{*}(\frac{Q}{2}-c_{3})^{-1}](\zeta^{*})^{2}dw\\\
=&\sqrt{3}\xi^{*}\int_{B_{T}}({\zeta^{*}})^{2}dw+\int_{B_{T}}2c_{2}[1-\xi^{*}(\frac{Q}{2}-c_{3})^{-1}](\zeta^{*})^{2}dw\\\
=&\left\\{2c_{2}+[\sqrt{3}-c_{2}(\frac{\xi}{2}-c_{3})^{-1}]\xi\right\\}\int_{B_{T}}({\zeta^{*}})^{2}dw.\end{split}$

The last item is larger than 0, guaranteed by the condition
$\xi>2(c_{2}+c_{3}).$

∎

Through the above analysis, the main result of this paper is proven. We
conclude it as follows.

######  Theorem 3.6.

If the depth of blood flow satisfies the condition (3.22). Then there exists a
$C^{1}$-curve $C_{loc}$ of small-amplitude solution $h\in
C^{3,\alpha}_{per}(\widetilde{R}_{z})$ and $(u,v,\Omega)\in
C^{2,\alpha}_{per}(\widetilde{R}_{z})\times
C^{2,\alpha}_{per}(\widetilde{R}_{z})\times C^{3,\alpha}_{per}(R)$. Moreover,
the solution curve $C_{loc}$ contains precisely one function that is
independent of $w$.

## Conflict of Interest

The authors declare that they have no conflict of interest.

## Data Availability Statement

No data was used for the research in this article.

## Contributions

We declare that all the authors have same contributions to this paper.

## References

  * [1] A. Ali, M. Hussain, M. Anwar, M. Inc, Mathematical modeling and parametric investigation of blood flow through a stenosis artery, Appl. Math. Mech. Engl. Ed., 42(11) (2021) 1675-1684.
  * [2] B. Basu, A flow force reformulation for steady irrotational water waves, J. Differential Equations, 268 (2020) 7417-7452.
  * [3] B. Basu, C. Martin, Capillary-gravity water waves: modified flow force formulation, J. Differential Equations, 269 (2020) 11231–11251.
  * [4] Y. Bi, Z. Zhang, Q. Liud, T. Liu, Research on nonlinear waves of blood flow in arterial vessels, Commun. Nonlinear Sci. Numer. Simul., 102 (2021) 105918.
  * [5] A. Constantin and J. Escher, Analyticity of periodic traveling free surface water waves with vorticity, Ann. Math., (2) 173 (2011) 559-568.
  * [6] A. Constantin, W. Strauss, Exact steady periodic water waves with vorticity, Commun. Pure Appl. Math., 57 (2004) 481-527.
  * [7] A. Constantin, W. Strauss, E. V$\breve{a}$rv$\breve{a}$ruc$\breve{a}$, Global bifurcation of steady gravity water waves with critical layers, Acta Math., 217 (2016) 195-262.
  * [8] A. Constantin, E. V$\breve{a}$rv$\breve{a}$ruc$\breve{a}$, Steady periodic water waves with constant vorticity: regularity and local bifurcation, Arch. Ration. Mech. Anal., 199 (2011) 33–67.
  * [9] A. Constantin, Some three-dimensional nonlinear equatorial flows, J. Phys. Oceanogr., 43 (2013) 165-175.
  * [10] A. Constantin, M. Ehrnström, E. Wahlén, Symmetry of steady periodic gravity water waves with vorticity, Duke Math. J., 140 (2007) 591-603.
  * [11] A. Constantin, Mean velocities in a Stokes wave, Arch. Ration. Mech. Anal., 207 (2013) 907-917.
  * [12] J. Chu, X. Wang, L. Wang, Z. Zhang, A flow force reformulation of steady periodic fixed-depth irrotational equatorial flows, J. Differential Equations, 292 (2021) 220-246.
  * [13] J. Chu, Q. Ding, Y. Yang, Steady periodic waves and formal stability for fixed-depth rotational equatorial flows, J. Differ. Equ., 269 (2020) 4192-4214.
  * [14] J. Chu, J. Escher, Steady periodic equatorial water waves with vorticity, Discrete Contin. Dyn. Syst. 39 (2019) 4713–4729.
  * [15] M. Crandall, P. Rabinowitz, Bifurcation from simple eigenvalues, J. Funct. Anal., 8 (1971) 321-340.
  * [16] B. Dong, Z. Zhang, Global regularity of the 2D micropolar fluid flows with zero angular viscosity, J. Differential Equations, 249 (2010) 200-213.
  * [17] L. Euler, Principia pro motu sanguinis per arterias determinando. Opera posthuma mathematica et physica anno 1844 detecta, 2 (1775) 814-823.
  * [18] C. Foias, J. Saut, L. Num$\acute{e}$rique, Linearization and normal form of the Navier-Stokes equations with potential forces, Ann. Inst. Henri Poincar$\acute{e}$, 4 (1987) 1-47.
  * [19] A. Friedman, W. hao, B. Hu, A free boundary problem for steady small plaques in the artery and their stability, J. Differential Equations, 259 (2015), 1227-1255.
  * [20] R. Ghigo, J. Fullana, P. Lagrée, A 2D nonlinear multiring model for blood flow in large elastic arteries, J. Comput. Phys., 350 (2017) 136-165.
  * [21] D. Henry, Steady periodic waves bifurcating for fixed-depth rotational flows, Q. Appl. Math., 71 (2013) 455-487.
  * [22] D. Henry, G. Villari, Flow underlying coupled surface and internal waves, J. Differential Equations, 310 (2022) 404-442.
  * [23] L. Hoang, V. Martinez, Asymptotic expansion for solutions of the Navier-Stokes equations with non-potential body forces, J. Math. Anal. Appl., 462 (2018) 84-113.
  * [24] Q. Li, M. Fečkan, J. Wang, Velocity field and pressure analysis of equatorial flows beneath solitary water waves, Nonlinear Anal-Real, 69 (2023) 103754.
  * [25] L. Müller, C. Parés, E. Toro, Well-balanced high-order numerical schemes for one-dimensional blood flow in vessels with varying mechanical properties, J. Comput. Phys., 242 (2013) 53-85.
  * [26] B. Riemann, Gesammelte mathematische Werke und wissenschaftlicher Nachlass. Leipzig, B.G. Teubner (ed.), (1876) 145-164.
  * [27] J. Stettler, P. Niederer and M. Anliker, Theoretical analysis of arterial hemodynamics including the in- fluence of bifurcations. part i: Mathematical model and prediction of normal pulse patterns. Ann. Biomed. Engng., 9 (1981) 145-164.
  * [28] R. Skalak, The synthesis of a complete circulation, In: D. Bergel (ed.), Cardiovascular Fluid Dynamics (vol. 2). London: Academic Press, (1972) 341-376.
  * [29] V. Srivastava, R. Rastogi, R. Vishnoi. A two-layered suspension blood flow through an overlapping stenosis, Comput. Math. Appl., 60 (2010) 432-441.
  * [30] A. Quarteroni, A. Veneziani, C. Vergara, Geometric multiscale modeling of the cardiovascular system, between theory and practice, Comput. Methods Appl. Mech. Eng., 302 (2016) 193-252.
  * [31] A. Quarteroni and L. Formaggia. Modelling of Living Systems, chapter Mathematical Modelling and Numerical Simulation of the Cardiovascular System. Handbook of Numerical Analysis. Elsevier Science, Amsterdam, 2002.
  * [32] C. Tabi, T. Motsumi, C. Kamdema, A. Mohamadoue, Nonlinear excitations of blood flow in large vessels under thermal radiations and uniform magnetic field, Commun. Nonlinear Sci. Numer. Simul., 49 (2017) 1-8.
  * [33] S. Walsh, Stratified steady periodic water waves, SIAM J. Math. Anal., 41 (2009) 1054–1105.
  * [34] S. Walsh, Steady stratified periodic gravity waves with surface tension I: local bifurcation, Discrete Contin. Dyn. Syst., 34 (2014) 3241–3285.
  * [35] S. Walsh, Steady stratified periodic gravity waves with surface tension II: global bifurcation, Discrete Contin. Dyn. Syst., 34 (2014) 3287–3315.
  * [36] N. Xiao, J. Humphrey, C. Figueroa, Multi-scale computational model of three-dimensional hemodynamics within a deformable full-body arterial network, J. Comput. Phys., 244 (2013) 22-40.
  * [37] T. Young, Hydraulic investigations, subservient to an intended Croonian lecture on the motion of the blood. Phil. Trans. R. Soc. London, 98 (1808), 164-186.


& \leq \sum_{v \in P_\ell}\frac{\dtsize(\T_v)}{\dtsize(\T_{N(v)})}\nonumber\\
& = \sum_{v \in P_\ell} 1 + \frac{\dtsize(\T_v) - \dtsize(\T_{N(v)})}{\dtsize(\T_{N(v)})}\nonumber\\
& \leq \sum_{v \in P_\ell} 1 + \dtsize(\T_v) - \dtsize(\T_{N(v)})\nonumber\\
& = |P_\ell| + \dtsize(\T) - 1\nonumber\\
& \leq 2\dtsize(\T)\nonumber.
\end{align}

Thus, this setting of weights, along with the settings

\[
\alpha = 2\dtsize(\T), \qquad \beta = \log \dtsize(\T),
\]

gives an objective value of $O(\sqrt{\dtsize(\T) \log \dtsize(\T)})$ Program <ref> for $\T$.