#  Gilbarg-Serrin Equation and Lipschitz Regularity

Vladimir Maz’ya  Department of Mathematics, Linköping University, SE-581 83
Linköping, Sweden and RUDN, 6 Miklukho-Maklay St, Moscow, 117198, Russia
<EMAIL_ADDRESS>and  Robert McOwen  Department of Mathematics,
Northeastern University, Boston, MA 02115<EMAIL_ADDRESS>Dedicated
in memory of James Serrin who passed away ten years ago.

(Date: September 11, 2021)

###### Abstract.

We discuss conditions for Lipschitz and $C^{1}$ regularity of solutions for a
uniformly elliptic equation in divergence form. We focus on coefficients
having the form that was introduced by Gilbarg & Serrin. In particular, we
find cases where Lipschitz or $C^{1}$ regularity holds but the coefficients
are not Dini continuous, or do not even have Dini mean oscillation. The form
of the coefficients also enables us to obtain specific conditions and examples
for which there exists a weak solution that is not Lipschitz continuous.

###### Key words and phrases:

Lipschitz continuity, Dini condition, square-Dini condition, Dini mean
oscillation, weak solutions, asymptotic analysis

###### 1991 Mathematics Subject Classification:

35B40, 35B45, 35B65, 35J15

##  1\. Introduction and General Theory

The topic of this paper is Lipschitz and $C^{1}$ regularity for weak solutions
of a uniformly elliptic equation in divergence form:

(1) $\partial_{j}(a_{ij}(x)\,\partial_{i}u)=0,\quad\hbox{for $x\in\Omega$},$

where the coefficients $a_{ij}=a_{ji}$ are bounded, measurable functions in an
open set $\Omega\subset{\mathbb{R}}^{n}$, $n\geq 2$, which contains the origin
$x=0$. For a restricted class of coefficients $a_{ij}$, we seek conditions
that guarantee that any weak solution $u\in
H_{{\ell\mathrm{o}c}}^{1,2}(\Omega)$, i.e. $u$ and $\nabla u$ are locally in
$L^{2}$, must be Lipschitz continuous at $x=0$, i.e. $|u(x)-u(0)|\leq c\,|x|$
holds for $x\in B_{\rho}$ with $\rho$ sufficiently small, or Lipschitz
continuous in a neighborhood of $x=0$, i.e. $|u(x)-u(y)|\leq c\,|x-y|$ for
$x,y\in B_{\rho}$ with $\rho$ sufficiently small. (Here and throughout this
paper, $c$ denotes a positive constant and $B_{\rho}$ denotes the ball of
radius $\rho$ centered at $x=0$ such that $\overline{B_{\rho}}\subset\Omega$;
the values of $c$ and $\rho$ may change with each occurrence.) Recall that
Lipschitz continuous functions are differentiable almost everywhere (by
Rademacher’s theorem) and so our solution will satisfy $\nabla u\in
L^{\infty}_{{\ell\mathrm{o}c}}(\Omega)$. We also seek conditions implying that
any weak solution is differentiable at $x=0$, or is $C^{1}$ near $x=0$.

###  1.1. Background

Let us review a little of what is known about the regularity of weak solutions
of (1). Any weak solution of (1) is known to be Hölder continuous, i.e.
$|u(x)-u(y)|\leq c\,|x-y|^{\alpha}$ for some $\alpha\in(0,1)$, at least
locally in $\Omega$ (cf. De Giorgi [11], Nash [35], Moser [34], Landis [25]).
When the coefficients are continuous in $\Omega$, then it is well-known (cf.
Agmon, Douglis, Nirenberg [1]) that $\nabla u\in
L^{p}_{{\ell\mathrm{o}c}}(\Omega)$ for $1<p<\infty$; in fact, this is even
true when the coefficients are in VMO (cf. Chiarenza, Frasca, Longo [4], Di
Fazio [12]), or in BMO but sufficiently close (depending on $pp^{\prime}$) to
VMO (cf. Maz’ya, Mitrea, Shaposhnikova [31]).111Here, as usual, BMO stands for
bounded mean oscillation and VMO stands for vanishing mean oscillation. Note
that continuity of the coefficients is not sufficient to show $\nabla u\in
L^{\infty}_{{\ell\mathrm{o}c}}(\Omega)$ (cf. Jin, Maz’ya, Van Schaftingen
[20]). However, if the coefficients are Dini continuous in $\Omega$, then a
weak solution is known to be not only Lipschitz continuous but $C^{1}$ (cf.
Hartman, Wintner [18], Burch [3], Taylor [37]). Special attention has been
given to pointwise bounds on the gradient of the solutions for linear and
nonlinear elliptic equations and systems in [28], [5], [23], [6], [24], [7],
[8]. The Dini condition also provides regularity in parabolic problems: [26],
[19], [2].

The classical Dini condition on the coefficients involves the modulus of
continuity $\omega$:

(2) $|a_{ij}(x)-a_{ij}(y)|\leq\omega(|x-y|)\quad\hbox{for $x,y\in\Omega$,}$

where $\omega(r)$ is a continuous, nondecreasing function
$[0,\infty)\to[0,\infty)$ with $\omega(0)=0$. The coefficients are said to be
Dini continuous if $\omega$ satisfies

(3) $\int_{0}^{1/2}\frac{\omega(r)}{r}\,dr<\infty.$

When $\omega(r)=r^{\alpha}$ for some $\alpha\in(0,1)$ then this regularity
coincides with Hölder continuity. There are weaker moduli of continuity that
satisfy the Dini condition. If we let $\omega(r)=|\log r|^{-\gamma}$ for
$0<r<1/2$ and $\omega(0)=0$, then $\omega$ satisfies (3) if $\gamma>1$.
However, if $\omega(r)=|\log r|^{-\gamma}$ for $0<\gamma\leq 1$, then the
coefficients are not Dini continuous, and the usual regularity properties of
solutions may not hold.

More recently, generalizations of the Dini condition have been used. For
example, consider the “square-Dini condition” on the modulus of continuity
$\omega$:

(4) $\int_{0}^{1/2}\frac{\omega^{2}(r)}{r}\,dr<\infty.$

A function is square-Dini continuous if its modulus of continuity $\omega$
satisfies (4). Clearly the condition (4) is more general than (3); for
example, $\omega(r)=|\log r|^{-\gamma}$ satisfies (4) when $\gamma>1/2$. In
[28], the square-Dini condition was only required at a point (see Section 1.2)
and was used to reduce the regularity of a weak solution at that point to
stability properties of a dynamical system obtained from the coefficients.
This dynamical systems approach to studying the properties of solutions of
elliptic equations has been used in many other publications: [21], [22], [33],
[29], [30]. We shall have more to say about this dynamical system method
below, but for now let us mention another generalization of the Dini
condition.

Several authors have considered a “Dini mean oscillation” condition on the
coefficients. In the comprehensive paper [23], Kuusi & Mingione used the
notion of Dini mean oscillation in their unified linear/nonlinear potential
theory to obtain gradient bounds for solutions of quasilinear elliptic
equations. For linear equations the condition was used by Dong & Kim [15] to
show that weak solutions of an elliptic equation in divergence form are
$C^{1}$ and for nondivergence form they are $C^{2}$. Let us formulate the
condition on the coefficients as follows: a function $a$ that is integrable in
$\Omega$ has Dini mean oscillation if the following function $\omega_{a}$
satisfies the Dini condition (3):

(5) $\omega_{a}(r):=\sup_{B_{r}(x)\subset\Omega}\
{\diagup\hskip-11.9501pt\int}_{B_{r}(x)}|a(y)-\widetilde{a}_{x,r}|\,dy\quad\hbox{where}\
\widetilde{a}_{x,r}:={\diagup\hskip-11.9501pt\int}_{B_{r}(x)}a(y)\,dy.$

(Here $B_{r}(x)=\\{y:|x-y|<r\\}$ and the slashed integral denotes mean value.)
There are functions that are not Dini continuous but do have Dini mean
oscillation; for example, if $a(x)=(-\log|x|)^{-\gamma}$ in $\Omega=B_{1/2}$,
then $\omega_{a}(r)=c\,|\log r|^{-\gamma-1}$, so $a$ has Dini mean oscillation
for all $\gamma>0$. The paper [15] has spawned a series of applications to
various linear elliptic and parabolic equations: [13], [14], [16]. As useful
as the concept of Dini mean oscillation has proven to be, however, we will see
that there are equations of the form (1) with coefficients that do not have
Dini mean oscillation, yet Lipschitz and $C^{1}$ regularity can be obtained by
other methods.

###  1.2. Dynamical Systems Method

Now let us describe the method in [28] of using dynamical systems to study the
regularity at a given point of weak solutions of (1). We remark that the
conditions on the coefficients need only occur at the given point. These
conditions and the dynamical system are most conveniently formulated when the
point in question is $x=0$ and $a_{ij}(0)=\delta_{ij}$; this can always be
achieved by a change of coordinates at $0$ since we have assumed the
coefficients are symmetric. Under these assumptions, the condition on the
coefficients is

(6) $\sup_{|x|=r}|a_{ij}(x)-\delta_{ij}|\leq\omega(r)\quad\hbox{as}\ r\to 0,$

where $\omega$ satisfies the square-Dini condition (4); we describe this by
saying the coefficients are square-Dini continuous at $x=0$.

The basic idea in [28] is to write the solution as

(7) $u(x)=u_{0}(|x|)+\vec{v}(|x|)\cdot x+w(x),$

where

$u_{0}(r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}u(r\theta)\,d\theta\quad\hbox{and}\quad\vec{v}_{k}(r)=\frac{n}{r}{\diagup\hskip-11.9501pt\int}_{S^{n-1}}u(r\theta)\theta_{k}\,d\theta;$

here $r=|x|$, $\theta=x/|x|\in S^{n-1}$, and $d\theta$ is standard surface
measure on $S^{n-1}$. As shown in [28], if the coefficients are square-Dini
continuous at $x=0$, then the behaviors of both $\vec{v}(r)$ and
$r\,\vec{v}\,^{\prime}(r)$ are controlled by the asymptotic properties of the
solutions of the $(n\times n)$ dynamical system

(8) $\frac{d\phi}{dt}+R(e^{-t})\,\phi=0\quad\hbox{for}\ T<t<\infty,$

where $t=-\log r$, $T$ is sufficiently large, and $R$ is the ($n\times
n$)-matrix

(9)
$R(r):={\diagup\hskip-11.9501pt\int}_{S^{n-1}}\left(A(r\theta)-nA(r\theta)\theta\otimes\theta\right)\,d\theta,\quad\hbox{for
$0<r<1$,}$

with $A=(a_{ij})$ and $A\theta\otimes\theta$ the outer product of the vectors
$A\theta$ and $\theta$. Note that $\boldsymbol{|}R(r)\boldsymbol{|}\leq
c\,\omega(r)$, where $\boldsymbol{|}\cdot\boldsymbol{|}$ denotes the matrix
norm; also note that $R$ need not be symmetric. The asymptotic property that
we require for solutions to be Lipschitz is uniform stability: we say that (8)
is uniformly stable as $t\to\infty$ if for every $\varepsilon>0$ there exists
$\delta=\delta(\varepsilon)>0$ such that any solution of (8) satisfying
$|\phi(t_{1})|<\delta$ for some $t_{1}>T$ satisfies $|\phi(t)|<\varepsilon$
for all $t>t_{1}$. We show in [28] that the uniform stability of (8) implies
not only that $\vec{v}(r)$ and $r\,\vec{v}\,^{\prime}(r)$ are bounded as $r\to
0$, but with additional analysis that
$|u_{0}(r)-u_{0}(0)|,|w(x)|=O(r\,\omega(r))$ as $r\to 0$. This means that

(10) $u(x)=u_{0}(0)+\vec{v}(|x|)\cdot x+O(r\,\omega(r)),$

which confirms that $u$ is Lipschitz at $x=0$.

For differentiabilty, we also require that every solution of (8) be
asymptotically constant, i.e. $\phi(t)\to\phi_{\infty}$ as $t\to\infty$. In
that case, we have $\vec{v}(r)=\vec{v}(0)+o(1)$ as $r\to 0$, which shows that
$u$ is differentiable at $x=0$. Here is the complete result:

######  Theorem 1 (Theorem 1 in [28]).

Suppose the $a_{ij}$ are bounded, measurable functions in $\Omega$ that are
square-Dini continuous at $x=0$. Suppose also that the dynamical system (8) is
uniformly stable. Then every weak solution $u\in H^{1,2}(\Omega)$ of (1) is
Lipschitz continuous at $x=0$ and satisfies

(11)
$|u(x)-u(0)|\leq\frac{c\,|x|}{\rho}\left({\diagup\hskip-11.9501pt\int}_{|y|<\rho}|u(y)|^{2}\,dy\right)^{1/2}\quad\hbox{for}\
|x|<\rho/2,$

for $\rho$ sufficiently small. In addition, if every solution of (8) is
asymptotically constant, then $u$ is differentiable at $x=0$ and

(12) $\partial_{j}u(0)=\lim_{r\to
0}\frac{n}{r}\,{\diagup\hskip-11.9501pt\int}_{S^{n-1}}u(r\theta)\theta_{j}d\theta,\quad\hbox{where}\
\theta_{j}=x_{j}/r.$

We also gave in [28] several simpler sufficient conditions for the uniform
stability of (8), such as

(13) $\frac{R(r)}{r}\int_{0}^{r}\frac{R(\rho)}{\rho}\,d\rho\in
L^{1}(0,\varepsilon),$

or

(14)
$\displaystyle\int_{r_{1}}^{r_{2}}\frac{\mu({\mathcal{S}}(\rho))}{\rho}\,d\rho\leq
K<\infty\quad\hbox{for all}\ 0<r_{1}<r_{2}<\varepsilon,$ where
$\mu({\mathcal{S}})=$ largest eigenvalue of ${\mathcal{S}}:=-(R+R^{t})/2$.

The condition that all solutions of a dynamical system are asymptotically
constant is independent of its uniform stability (cf. Section 5 in [28]); when
(8) is a scalar equation, as we shall have in Section 2.1, it is elementary to
completely characterize this condition.

Also note that, if the $a_{ij}$ are radial functions, then $R=0$ and Theorem 1
implies Lipschitz regularity and differentiability at $x=0$ only assuming the
$a_{ij}$ are square-Dini continuous at $x=0$. More generally, if we let
$a^{0}_{ij}(r)$ denote the mean of $a_{ij}$ over the sphere $|x|=r$, then we
can uniquely write $a_{ij}(x)=a^{0}_{ij}(|x|)+a^{1}_{ij}(x)$; if
$a^{1}_{ij}(x)$ is Dini continuous and $a_{ij}^{0}(|x|)$ is only square-Dini
continuous at $x=0$, then every solution is Lipschitz continuous and
differentiable at $x=0$.

###  1.3. From a Point to a Neighborhood

Theorem 1 above concerns the Lipschitz continuity and differentiability at
$x=0$, but under certain conditions, these regularity properties extend
locally. Let us consider a simple case where the coefficients are almost
everywhere differentiable with an estimate for the growth of $\nabla
a_{ij}(x)$ as $|x|\to 0$ and obtain a uniform bound near $x=0$ on the gradient
of a weak solution; this will also show the solution is $C^{1}$ near $x=0$.

######  Theorem 2.

Suppose the $a_{ij}$ are differentiable in $\Omega\backslash\\{0\\}$ with

(15) $|\nabla a_{ij}(x)|\leq\frac{c}{|x|}\quad\hbox{for $0<|x|<2\rho$.}$

Suppose the $a_{ij}$ are square-Dini continuous at $x=0$ and the dynamical
system (8) is uniformly stable. Then every weak solution $u\in
H^{1,2}(\Omega)$ of (1) is Lipschitz continuous in $B_{\rho}$ and there is a
positive constant $c$ such that

(16) $|\nabla
u(x)|\leq\frac{c}{\rho}\left({\diagup\hskip-11.9501pt\int}_{|y|<2\rho}|u(y)|^{2}\,dy\right)^{1/2}\quad\hbox{for
$0<|x|<\rho$.}$

Moreover, $u\in C^{1}(\overline{B_{\rho}})$.

Proof. We will require standard elliptic interior estimates in an annulus. For
$\rho>0$, let us introduce the annulus

(17) $A_{r}=\\{x\in{\mathbb{R}}^{n}:\frac{r}{2}<|x|<2r\\}.$

We first consider $r=1$. We know $|\nabla a_{ij}(x)|\leq c$ for $x\in A_{1}$,
so if $A_{1}^{\prime}$ is a concentric open sub-annulus of $A_{1}$ such that
$\overline{A_{1}^{\prime}}\subset A_{1}$ and $1<p<\infty$, then there is a
positive constant $c$ such that

(18) $\|u\|_{H^{2,p}(A_{1}^{\prime})}\leq c\,\|u\|_{L^{1}(A_{1})}.$

(The above can be obtained from the standard $L^{p}$ elliptic estimates: cf.
Theorem 15.1” and Remark 1 following Theorem 7.3 in [1]; also (14.3.6) in
Chapter 14 of [32].) Taking $p>n$ so that we can use the Sobolev embedding
$H^{2,p}(A_{1}^{\prime})\subset C^{1}(\overline{A_{1}^{\prime}})$, we obtain
the following:

(19) $|\nabla u(x)|\leq c\,\|u\|_{L^{1}(A_{1})}\quad\hbox{for}\
\frac{3}{4}\leq|x|\leq\frac{3}{2}.$

Applying a dilation to the annulus, we find that if

(20) $|\nabla a_{ij}(x)|\leq\frac{c}{r}\quad\hbox{for $x\in A_{r}$,}$

where $c$ is independent of $r$, then there is a constant $c$ independent of
$r$ such that

(21) $|\nabla
u(x)|\leq\frac{c}{r}\,{\diagup\hskip-11.9501pt\int}_{A_{r}}|u(y)|\,dy\quad\hbox{for
$\frac{3}{4}r\leq|x|\leq\frac{3}{2}r$.}$

For $x\in B_{\rho}\backslash\\{0\\}$, choose $r>0$ so that
$\frac{3}{4}r\leq|x|\leq\frac{3}{2}r$. Since $u(x)-u(0)$ is also a solution of
(1), we can apply (21) to it and then use (11):

(22) $\displaystyle|\nabla u(x)|$
$\displaystyle\leq\frac{c}{r}\,{\diagup\hskip-11.9501pt\int}_{A_{r}}|u(y)-u(0)|\,dy\leq\frac{c}{r}\,{\diagup\hskip-11.9501pt\int}_{A_{r}}\frac{c|y|}{\rho}dy\left({\diagup\hskip-11.9501pt\int}_{|z|<2\rho}|u(z)|^{2}\,dz\right)^{1/2}$
$\displaystyle\leq\frac{c}{\rho}\left({\diagup\hskip-11.9501pt\int}_{|z|<2\rho}|u(z)|^{2}\,dz\right)^{1/2},$

where we have used $|y|<2r$ for $y\in A_{r}$.

The uniform bound (16) also enables us to conclude that $u\in
C^{1}(\overline{B_{\rho}})$. In fact, we can use a mollifier to obtain
coefficients $a_{ij}^{h}(x)$ that are smooth on $B_{2\rho}$ and $a_{ij}^{h}\to
a_{ij}$ in $L^{p}(B_{2\rho})$ as $h\to 0$ for $1\leq p<\infty$. If we let
$u^{h}$ denote the solution of (1) in $B_{2\rho}$ with coefficients
$a_{ij}^{h}$ and Dirichlet condition $u^{h}=u$ on $\partial B_{2\rho}$, then
we know that $u^{h}\in C^{\infty}(B_{2\rho})$ and $u^{h}\to u$ in
$H^{1,2}(B_{2\rho})$ (see pages 150-151 in [27] or Section 4 of Chapter 2 in
[25]). Moreover the continuity modulus of the $a_{ij}^{h}$ does not increase,
and the estimate (15) continues to hold with a constant $c$ independent of $h$
near zero, so the estimate (16) applies to $u^{h}$ with $c_{1}$ independent of
$h$ near zero. If we consider a sequence $v^{m}=u^{h_{m}}$ with $h_{m}\to 0$
as $m\to\infty$, we know $\|v^{\ell}-v^{m}\|_{L^{2}(B_{2\rho})}\to 0$ as
$\ell,m\to\infty$ and can use (16) to conclude
$\|\nabla(v^{\ell}-v^{m})\|_{L^{\infty}(B_{\rho})}\to 0$. By completeness of
the Banach space $C^{1}(\overline{B_{\rho}})$, this shows $u\in
C^{1}(\overline{B_{\rho}})$. $\Box$

######  Remark 1.

Instead of assuming differentiability of the coefficients for $x\not=0$ in
Theorem 1, we could have assumed they satisfy
$|a_{ij}(x)-a_{ij}(y)|\leq\omega(|x-y|/|x|)$ for $|x-y|<\frac{1}{2}|x|$, where
$\omega$ satisfies the Dini condition (2). Let us also mention the recent work
of De Filippis & Mingione [10] where (15) and our additional conditions on the
coefficients at $x=0$ are replaced by the gradient of the coefficients
belonging to the Lorentz space $L(n,1)$.

##  2\. Gilbarg-Serrin Equations

In [17] and [36], Gilbarg and Serrin considered coefficients of the form

(23) $a_{ij}(r,\theta)=\delta_{ij}+g(r)\,\theta_{i}\theta_{j},\quad\hbox{where
$r=|x|$ and $\theta_{i}=x_{i}/r$.}$

They assumed $g(r)$ is a bounded function satisfying $g(r)>-1+\varepsilon$,
which guarantees uniform ellipticity, and studied solutions of equations
involving the associated operators in both divergence and nondivergence form.
In this section, we consider conditions on $g(r)$ that guarantee that any weak
solution of (1) must be Lipschitz continuous or differentiable at $x=0$, or
when there exists a weak solution that is not Lipschitz continuous at $x=0$.
In cases where the solution is Lipschitz continuous in a neighborhood of
$x=0$, we will use Theorem 2 to conclude that the solution is $C^{1}$. Let us
state several results in this direction.

###  2.1. Application of the Dynamical Systems Method

Let us apply the results of the dynamical systems method to coefficients of
the form (23). It is shown in Section 5 of [28] that the dynamical system (8)
reduces in this case to the (scalar) ordinary differential equation

(24) $\frac{d\phi}{dt}=\frac{n-1}{n}\,{\rm g}(t)\,\phi,$

where ${\rm g}(t)=g(e^{-t})$. We can integrate this ODE to see that it is
uniformly stable if and only if there is a real constant $K$ such that

$\int_{s}^{t}{\rm g}(\tau)\,d\tau\leq K<\infty\quad\hbox{for $t>s>T$,}$

and all solutions are asymptotically constant if and only if the improper
integral

$\int_{T}^{\infty}{\rm g}(\tau)\,d\tau\quad\hbox{converges to an extended real
number $<\infty$.}$

We are now able to combine Theorems 1 and 2 to obtain:

######  Proposition 1.

Let the coefficients $a_{ij}$ be given by (23), where $g(r)$ is square-Dini
continuous at $0$, and let $u\in H^{1,2}(\Omega)$ be a weak solution of (1).
If, for some real constant $K$, we have

(25) $\int_{r_{1}}^{r_{2}}\frac{g(r)}{r}\,dr\leq K<\infty\quad\hbox{for \
$0<r_{1}<r_{2}<\rho$,}$

then $u$ is Lipschitz continuous at $x=0$ and satisfies (11). If, in addition,

(26) $\int_{0}^{1}\frac{g(r)}{r}\,dr\quad\hbox{converges to an extended real
number $<\infty$,}$

then $u$ is differentiable at $x=0$ and satisfies (12). On the other hand, if
we just assume (25) and that $g(r)$ is differentiable for $r>0$ with
derivative satisfying $|r\,g^{\prime}(r)|\leq c$ for $0<r<\rho$, then $u\in
C^{1}(\overline{B_{\rho}})$ and satisfies the estimate (16).

######  Corollary 1.

If $g(r)$ is square-Dini continuous at $r=0$ and satisfies both (25) and

(27) $\int_{\varepsilon}^{1}\frac{g(r)}{r}\,dr\to-\infty\quad\hbox{as}\
\varepsilon\to 0,$

then every weak solution of (1) is differentiable at $r=0$ with $\nabla
u(0)=0$.

###  2.2. Application of a Comparison Principle

We will discuss a comparison principle based upon solutions of the following
ordinary differential equation:

(28)
$\frac{1}{r^{n-1}}\frac{d}{dr}\left(r^{n-1}(1+g(r))\frac{dv}{dr}\right)-\frac{n-1}{r^{2}}\,v=0,\quad\hbox{for\
$0<r<1$.}$

######  Proposition 2.

Let $Z(r)$ denote a solution of (28) satisfying the finite energy condition

(29)
$\int_{0}^{1}\left(|Z^{\prime}(r)|^{2}+r^{-2}|Z(r)|^{2}\right)r^{n-1}\,dr<\infty.$

If $Z(r)$ satisfies $0<Z(r)\leq c\,r$ as $r\to 0$, then every weak solution
$u\in H^{1,2}(\Omega)$ of (1) with coefficients given by (23) is Lipschitz
continuous at $x=0$. In addition, if $g(r)$ in (23) is differentiable for
$r>0$ and $|rg^{\prime}(r)|\leq c$ for $0<r<2\rho$, then $u\in
C^{1}(\overline{B_{\rho}})$ and satisfies the estimate (16).

Proof. We let $Z(r)$ be a solution of (28) satisfying $Z(1)>0$ and the finite
energy condition (29). By the non-oscillatory character of (28), $Z(r)$ cannot
change sign for $0<r<1$, so we have $Z(r)>0$. With coefficients $a_{ij}$ given
by (23), the equation (1) can be written as

(30) $\frac{1}{r^{n-1}}\frac{\partial}{\partial
r}\left(r^{n-1}(1+g(r))\frac{\partial u}{\partial
r}\right)+\frac{1}{r^{2}}\Delta_{\theta}u=0,$

where $\Delta_{\theta}$ denotes the Laplace-Beltrami operator on the sphere
$S^{n-1}$. If $u\in H^{1,2}(B_{1})$ is a weak solution, we can integrate this
equation over $S^{n-1}$ to find that $u$ satisfies

${\diagup\hskip-11.9501pt\int}_{S^{n-1}}u(r\,\theta)\,d\theta=constant=u(0),$

where $u(0)$ denotes the value of $u(x)$ at $x=0$. But since we can subtract a
constant from a weak solution of (1), we may assume that $u(0)=0$, so
$u(r,\theta)$ is orthogonal to 1 on every sphere centered at $0$:

(31) $\int_{S^{n-1}}u(r\,\theta)\,d\theta=0\quad\hbox{for $0<r<1$},$

We use the following comparison principle between $Z$ and $u$:

######  Lemma 1.

Let the coefficients $a_{ij}$ be given by (23) and $u\in H^{1,2}(B_{1})$ is a
weak solution of (1) satisfying (31). Then

(32)
$\left({\diagup\hskip-11.9501pt\int}_{S^{n-1}}|u(r\,\theta)|^{2}\,d\theta\right)^{1/2}\leq\frac{Z(r)}{Z(1)}\left({\diagup\hskip-11.9501pt\int}_{S^{n-1}}|u(\theta)|^{2}\,d\theta\right)^{1/2}\quad\hbox{for\
$0<r<1$.}$

We will prove Lemma 1 below. But to complete the proof of Proposition 2,
recall from Moser [34] that a weak solution of (1) satisfies the estimate

(33) $\sup_{|x|<\delta}|u(x)|\leq
c\left({\diagup\hskip-11.9501pt\int}_{B_{2\delta}}|u(y)|^{2}\,dy\right)^{1/2}\quad\hbox{for}\
0<\delta\leq 1/2.$

We may combine this with (32) and the assumption $Z(r)\leq c\,r$ to conclude
that $u(x)$ is Lipschitz continuous at $x=0$. If $g$ is differentiable with
$|rg^{\prime}(r)|\leq c$ for $0<r<2\rho$, then we can use Theorem 2 to
conclude $u\in C^{1}(\overline{B_{\rho}})$ and obtain the estimate (16).
$\Box$

######  Remark 2.

The differentiability at $x=0$ of the solution $Z(r)\Theta_{1}(\theta)$ of (1)
means $Z(r)\Theta_{1}(\theta)=V\cdot x+o(r)$ where $V$ is a constant vector.
Multiplying both sides by $\Theta_{1}$ and integrating over the unit sphere we
obtain $Z(r)=c\,r+o(r)$ for some constant $c$.

Proof of Lemma 1. Considering $u$ as a weak solution of (30) in the ball
$B_{\rho}$ for $0<\rho<1$, we have

(34) $\int_{B_{\rho}(0)}\left((1+g(r))\left|\frac{\partial u}{\partial
r}\right|^{2}+\frac{1}{r^{2}}\left|\nabla_{\theta}u\right|^{2}\right)\,dx=(1+g(\rho))\rho^{n-1}\int_{S^{n-1}}u(\rho\,\theta)\frac{\partial
u}{\partial r}(\rho\,\theta)\,d\theta.$

Using (31), minimizing the integration over the sphere $S^{n-1}$ on the left
hand side, and rewriting the expression on the right hand side, we find

(35) $\int_{S^{n-1}}\\!\int_{0}^{\rho}\left((1+g(r))\left|\frac{\partial
u}{\partial
r}\right|^{2}+\frac{n-1}{r^{2}}\left|u\right|^{2}\right)r^{n-1}drd\theta\leq(1+g(\rho))\frac{\rho^{n-1}}{2}\frac{d}{d\rho}\int_{S^{n-1}}(u(\rho\,\theta))^{2}d\theta.$

Now consider the variational functional associated with (28):

(36a)
$J(v)=\int_{0}^{\rho}\left((1+g(r))|v^{\prime}(r)|^{2}+\frac{n-1}{r^{2}}|v(r)|^{2}\right)r^{n-1}\,dr$
with the boundary condition (36b)
$v(\rho)=\left({\diagup\hskip-11.9501pt\int}_{S^{n-1}}|u(\rho\,\theta)|^{2}\,d\theta\right)^{1/2}.$

If we let $v_{*}$ denote the minimizer for (36), then we can express it as a
scalar multiple of $Z(r)$ by

(37)
$v_{*}(r)=\frac{Z(r)}{Z(\rho)}\left({\diagup\hskip-11.9501pt\int}_{S^{n-1}}|u(\rho\,\theta)|^{2}\,d\theta\right)^{1/2},$

and the minimum value of $J$ is

(38)
$J(v_{*})=(1+g(\rho))\rho^{n-1}\frac{Z^{\prime}(\rho)}{Z(\rho)}{\diagup\hskip-11.9501pt\int}_{S^{n-1}}|u(\rho\,\theta)|^{2}\,d\theta.$

Now if we fix $\theta\in S^{n-1}$ and evaluate $J$ on $v(r)=u(r\theta)$, we
obtain

$\int_{0}^{\rho}\left((1+g(r))\left|\frac{\partial v}{\partial
r}\right|^{2}+\frac{n-1}{r^{2}}\left|v\right|^{2}\right)\,r^{n-1}dr\geq
J(v_{*}).$

If we integrate both sides over $\theta\in S^{n-1}$ and use (35) and (38), we
obtain

(39)
$\frac{1}{2}\frac{d}{d\rho}\|u(\rho)\|^{2}\geq\frac{Z^{\prime}(\rho)}{Z(\rho)}\|u(\rho)\|^{2}\quad\hbox{where}\
\|u(\rho)\|^{2}:=\int_{S^{n-1}}|u(\rho\,\theta)|^{2}\,d\theta.$

But (39) implies

$\frac{d}{d\rho}\log\frac{\|u(\rho)\|}{Z(\rho)}\geq 0,$

which implies that $\|u(\rho)\|/Z(\rho)$ is increasing, so (32) follows.
$\Box$

###  2.3. Conditions for Loss of Lipschitz Regularity

Propositions 1 and 2 provide conditions under which Lipschitz regularity
holds. Now let us consider conditions under which a weak solution exists which
is not Lipschitz continuous.

######  Proposition 3.

Let the coefficients $a_{ij}$ be given by (23), where $g(r)$ is square-Dini
continuous at $0$ and

(40) $\limsup_{r\to 0}\int_{r}^{1/2}\frac{g(\rho)}{\rho}\,d\rho=+\infty.$

Also assume $g(r)$ has finite total variation, i.e.
$\int_{0}^{1/2}|g^{\prime}(r)|\,dr<\infty$. Then there exists a weak solution
$u$ of (1) in $B_{\rho}$ that is not Lipschitz continuous at $x=0$. If we also
have $g(r)>0$ for $r<\varepsilon$, then we do not require $g(r)$ to be square-
Dini continuous at $r=0$ for the conclusion to hold.

Proof. We want to find a weak solution $u\in H^{1,2}(B_{\rho})$ of (1) that is
bounded but not Lipschitz continuous. Without loss of generality, we may
assume $\rho=1/2$ and $u(0)=0$. We shall look for a solution in the form
$u(r,\theta)=v(r)\,\Theta_{1}(\theta)$, where $v(0)=0$ and $\Theta_{1}$ is an
eigenfunction for the first eigenvalue $\lambda=n-1$ for the Laplace-Beltrami
operator on the unit sphere. We find that $v(r)$ should satisfy the ordinary
differential equation (28). With a change of independent variable to $t=-\log
r$, the equation for $v$ can be written as

(41) $\frac{d}{dt}\left((1+{\rm g})\frac{dv}{dt}\right)-(n-2)(1+{\rm
g})\frac{dv}{dt}-(n-1)v=0,\quad\hbox{for $0<t<\infty$,}$

where ${\rm g}(t)=g(e^{-t})$, and this can be written as a 1st-order system

(42) $\frac{d}{dt}\begin{bmatrix}v\\\ w\end{bmatrix}=\begin{bmatrix}0&(1+{\rm
g})^{-1}\\\ n-1&n-2\end{bmatrix}\begin{bmatrix}v\\\
w\end{bmatrix}\quad\hbox{where}\ w=(1+{\rm g})\frac{dv}{dt}.$

Let $M_{\rm g}$ denote the $2\times 2$ matrix on the right hand side of (42).
For ${\rm g}=0$, $M_{0}$ has eigenvalues $\lambda_{1}=-1$ and
$\lambda_{2}=n-1$ (corresponding to the solutions $r$ and $r^{1-n}$ of (28)
when $g=0$). If we let $E$ denote the matrix of eigenvectors for $M_{0}$, then
we can rewrite (42) as

(43a) $\frac{d}{dt}V=(A_{0}+B({\rm g}))\,V,\quad\hbox{where}\
\begin{bmatrix}v\\\ w\end{bmatrix}=E\,V$ and (43b)
$A_{0}=\begin{bmatrix}-1&0\\\ 0&n-1\end{bmatrix},\quad B({\rm g})=\frac{\rm
g}{n(1+{\rm g})}\begin{bmatrix}n-1&-(n-1)^{2}\\\ 1&-(n-1)\end{bmatrix}.$

The eigenvalues of $A_{0}+B({\rm g})$ may be denoted $\Lambda_{1}({\rm g})$
and $\Lambda_{2}({\rm g})$, where $\Lambda_{1}({\rm g})\to-1$ and
$\Lambda_{2}({\rm g})\to n-1$ as ${\rm g}\to 0$. We can calculate
$\Lambda^{\prime}:=d\Lambda/d{\rm g}$ and evaluate $\Lambda^{\prime}(0)$ to
find

(44) $\Lambda_{1}^{\prime}(0)=\frac{n-1}{n}.$

Consequently

(45) $\Lambda_{1}({\rm g})=-1+\frac{n-1}{n}{\rm g}+O(|{\rm
g}|^{2})\quad\hbox{as}\ {\rm g}\to 0.$

Similarly, we can use $\Lambda(0)=n-1$ and obtain

(46) $\Lambda_{2}({\rm g})=n-1+\frac{n-1}{n}{\rm g}+O(|{\rm
g}|^{2})\quad\hbox{as}\ {\rm g}\to 0.$

At this point we have not used the fact that $g(r)$ is square-Dini continuous
at $r=0$, so we do not assert that $\int_{0}^{\infty}|{\rm
g}(t)|^{2}\,dt<\infty$. Nevertheless, because we have assumed that $g(r)$ has
finite total variation for $0<r<1/2$, we know that the matrix function
$B(t):=B({\rm g}(t))$ has finite total variation for $t>\log 2$, so we can use
the asymptotic theory of linear systems of ordinary differential equations.
Denoting the eigenvalues as $\lambda_{1}(t):=\Lambda_{1}({\rm g}(t))$ and
$\lambda_{2}(t):=\Lambda_{2}({\rm g}(t))$, we have

(47) $\lambda_{1}(t)=-1+\frac{n-1}{n}\,{\rm g}(t)+O(|{\rm
g}(t)|^{2})\quad\hbox{as}\ t\to\infty.$

We can apply Theorem 8.1 in [9] (see also Theorem 16.3.1 in [21]) to conclude
that there exists a solution $V_{1}(t)$ with the asymptotic behavior

(48)
$\lim_{t\to\infty}V_{1}(t)\,\exp\left[-\int_{0}^{t}\lambda_{1}(\tau)\,d\tau\right]=V_{0}\,\quad\hbox{as}\
t\to\infty,$

where $V_{0}=(1,0)^{t}$ is the eigenvector for the eigenvalue
$\Lambda_{1}(0)=-1$.

Now we shall invoke the sqaure-Dini assumption $\int_{0}^{\infty}|{\rm
g}(t)|^{2}\,dt<\infty$. This mean that (43) can be written

(49a) $\frac{d}{dt}V=(A_{0}+{\rm g}(t)\,B_{0})\,V+C({\rm g}(t))\,V,$ where
(49b) $A_{0}=\begin{bmatrix}-1&0\\\ 0&n-1\end{bmatrix},\quad
B_{0}:=\frac{1}{n}\begin{bmatrix}n-1&-(n-1)^{2}\\\
1&-(n-1)\end{bmatrix},\quad\hbox{and}\quad\int_{0}^{\infty}|C({\rm
g}(t))|\,dt<\infty.$

Now Theorem 8.1 in [9] shows that (48) can be sharpened somewhat to

(50) $V_{1}(t)=V_{0}\,\exp\left[-t+\frac{n-1}{n}\int_{0}^{t}{\rm
g}(\tau)\,d\tau\right](1+o(1))\quad\hbox{as}\ t\to\infty.$

If we transform back to our original dependent variables $v,w$ using
$(v,w)^{t}=E^{-1}V$ and return to the original independent variable
$r=e^{-t}$, then we find there is a solution $v(r)$ of (28) with the
asymptotic behavior

(51)
$v(r)=c\,r\exp\left[\frac{n-1}{n}\int_{r}^{1/2}\frac{g(\rho)}{\rho}\,d\rho\right](1+o(1))\quad\hbox{as}\
r\to 0.$

But the condition (40) shows that this solution of (28) is not Lipschitz at
$r=0$. In other words, we have found a bounded solution
$u(r,\theta)=v(r)\Theta_{1}(\theta)$ of (1) that is not Lipschitz at $r=0$.

If we do not assume that $g(r)$ is square-Dini at $r=0$, but we know $g(r)>0$
for $r<\varepsilon$, then we still can conclude

$\frac{n-1}{n}\int_{0}^{t}\left({\rm g}(\tau)+O(|{\rm
g}(\tau)|^{2})\right)\,d\tau\to+\infty\quad\hbox{as}\ t\to\infty,$

which shows that our solution $v(r)$ is not Lipschitz continuous at $r=0$.
$\Box$

######  Remark 3.

Propositions 2 and 3 both depend on the asymptotic behavior of the positive
finite energy solution $Z(r)$ of (28) in such a way that $Z(r)\leq c\,r$
becomes a necessary and sufficient condition for every weak solution of (1) to
be Lipschitz continuous at $x=0$.

######  Remark 4.

As we shall see in the following Examples, it is possible to find functions
$g(r)$ that do not have Dini mean oscillation, yet satisfy the conditions of
Proposition 1 or 2 that show all weak solutions are $C^{1}$ in a neighborhood
of $x=0$. Proposition 3, on the other hand, can be used to find examples of
coefficients $a_{ij}$ that do not satisfy the Dini mean oscillation condition:
if $g(r)$ satisfies the square-Dini condition and (40), then (1) has weak
solutions $u\in H^{1,2}(\Omega)$ that are not Lipschitz continuous, so by the
results of [15] the coefficients (23) cannot have Dini mean oscillation.

###  2.4. Examples

In this section we consider several specific functions $g(r)$ and what our
preceding results tell us about Lipschitz regularity at $x=0$ or in a
neighborhood of $x=0$. We frequently find that Proposition 2 enables us to
extend results beyond what Proposition 1 is able to provide due to the
condition that $g(r)$ be square-Dini continuous at $r=0$. But both
Propositions 1 and 2 apply in some cases that do not have Dini mean
oscillation.

Example 1. Let us consider (23) with

(52) $g(r)=|\log r|^{-\gamma}\quad\hbox{for $\gamma>0$.}$

For $\gamma>1$, $g(r)$ is Dini continuous at $r=0$, so we know weak solutions
of (1) are Lipschitz continuous (in fact, $C^{1}$). For $0<\gamma\leq 1$,
$g(r)$ is not Dini continuous: in fact

$\int_{0}^{1/2}\frac{g(r)}{r}\,dr=\int_{0}^{1/2}\frac{|\log
r|^{-\gamma}}{r}\,dr=+\infty,$

so condition (40) holds. Moreover, $g(r)$ has finite total variation:

$\int_{0}^{1/2}|g^{\prime}(r)|\,dr=\gamma\int_{0}^{1/2}\frac{|\log
r|^{-\gamma-1}}{r}\,dr<\infty.$

Now for $1/2<\gamma\leq 1$, $g(r)$ is square-Dini continuous, so we can use
Proposition 3 to conclude there exists a weak solution of (1) that is not
Lipschitz continuous at $x=0$. In fact, since $g(r)>0$ for $0<r<1/2$ we can
use Proposition 3 to extend this conclusion to $0<\gamma\leq 1$ . Summary for
(52): For $\gamma>1$, all weak solutions of (1) are Lipschitz continuous; but
for $0<\gamma\leq 1$, there exists a weak solution of (1) that is not
Lipschitz continuous at $x=0$.

Now let us consider (23) with

(53) $g(r)=-|\log r|^{-\gamma}\quad\hbox{for $\gamma>0$.}$

For $\gamma>1$, as with (52), $g(r)$ is Dini continuous at $r=0$, so we know
weak solutions are Lipschitz continuous; but for all $\gamma>0$ we have
$|rg^{\prime}(r)|\leq c$ and

$-\infty\leq\int_{0}^{1/2}\frac{g(r)}{r}\,dr=-\int_{0}^{1/2}\frac{|\log
r|^{-\gamma}}{r}<0.$

For $1/2<\gamma\leq 1$, $g(r)$ is square-Dini continuous and the condition
(25) is trivially satisfied with $K=0$, so we could use Proposition 1 to
conclude that all weak solutions of (1) are Lipschitz continuous. However, we
can instead use Proposition 2 to handle all values $0<\gamma\leq 1$. We need
to find the solution $Z(r)$ of (28) with the desired properties. From (47) and
(48) we see that there is a solution $V(t)$ of (43) such that
$|V(t)|=o(e^{-t})$, which means we have our solution $v(t)=Z(r)$ satisfying
$|Z(r)|\leq c\,r$. If we rescale so that $Z(1)=1$, then the non-oscillatory
character of (28) shows $Z(r)$ remains positive. We can also confirm that $Z$
satisfies the finite energy condition (29), so we can apply Proposition 2 to
reach our conclusion. Summary for (53): For $\gamma>0$, all weak solutions of
(1) are Lipschitz continuous.

Note: For $0<\gamma\leq 1$, these results for (52) and (53) cannot be obtained
using Dini mean oscillation for $a_{ij}$ because that condition is satisfied
only for $\gamma>1$: this can be confirmed by direct calculation at $x=0$ (see
the Appendix), but for (52) it also follows from Remark 4.

Example 2. Let us consider an example where $g(r)$ oscillates as $r\to 0$:

(54) $g(r)=\frac{\sin(|\log r|)}{|\log r|^{\beta}},\quad\hbox{where
$\beta>0$.}$

The condition $|rg^{\prime}(r)|\leq c$ in Propositions 1 and 2 is satisfied.
Also note that

(55)
$\int_{0}^{e^{-1}}\frac{g(r)}{r}\,dr=\int_{1}^{\infty}\frac{\sin\tau}{\tau^{\beta}}\,d\tau\quad\hbox{converges}$

since $\tau^{-\beta}$ is strictly decreasing and $\sin\tau$ alternates sign on
intervals $[k\pi,(k+1)\pi]$ for $k=0,1,\dots$. In fact, for any interval
$[t_{1},t_{2}]\subset(2\pi,\infty)$ we have

(56)
$\int_{t_{1}}^{t_{2}}\frac{\sin\tau}{\tau^{\beta}}\,d\tau\leq\int_{2\pi}^{3\pi}\frac{\sin\tau}{\tau^{\beta}}\,d\tau<\infty.$

For $\beta>1/2$, this $g(r)$ is square-Dini continuous and (56) confirms
condition (25), so by Proposition 1 we conclude that all weak solutions of (1)
are Lipschitz continuous. To obtain this conclusion for all $\beta>0$, let us
appeal to Proposition 2. Using asymptotic analysis as above with (55), we can
find a finite energy solution $Z(t)$ of (28) satisfying $0<Z(r)\leq c\,r$ for
$0<r<1$, so Proposition 2 shows that all weak solutions are Lipschitz
continuous.

Note: These results for $0<\beta\leq 1$ in (54) cannot be obtained using Dini
mean oscillation since one can show by direct calculation (see the Appendix)
that the coefficients $a_{ij}$ have Dini mean oscillation only for $\beta>1$.

Example 3. Consider the following function:

(57) $g(r)=\frac{-C_{1}\,\sin(|\log r|)-C_{2}\cos(|\log r|)}{A+\sin(|\log
r|)-\cos(|\log r|)}$

where

$C_{1}=\frac{(n-1)^{2}}{(n-1)^{2}+1}+1,\quad C_{2}=\frac{n-1}{(n-1)^{2}+1}-1,$

and $A>1$ is chosen sufficiently large. This complicated expression for $g(r)$
arises because we actually want the following function $Z(r)$ to be a solution
of (28):

(58) $Z(r)=r\,\left(A+\sin(|\log r|)\right).$

Using $Z(r)$ in (28), we can work backwards to find $g(r)$ and confirm (57).
This means that, if we use this $g(r)$ in our coefficients $a_{ij}$ as in
(23), then by Proposiiton 2 any weak solution $u\in H^{1,2}(B_{1/2}(0))$ of
(1) must be Lipschitz continuous in $B_{1/2}(0)$. This is in spite of the fact
that the coefficients are not square-Dini continuous or Dini mean oscillation
(see Appendix); in fact, $g$ does not even vanish as $r\to 0$, and we claim
that

(59) $\int_{0}^{1/2}\frac{g(r)}{r}\,dr=+\infty.$

This is interesting since the condition (59) is associated with the existence
of a weak solution that is not Lipschitz continuous (cf. Proposition 3).
However, since $g(r)$ does not even vanish as $r\to 0$, any intuition from the
asymptotic analysis used in the proof of Proposition 3 is not relevant for
this example.

To confirm (59), let us compute

$\int_{0}^{1}\frac{g(r)}{r}\,dr=\int_{0}^{\infty}{\rm
g}(t)\,dt=\int_{0}^{\infty}\frac{-C_{1}\,\sin t-C_{2}\cos t}{A+\sin t-\cos
t}\,dt$

By periodicity, it suffices to compute

$\begin{aligned} \int_{0}^{2\pi}\frac{-C_{1}\,\sin t-C_{2}\cos t}{A+\sin
t-\cos t}\,dt&=\frac{1}{A}\int_{0}^{2\pi}\frac{-C_{1}\,\sin t-C_{2}\cos
t}{1+\frac{1}{A}(\sin t-\cos t)}\,dt\\\ &=A^{-2}\int_{0}^{2\pi}(\sin t-\cos
t)\left[C_{1}\sin t+C_{2}\cos t\right]\,dt+O(A^{-3})\end{aligned},$

provided $A$ is sufficiently large. By orthogonality, the cross terms vanish
and we are left with computing

$\int_{0}^{2\pi}\left[C_{1}\sin^{2}t-C_{2}\cos^{2}t\right]\,dt=\int_{0}^{2\pi}\frac{C_{1}-C_{2}}{2}\,dt,$

where we have used half-angle trigonometric identities for $\sin^{2}$ and
$\cos^{2}$ as well as $\int_{0}^{2\pi}\cos 2t\,dt=0=\int_{0}^{2\pi}\sin
2t\,dt$. But it is clear that $C_{1}-C_{2}>0$, so we obtain

$\int_{0}^{\infty}\frac{-C_{1}\,\sin t-C_{2}\cos t}{A+\sin t-\cos
t}\,dt=+\infty.$

##  Appendix A Computations of mean oscillation at $x=0$

All of the functions that we consider are smooth except at $x=0$, so we will
restrict our attention to the mean oscillation at $x=0$, i.e. we take the sup
over $B_{r}=\\{y:|y|<r\\}$ in (3). However, if we want to show for certain
examples that Dini mean oscillation over a domain fails, it suffices to show
that it fails at $x=0$.

First let us confirm that $a(x)=|\log r|^{-\gamma}$ has Dini mean oscillation
at $x=0$ for all $\gamma>0$. Since $a(x)=a(r)$ is radial, we easily compute

(60) $\displaystyle\widetilde{a}(r):=\widetilde{a}_{0,r}$
$\displaystyle={\diagup\hskip-11.9501pt\int}_{B_{r}}a(x)\,dx=\frac{n}{r^{n}}\int_{0}^{r}|\log\rho|^{-\gamma}\,\rho^{n-1}\,d\rho$
$\displaystyle=|\log r|^{-\gamma}-\gamma\,|\log
r|^{-\gamma-1}+\frac{\gamma(\gamma+1)n}{r^{n}}\int_{0}^{r}|\log\rho|^{-\gamma-2}\,\rho^{n-1}\,d\rho.$

So $a(x)-\widetilde{a}(r)=\gamma\,|\log r|^{-\gamma-1}+O(|\log
r|^{-\gamma-2})$ and hence the mean oscillation is

(61)
$\omega_{a}(r)={\diagup\hskip-11.9501pt\int}_{B_{r}}|a(x)-\widetilde{a}(r)|\,dx=\gamma\,|\log
r|^{-\gamma-1}+O(|\log r|^{-\gamma-2}).$

Since $\omega_{a}(r)$ satisfies (2), this confirms that $a(x)$ has Dini mean
oscillation at $x=0$ for all $\gamma>0$.

In the above calculation, there was no angular oscillation, and this
contributed to $a$ satisfying the Dini mean oscillation condition for
$\gamma>0$. The situation is quite different for functions with angular
dependence. In order to estimate the mean oscillation for the coefficients
$a_{ij}$ in (23), let us introduce the matrix

(62) $\Theta=\left(\theta_{i}\theta_{j}\right)_{i,j=1,\dots,n},$

and denote its mean value over the unit sphere $S^{n-1}$ by
$\widetilde{\Theta}$. Then a calculation shows

(63) $\widetilde{\Theta}=\frac{1}{n}\left(\delta_{ij}\right)_{i,j=1,\dots,n}.$

Now let us multiply the matrix $\Theta$ by a scalar function $g(r)$ that is
smooth for $0<r<1/2$. If we denote by $\widetilde{g\Theta}(r)$ the mean value
of $g(r)\Theta$ over the ball $B_{r}(0)$, then we have

$\widetilde{g\Theta}(r)=\widetilde{g}(r)\,\widetilde{\Theta},\quad\hbox{where}\
\widetilde{g}(r)=\frac{n}{r^{n}}\int_{0}^{r}g(\rho)\rho^{n-1}\,d\rho.$

Now we may estimate the matrix norm of $g\,\Theta-\widetilde{g\,\Theta}$ from
below by

(64)
$\boldsymbol{|}g(r)\Theta-\widetilde{g\,\Theta}(r)\boldsymbol{|}\geq|g(r)|\,\boldsymbol{|}\Theta-\widetilde{\Theta}\boldsymbol{|}-|g(r)-\widetilde{g}(r)|\,\boldsymbol{|}\widetilde{\Theta}\boldsymbol{|}.$

If we can show $|g(r)-\widetilde{g}(r)|=o(|g(r)|)$ as $r\to 0$, then we know
that the matrix $g(r)\Theta$ has Dini mean oscillation at $x=0$ if and only if
the function $|g(r)|$ itself satisfies the Dini condition at $r=0$.

For example, consider the coefficients in Example 1, which we can write in
matrix form as

(65) $A(x)=I+|\log r|^{-\gamma}\Theta,\quad\hbox{where $I$ denotes the
identity matrix}.$

The mean value of $I$ is itself, so we have

$A(x)-\widetilde{A}(r)=g(r)\,\Theta-\widetilde{g}(r)\,\widetilde{\Theta}\quad\hbox{where}\
g(r)=|\log r|^{-\gamma}.$

As we found in (60), $|g(r)-\widetilde{g}(r)|=o(|g(r)|)$ as $r\to 0$, so by
(64), $A(x)$ has Dini mean oscillation at $x=0$ if and only if $g(r)=|\log
r|^{-\gamma}$ satisfies the Dini condition. So Example 1 has Dini mean
oscillation at $x=0$ for $\gamma>1$, but not for $0<\gamma\leq 1$.

Next let us consider the coefficients in Example 2. Again we need only concern
ourselves with the matrix

$g(r)\Theta\quad\hbox{where}\ g(r)=\frac{\sin(|\log r|)}{|\log r|^{\beta}}.$

For $0<r<1$ we use integration by parts and a geometric series to evaluate

$\displaystyle\widetilde{g}(r)$
$\displaystyle=\frac{n}{r^{n-1}}\int_{0}^{r}\frac{\sin(-\log\rho)}{(-\log\rho)^{\beta}}\rho^{n-1}\,d\rho=ne^{nt}\int_{t}^{\infty}\frac{\sin\tau}{\tau^{\beta}}e^{-n\tau}d\tau$
$\displaystyle=\frac{\sin
t}{t^{\beta}}+e^{nt}\int_{t}^{\infty}\frac{\cos\tau}{\tau^{\beta}}e^{-n\tau}d\tau=\frac{\sin
t}{t^{\beta}}+\frac{\cos
t}{nt^{\beta}}-\frac{1}{n}e^{nt}\int_{t}^{\infty}\frac{\sin\tau}{\tau^{\beta}}e^{-n\tau}d\tau$
$\displaystyle=\frac{n}{n^{2}+1}\left(\frac{n\sin t+\cos
t}{t^{\beta}}\right).$

Consequently,

(66) $g(r)-\widetilde{g}(r)=\frac{1}{n^{2}+1}\,\frac{\sin(|\log
r|)-n\cos(|\log r|)}{|\log r|^{\beta}}.$

Now we can estimate

${\diagup\hskip-11.9501pt\int}_{S^{n-1}}\boldsymbol{|}g(r)\Theta-\widetilde{g}(r)\widetilde{\Theta}\boldsymbol{|}\,ds\geq\left|\,{\diagup\hskip-11.9501pt\int}_{S^{n-1}}(g(r)\Theta-\widetilde{g}(r)\widetilde{\Theta})\,ds\,\right|=|g(r)-\widetilde{g}(r)|\,\boldsymbol{|}\widetilde{\Theta}\boldsymbol{|}.$

But if $0<\beta\leq 1$, then

$\int_{0}^{1/2}\frac{|g(r)-\widetilde{g}(r)|}{r}\,dr=\frac{1}{n^{2}+1}\int_{\log
2}^{\infty}\frac{|\sin t-n\cos t|}{t^{\beta}}\,dt=\infty$

since the numerator is periodic. This shows that the coefficients in Example 2
have Dini mean oscillation at $x=0$ only for $\beta>1$.

Finally, we consider the coefficients in Example 3. Recall

$g(r)=\frac{-C_{1}\sin(|\log r|)-C_{2}\cos(|\log r|)}{A+\sin(|\log
r|)-\cos(|\log r|)}=A^{-1}(-C_{1}\sin(|\log r|)-C_{2}\cos(|\log
r|))+O(A^{-2})$

for large $A$. We can use this to approximate $\widetilde{g}(r)$, the mean
value over the ball $B_{r}$, and then express the answer in terms of $t=-\log
r$:

$\widetilde{g}(e^{-t})=-\frac{n}{A(n^{2}+1)}(C_{1}(n\sin t-\cos t)+C_{2}(n\cos
t-\sin t))+O(A^{-2}).$

Hence

${g}(e^{-t})-\widetilde{g}(e^{-t})=-\frac{1}{A(n^{2}+1)}\left[(C_{1}+nC_{2})\sin
t+(C_{2}+nC_{1})\cos t\right]+O(A^{-2}).$

But

$\int_{\log r}^{\infty}|(C_{1}+nC_{2})\sin t+(C_{2}+nC_{1})\cos
t|\,dt=+\infty,$

since the integrand is periodic, so $|g(r)-\widetilde{g}(r)|$ does not satisfy
the Dini condition. Hence Example 3 does not have Dini mean oscillation.

Acknowledgement: This paper has been supported by the RUDN University
Strategic Academic Leadership Program.

## References

  * [1] S. Agmon, A. Douglis and L. Nirenberg, Estimates near the boundary for solutions of elliptic partial dif- ferential equations satisfying general boundary conditions I, Comm. Pure Appl. Math., 12 (1959), 623-727.
  * [2] Yu. A. Alkhutov, Smoothness and limiting properties of solutions of a second-order parabolic equation, Matematicheskie Zametki, 50, 150-152.
  * [3] C. Burch, The Dini condition and regularity of weak solutions of elliptic equations, J. Diff. Eq. 30 (1978), 308-323.
  * [4] F. Chiarenza, M. Frasca and P. Longo, Solvability of the Dirichlet problem for nondivergence elliptic equa- tions with VMO coefficients, Trans. Amer. Math. Soc., 336 (1993), 841–853.
  * [5] A. Cianchi, V. Maz’ya, Global Lipschitz regularity for a class of quasilinear elliptic  Commun. PDE 36, (2011), 100–133.
  * [6] A. Cianchi, V. Maz’ya, Gradient regularity via rearrangements for p-Laplacian type elliptic boundary value problems, Eur. Math. Soc. (JEMS), (2014), 16, p.571-595
  * [7] A. Cianchi, V. Maz’ya, Global boundedness of the gradient for a class of nonlinear elliptic systems, Arch. Ration. Mech. Anal., 212 (2014), p.129-177.
  * [8] A. Cianchi, V. Maz’ya, Global gradient estimates in elliptic problems under minimal data and domain regularity, Commun. Pure Appl. Anal. 14 (2015), no. 1, 285–311.
  * [9] E. Coddington, N. Levinson, Theory of Ordinary Differential Equations, McGraw-Hill, New York, 1955.
  * [10] C. De Filippis, G. Mingione, Lipschitz bounds and nonautonomous integrals, Arch. Rat. Mech. Anal, doi: 10.1007/s00205-021-01698-5
  * [11] E. De Giorgi, Sulla differenziabilità e l’analiticità delle estremali degil integrali multipli regolari, Mem. Accad. Sci. Torino Cl. Sci. Fis. Mat. Natur. (3) 3 (1957), 25-43.
  * [12] G. Di Fazio, Estimates for divergence form elliptic equations with discontinuous coefficients, Boll. Un. Mat. Ital A(7), 10 (1996), 409-420.
  * [13] H. Dong, L. Escauriaza, S. Kim, On $C^{1}$, $C^{2}$, and weak type $(1,1)$ estimates for linear elliptic operators, Part II, Math. Ann. 370 (2018), no. 1-2, 447-489.
  * [14] H. Dong, L. Escauriaza, S. Kim, On $C^{1/2,1}$, $C^{1,2}$, and $C^{0,0}$ estimates for linear parabolic operators, J. Evol. Equ. (2021) https://doi.org/10.1007/s00028-021-00729-8
  * [15] H. Dong, S. Kim, On $C^{1}$, $C^{2}$, and weak type $(1,1)$ estimates for linear elliptic operators, Comm. Partial Differential Equations, 42 (2017), 417-435.
  * [16] H. Dong, J. Lee, S. Kim, On conormal and oblique derivative problem for elliptic equations with Dini mean oscillation coefficients, Indiana Univ. Math. J. 69 (2020), no. 6, 1815–1853.
  * [17] D. Gilbarg, J. Serrin, On isolated singularities of solutions of second-order elliptic equations, J. Analyse Math.  4 (1955/56), 309-340.
  * [18] P. Hartman, A. Wintner, On uniform Dini conditions in the theory of linear partial differential equations of elliptic type, American J. Math. 77 (1955), 329-354.
  * [19] A.M. ll’in, On parabolic equations whose coefficients satisfy Dini’s condition, Mat. Zametki, I (1967), 71-80 = Math. Notes 1 (1967), 46- 51.
  * [20] T. Jin, V.G. Maz’ya, J. Van Schaftingen, Patholgical solutions to elliptic problems in divergence form with continuous coefficients, Comptes rendus. Mathématique, 47 (2009), p.773-778.
  * [21] V.A. Kozlov, V.G. Maz’ya, Differential Equations with Operator Coefficients, Springer-Verlag, New York, 1999.
  * [22] V.A. Kozlov, V.G. Maz’ya, Asymptotic formula for solutions to the Dirichlet problem of elliptic equations with discontinuous coefficients near the boundary, Ann. Scuola Norm. Sup. Pisa Cl. Sci. (5), 2 (2003), 551-600.
  * [23] T. Kuusi, G. Mingione, Universal potential estimates, J. Funct. Anal. 262 (2012), 4205-4269.
  * [24] T. Kuusi, G. Mingione, A nonlinear Stein theorem, Calc. Var. 51 (2014), 45–86.
  * [25] E. M. Landis, Second order equations of elliptic and parabolic type, Translations of Mathematical Monographs, 171. American Mathematical Society, Providence, RI, 1998.
  * [26] M. I. Matıchuk, S. D. Eıdel’man, On parabolic systems with coefficients satisfying Dini’s condition. Dokl. Akad. Nauk SSSR 165 (1965), 482-485.
  * [27] V.G. Maz’ya, On weak solutions of the Dirichlet and Neumann problems, Trans. Moscow Math. Soc. 20 (1969), 135-172.
  * [28] V.G. Maz’ya, R.McOwen, Differentiablilty of solutions to second-order elliptic equations via dynamical systems, J. Differential Equations, 250 (2010), 1137-1168.
  * [29] V.G. Maz’ya, R.McOwen, Second-order differentiability for solutions of elliptic equations in the plane, J. Mathematical Sciences, 191 (2013), 243-253.
  * [30] V.G. Maz’ya, R.McOwen, Differentiability of solutions to the Neumann problem with low-regularity data via dynamical systems, Operator Theory: Advances and Applications, 261 (2017), 343-385.
  * [31] V.G. Maz’ya, M. Mitrea, T. Shaposhnikova, The Dirichlet problem in Lipschitz domains for higher order elliptic systems with rough coefficients, J. Anal. Math. 110 (2010), 167–239.
  * [32] V.G. Maz’ya, T. Shaposhnikova, Theory of Sobolev Multipliers, Springer, 2009.
  * [33] R. McOwen, On elliptic operators in nondivergence and double divergence form, Operator Theory: Advances and Applications, 193 (2009), 159-169.
  * [34] J. Moser, A new proof of De Giorgi’s theorem concerning the regularity problem for elliptic differential equations, Comm. Pure Appl. Math. 13 (1960), 457–468.
  * [35] J. Nash, Continuity of solutions of parabolic and elliptic equations, Amer. J. Math. 80 (1958), 931-954.
  * [36] J. Serrin,  Pathological solutions of elliptic differential equations, Ann. Scuola Norm. Sup. Pisa (3) 18 (1964), 385–387.
  * [37] M. Taylor, Tools for PDE, Mathematical Surveys and Monographs 81, AMS, 2000.


#  LIGO–Virgo correlations between mass ratio and effective inspiral spin:
testing the active galactic nuclei channel

B. McKernan1,2,3,4, K.E.S. Ford1,2,3,4, T. Callister2, W.M. Farr2,5,
R.O’Shaughnessy6, R.Smith7,8, E. Thrane7,8 & A.Vajpeyi7,8
1Department of Astrophysics, American Museum of Natural History, New York, NY
10024, USA
2Center for Computational Astrophysics, Flatiron Institute, New York, NY
10010, USA
3Graduate Center, City University of New York, 365 5th Avenue, New York, NY
10016, USA
4Department of Science, BMCC, City University of New York, New York, NY 10007,
USA
5Department of Physics and Astronomy, Stony Brook University, Stony Brook NY
11794, USA
6 Center for Computational Relativity and Gravitation, Rochester Institute of
Technology, Rochester, NY 14623, USA
7School of Physics and Astronomy, Monash University, Clayton, VIC 3800,
Australia
8 OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery,
Clayton, VIC 3800, Australia
E-mail:bmckernan at amnh.org (BMcK)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Observations by LIGO–Virgo of binary black hole mergers suggest a possible
anti-correlation between black hole mass ratio ($q=m_{2}/m_{1}$) and the
effective inspiral spin parameter $\chi_{\rm eff}$, the mass-weighted spin
projection onto the binary orbital angular momentum (Callister et al., 2021).
We show that such an anti-correlation can naturally occur for binary black
holes assembled in active galactic nuclei (AGN) due to spherical and planar
symmetry-breaking effects. We describe a phenomenological model in which: 1)
heavier black holes live in the AGN disk and tend to spin up into alignment
with the disk; 2) lighter black holes with random spin orientations live in
the nuclear spheroid; 3) the AGN disk is dense enough to rapidly capture a
fraction of the spheroid component. but small in radial extent to limit the
number of bulk disk mergers; 4) migration within the disk is non-uniform,
likely disrupted by feedback from migrators or disk turbulence; 5) dynamical
encounters in the disk are common and preferentially disrupt binaries that are
retrograde around their center of mass, particularly at stalling orbits, or
traps. This model may explain trends in LIGO–Virgo data while offering
falsifiable predictions. Comparisons of predictions in ($q,\chi_{\rm eff}$)
parameter space for the different channels may allow us to distinguish their
fractional contributions to the observed merger rates.

###### keywords:

accretion disks–accretion–galaxies: active –gravitational waves–black hole
physics

††pubyear: 2020††pagerange: LIGO–Virgo correlations between mass ratio and
effective inspiral spin: testing the active galactic nuclei
channel–LABEL:lastpage

##  1 Introduction

Active galactic nuclei (AGN) are powered by the accretion of gas onto a
supermassive black hole (SMBH). A dense population of stars and stellar
remnants is also expected to orbit the SMBH (Morris, 1993; Miralda-Escudé &
Gould, 2000; Hailey et al., 2018; Generozov et al., 2018). A fraction of this
orbiting spheroid population is captured by the AGN over the disk lifetime
(Fabj et al., 2020; MacLeod & Lin, 2020), leading to a large population of
embedded orbiters within AGN disks. The embedded population experiences gas
torques which leads to migration and close dynamical encounters. If binaries
form as a result of these encounters, they can be hardened or softened by gas
torques or further dynamical interactions (e.g. McKernan et al., 2012; Leigh
et al., 2018). Hardening of binaries to merger yields gravitational wave (GW)
signals detectable with LIGO–Virgo (e.g. McKernan et al., 2014; Bartos et al.,
2017; Stone et al., 2017). Since AGN live in deep potential wells, merger
products are easily retained in spite of $\sim 1000\,\mathrm{km\,s^{-1}}$ GW
kicks (Gerosa & Berti, 2019), yielding a mass hierarchy due to repeated
mergers, possibly over multiple AGN episodes.

The properties of AGN are poorly constrained by electromagnetic (EM)
observations, with multiple orders of magnitude variation in density and disk
aspect ratio allowed within common models (Sirko & Goodman, 2003; Thompson et
al., 2005) and orders of magnitude uncertainty in AGN lifetimes (Schawinski et
al., 2015). This necessarily yields a very large range of predicted compact
binary merger rates from this channel (McKernan et al., 2018; Gröbner et al.,
2020). Conversely, GW observations of mergers from this channel allow us to
‘reverse engineer’ properties of both the AGN disk and the nuclear star
cluster (e.g. McKernan et al., 2018).

Detailed simulations of merging populations in the AGN channel include
$N$-body simulations (e.g. Secunda et al., 2019; Secunda et al., 2020b), Monte
Carlo simulations (e.g. Yang et al., 2019; Tagawa et al., 2019; McKernan et
al., 2020a; Tagawa et al., 2020b; McKernan et al., 2020b; Yang et al., 2020;
Tagawa et al., 2021) and hydrodynamic simulations (e.g. Li et al., 2021;
Derdzinski et al., 2021). Broadly, BBH systems can merge in the bulk AGN disk,
outside the disk if ejected by dynamical encounters, and at special locations
in the AGN disk such as a migration trap (Bellovary et al., 2016), where most
migrators arrive from larger disk radii (Secunda et al., 2020a; McKernan et
al., 2020a). The predicted population properties of the AGN channel, including
mass and spin distributions as well as rates are broadly consistent with the
results of O3a (Abbott et al., 2021b), including the relative rates of highly
asymmetric mass ratio events such as GW190814 (Abbott et al., 2020b), or
intermediate mass black hole formation events such as GW190521 (Abbott et al.,
2020a, c). Uniquely, EM counterparts to BBH mergers may occur in this channel
(McKernan et al., 2019; Graham et al., 2020), although confirming candidates
is difficult (Ashton et al., 2020; Palmese et al., 2021). However,
distinguishing between different merger channels is difficult based on mass
and spin distributions alone.

Here we discuss the phenomenological implications for the AGN channel from the
recent claim that there is a negative correlation between the mass ratio
($q\equiv m_{2}/m_{1}$) and the effective inspiral spin parameter $\chi_{\rm
eff}$ found in the LIGO–Virgo O3a data (Callister et al., 2021).

##  2 An anti-correlation between mass ratio and effective inspiral spin
parameter in BBH mergers

With the conclusion of their O3a observing run, the Advanced LIGO (Aasi et
al., 2015) and Virgo (Acernese et al., 2015) experiments have observed
gravitational waves from 48 candidate BBH detections (Abbott et al., 2021a),
with additional candidates reported by independent re-analyses of LIGO–Virgo
data (Venumadhav et al., 2020; Nitz et al., 2021). Using hierarchical Bayesian
inference, Callister et al. (2021) explored the degree of correlation between
BBH mass ratio $q\equiv m_{2}/m_{1}$ and the effective inspiral spin parameter

$\chi_{\mathrm{eff}}=\frac{\vec{\chi}_{1}+q\vec{\chi}_{2}}{1+q}\cdot\hat{L}_{\rm
b},$ (1)

defined as the mass-weighted average of a binary’s dimensionless component
spins $\vec{\chi}_{i}$, projected onto the direction of the binary’s orbital
angular momentum $\hat{L}_{\rm b}$. Modeling the $\chi_{\mathrm{eff}}$
distribution at fixed $q$ as a truncated Gaussian with peak value $\mu$ and
width $\sigma$, Callister et al. (2021) measured the slopes $d\mu/dq$ and
$d\log\sigma/dq$ with which the distribution’s mean and (log) standard
deviation vary with $q$. They found $d\mu/dq=-0.46^{+0.29}_{-0.28}$ (median
and 90% credible uncertainties) and constrained $d\mu/dq<0$ at 99%
credibility. The data does not determine whether the width of the
$\chi_{\mathrm{eff}}$ distribution varies with BBH mass ratio, and remain
consistent with a $\chi_{\mathrm{eff}}$ distribution of width $\sigma\approx
0.1$ across all mass ratios.

The strong preference for negative $d\mu/dq$ may indicate that BBH mass ratios
and effective spins are anti-correlated, with unequal-mass events possessing
(on average) positive $\chi_{\mathrm{eff}}$ and equal-mass binaries possessing
small or vanishing $\chi_{\mathrm{eff}}$. The anticipated detection of several
tens of additional BBH detections in LIGO & Virgo’s O3b observing run,
however, will further clarify the relationship between $q$ and
$\chi_{\mathrm{eff}}$.

##  3 Phenomenological Considerations

Here we shall make the assumption that all BBH mergers come from AGN. We can
then determine what regions of parameter space—and what physical conditions—we
require to produce some kind of anti-correlation between ($q,\chi_{\rm eff}$).
In the future, a detailed parameter space study could determine the permitted
values for AGN channel parameters (and their degeneracies), as well as the
degeneracies with various possible mixing fractions for other channels;
however that is beyond the scope of this work. In particular, we should bear
in mind that the anti-correlation might involve only a fractional contribution
from the AGN channel. For example, AGN might dominate the production of
asymmetric mass mergers involving higher masses and spin but might be sub-
dominant at lower $\chi_{\rm eff}$ or $q\sim 1$ (Gayathri et al., 2021). In
this case, it may be the mixture between the AGN and other binary formation
channels that gives rise to the observed correlation between mass ratio and
spin.

Below, we discuss properties of binaries and their encounters. We define
‘retrograde’ (‘prograde’) binaries as those which rotate around their center
of mass such that $\vec{L}_{\rm b}\cdot\vec{L}_{\mathrm{disk}}<0$
($\vec{L}_{\rm b}\cdot L_{\mathrm{disk}}>0$)—see Figure 1. The hard-soft
boundary for a binary is when the binary orbital energy is approximately equal
to the average energy ($M_{3}\Delta v^{2}$) of encounters. Softening
(hardening) encounters tend to increase (decrease) binary separation. A binary
is ‘ionized’ when $KE+PE>0$ after an encounter.

Figure 1: Prograde versus Retrograde Binaries: Single black holes in an AGN
disk that are orbiting the SMBH with their angular momentum aligned to the
disk angular momentum $\vec{L}_{\rm d}$ will form binaries whose center of
mass will continue to orbit with the disk gas. However, the orbital angular
momentum of the binary around its own center of mass may be either prograde
(aligned with $\vec{L}_{\rm d}$) or retrograde (anti-aligned with
$\vec{L}_{\rm d}$). We use the terms ‘prograde’ and ‘retrograde’ binary,
respectively, for these two arrangements. Because single BH orbiting
retrograde with respect to the AGN disk gas should have lifetimes short
compared to the gas disk (Secunda et al., 2020a), we do not expect binaries
whose center of mass orbits in a retrograde sense with respect to the disk
gas, and we neglect such objects entirely in our analysis.

###  3.1 Symmetry breaking in AGN

Classical dynamical merger channels (e.g. globular clusters and nuclear star
clusters) generally assume spherical symmetry. AGN are distinctive in that
there are several sources of spherical symmetry breaking, some of which may
play a role in the possible anti-correlation discussed here. We outline some
of the basic concepts here and elaborate further below.

The AGN disk tends to break spherical mass symmetry in a galactic nucleus
since 1) more massive BH are preferentially captured by the disk (Fabj et al.,
2020), 2) star formation in AGN disks should come with a top-heavy IMF (e.g.
Levin & Beloborodov, 2003) and 3) a combination of mass segregation/star
formation in galactic nuclei may yield disky assemblies of massive BH (e.g.
Alexander et al., 2007).

The symmetry of dynamical encounters in AGN (between a binary and a singleton)
may be broken because: 1) spherical symmetry for encounters in the galactic
nucleus is broken by the preference for $\vec{L}_{\rm b}$ to form parallel or
anti-parallel to $\vec{L}_{\rm d}$, and disk capture of singleton orbiters
then biases encounters towards co-planar arrangements 2) there may be regions
in the inner disk where migration stalls, e.g. migration traps, so encounters
will typically arrive from the outer disk, breaking even the planar symmetry
of encounters (Bellovary et al., 2016; Secunda et al., 2020b), 3) BH orbiting
the SMBH on retrograde orbits may rapidly decay onto the SMBH, biasing orbits
to prograde (Secunda et al., 2020a), 4) Retrograde binaries have a different
hard-soft boundary than prograde binaries for encounters that are prograde, 5)
massive objects migrate faster in the disk so there is a mass bias in ‘catch
up’ encounters, 6) gas torques may preferentially harden/soften
retrograde/prograde binaries, depending on the details of gas flow around the
BBH (Baruteau et al., 2011; Li et al., 2021).

###  3.2 Phenomenology

McKernan et al. (2020a, b) performed Monte-Carlo simulations of large numbers
of BBH mergers in 1D models of AGN disks spanning a wide range of assumptions.
Fig. 2(a) shows ($q,\chi_{\rm eff}$) for BBH mergers with assumptions from R7
in (McKernan et al., 2020a) with no discernible correlation in the simulated
mergers. Here we revisit our assumptions to investigate what is required to
generate a ($q,\chi_{\rm eff}$) anti-correlation for the AGN channel. It helps
to guide our intuition to consider that an anti-correlation in ($q,\chi_{\rm
eff}$) can appear only if we suppress mergers in specific parts of
($q,\chi_{\rm eff}$) parameter space. These regions, in turn, correspond to
particular physical processes. We discuss three particular regions of
parameter space (regions 1,2 & 3 highlighted in Fig. 2(a)) below, where BBH
mergers must be suppressed to obtain a correlation in
$q$-$\chi_{\mathrm{eff}}$ similar to that observed by Callister et al. (2021).

Figure 2:  Left panel: Mass ratio ($q=M_{2}/M_{1}$) as a function of
$\chi_{\rm eff}$ for binary black hole mergers in the bulk disk from 100
realizations of the R7 model of McKernan et al. (2020a). Blue points
correspond to mergers in the bulk disk, red points correspond to mergers at a
trap. Horizontal structures correspond to artefacts of mass bins among early
mergers (e.g. mergers between $5.0M_{\odot}+10.0M_{\odot}$ BH before any
significant accretion has occurred). The results show no interesting
correlations. In order to produce the $q-\chi_{\rm eff}$ anti-correlation
observed in Callister et al. (2021), there must be physical mechanisms for
suppressing 3 types of mergers, highlighted in the plot and discussed in the
text: 1) negative $\chi_{\rm eff}$ mergers (purple region); 2) mergers with
$q\ga 0.5$ and $\chi_{\rm eff}\ga 0.3$ (orange region); and 3) mergers with
$q\la 0.3$ and $\chi_{\rm eff}\la 0.5$ (green region). Right panel: Mass ratio
($q=M_{2}/M_{1}$) as a function of $\chi_{\rm eff}$ for BBH mergers at the
trap (red points) and the remaining bulk disk (blue points) in 100
realizations of the previous model, modified by: 1) suppressing mergers of
binaries which form with $\vec{L}_{\rm b}$ anti-aligned with respect to
$\vec{L}_{\rm d}$ (retrograde binaries, see also Figure 1); 2) decreasing the
number of high mass BH relative to low mass BH, and modifying migration
torques with a luminosity-dependent prescription from Hankla et al. (2020);
and 3) segregating BH by mass such that more massive BH are more likely to be
disk-embedded, and have their spin axis aligned with the disk (while less
massive BH have lower spins that are initially isotropic). See text for
details.

First, BBH mergers with negative $\chi_{\rm eff}$ (purple region of Fig. 2(a))
are suppressed with respect to positive $\chi_{\rm eff}$ mergers, except
possibly around $q\sim 1$. In order to do this, we require binaries in AGN to
preferentially form either:1) biased towards both prograde spins _and_
prograde orbital angular momentum, or 2) biased towards both retrograde spins
_and_ retrograde orbital angular momentum. Here we shall simply assume that
‘retrograde’ binaries are ionized or softened at much higher rate than
‘prograde’ binaries due to symmetry breaking (see below). We do not include in
our model orbiters moving backwards through the disk gas (anti-parallel
orbital angular momentum to $\vec{L}_{\rm d}$), since these orbits decay on
short-timescales relative to the lifetime of an AGN disk (Secunda et al.,
2020a).

Second, BBH mergers with $q\ga 0.5$ but $\chi_{\rm eff}\ga 0.3$ (orange region
of Fig. 2 (a)) are suppressed. The simplest way of doing this in the AGN
channel is to suppress the number of massive BH mergers relative to less
massive BH mergers. This requires decreasing the relative number of more
massive BH in the disk (e.g., steepening the slope of the BH mass function),
or disrupting their mass-dependent migration, or both.

Third, mergers with $q\la 0.3$ and $\chi_{\rm eff}\la 0.3$ (green region in
Fig. 2(a)) are suppressed. The simplest way of doing this in the AGN channel
is to assume that spin orientations in encounters between more massive and
less massive BH are not isotropically distributed. Rather, the more massive BH
have spin orientations that prefer weak alignment with $L_{\rm disk}$, and
consequently will be weakly aligned with the orbital angular momentum of
binaries formed in the disk ($\vec{L}_{\rm b}$). This will tend to push $q\la
0.3$ mergers to higher positive $\chi_{\rm eff}$ values. We shall consider
each of these suppressive effects in turn below.

###  3.3 Softening or Ionizing retrograde binaries: suppressing $\chi_{\rm
eff}<0$

BBH mergers with negative $\chi_{\rm eff}$ imply that the projection of spins
from masses $M_{1},M_{2}$ onto the binary orbital angular momentum about its
center of mass ($\vec{L}_{\rm b}$) is net negative. i.e. the spin orientations
are significantly anti-aligned with $\vec{L}_{\rm b}$. There is a relative
deficit of mergers (in O3a) with $\chi_{\rm eff}<0$ (Abbott et al., 2021b). In
a gas-free dynamics merger channel, such as globular or nuclear clusters, an
isotropic distribution of spin orientations should be expected and dynamical
encounters should be spherically symmetric. So approximately half of the
mergers in clusters would be expected to have negative $\chi_{\rm eff}$. In
AGN disks, while binaries can form with $\vec{L}_{b}$ aligned or anti-aligned
with the disk orbital angular momentum ($\vec{L}_{\rm d}$), the fate of the
binary depends on the details and rates of gas or dynamical hardening. The
preferential plane of the AGN gas disk provides a mechanism for breaking the
symmetries of a gas-free ‘pure’ dynamics channel, particularly at locations in
disks with stalling orbits like migration traps.

Gas-hardening efficiency should decrease at modest binary separations. Binary
hardening can then stall and the details of dynamical encounters with tertiary
objects in AGN disks, including other binaries, becomes important (e.g. Leigh
et al., 2018; Samsing et al., 2020; Tagawa et al., 2020a). Hard, stalled
binaries in AGN migrate essentially as a single object with mass $M_{\rm
bin}$, and encounter other migrators; only a few such close encounters are
needed to merge the binary (Leigh et al., 2018). Here we point out that there
are a few places where symmetry of dynamical encounters may be broken which
may contribute to a deficit of mergers with $\chi_{\rm eff}<0$. We caution
that the details of the dynamical encounters will be important for our
conclusions; there have not yet been sufficient large scale numerical
experiments of encounters in the regime where orbital angular momentum, non-
negligible energy of encounter $\Delta E$, and orbital element perturbation
are considered together.

First, and most obviously, the hard-soft boundary for retrograde and prograde
binaries is different for encounters on a given ’side’ of the binary. The
relative velocity, $\Delta v$, for encounters with migrators arriving from the
outer disk is larger for retrograde binaries (see Fig. 3). This effect is most
obvious at migration traps in the inner disk where most encounters arrive from
the outer disk. Thus, the symmetry of encounters between retrograde and
prograde binaries may be broken at a migration trap where most encounters
arrive from the outer disk.

Second, since more massive embedded objects migrate faster in disks, more
massive objects ‘catch up’ with less massive objects at smaller semi-major
axis. More massive encounters with a binary will have a higher (lower)
encounter velocity for retrograde (prograde) binaries. Conversely, if a binary
‘catches up’ to less massive interior orbiters, the encounter velocity is
lower (higher) for retrograde (prograde) binary orbits around its center of
mass—see Figure 3. As a result, the symmetry breaking provided by the disk on
the direction of typical disk encounter could bias hardening encounters to
those between prograde binaries and other migrators. However, the degree of
asymmetry of this effect depend on the details of the encounter. For example,
if a migrating tertiary spends an equal amount of time on both sides of a
binary then symmetry of interactions is restored.

Third, the (complicated) details of gas torques may yield a difference in the
hardening efficiency of retrograde versus prograde binaries (e.g. Baruteau et
al., 2011; Derdzinski et al., 2021; Li et al., 2021). Baruteau et al. (2011)
find that retrograde binaries harden faster than prograde binaries; Li et al.
(2021) find that gas torques soften prograde binaries and only retrograde
binaries will be preserved. If correct, the only way of preserving the bias
against negative $\chi_{\rm eff}$ mergers in the AGN channel, would be if the
BH in the disk have a bias towards negative spins. Interestingly, Jermyn et
al. (2021) suggest that stars embedded in AGN disks are driven to strongly
negative spins, which might yield a population of BH forming in the disk with
negative spins. However, we should be cautious about gas hardening results—in
all cases, the hardening torques examined are operating on binaries near the
hard-soft boundary, and while gas hardening may produce binaries with small
semi-major axes, such binaries will likely stall before GW efficient merger.
For stalled binaries, dynamical interactions must become involved to produce
mergers, but could also ionize the binaries (see below). In addition, while Li
et al. (2021) do not produce hardened prograde binaries, they emphasize the
critical effect of the detailed distribution of gas nearest to each of the
binary partners in determining torque sign. Given given the effect of going
from 2D to 3D simulations (see e.g. Zhu et al., 2012; Disk Dynamics
Collaboration et al., 2020, for work in the protoplanetary context), and the
lack of feedback in any binary hardening simulations to date (which will
certainly change the distribution of gas near each BH), we should be cautious
in taking present conclusions at face value.

One might also consider retrograde migrators interacting with prograde
binaries, but the gas disk provides a very strong symmetry breaking effect; as
shown in Secunda et al. (2020a), single retrograde migrators will rapidly
decay into the SMBH due to eccentricity pumping and will be unavailable to
interact with any other disk components after a relatively short period
$O(10^{5})~{}$yr.

###  3.4 Suppressing $q\ga 0.5$, $\chi_{\rm eff}\ga 0.3$ mergers

In general, near equal mass, high $\chi_{\rm eff}$ mergers in the AGN channel
come from the merger of hierarchical merger products, often at or near a
migration trap. In order to suppress these mergers, we must decrease the
relative number of high mass, high spin BH embedded in the disk (compared to
the spherical distribution), or prevent these objects from efficiently finding
each other, as generally happens if the migration trap is universal. To reduce
the number of of high-mass disk-embedded objects, we can give the disk
embedded objects a flatter mass function than those in the spheroid
population; to reduce the efficiency of trap mergers we can disrupt the smooth
migration torques provided by the disk gas. Migration disruption can occur if
AGN disks are typically highly turbulent or if feedback from accretion onto
the embedded population alters the migration torques. For example, Hankla et
al. (2020) showed that accretion onto orbiters embedded in a disk generates a
heating torque that can act to drive the orbiters outwards, i.e., against the
typical sense of Type I migration within disks; importantly, the strength of
the counter-directional heating torque is related to the mass of the migrator
through its accretion luminosity, which scales with mass.

###  3.5 Spin alignment among heavier BH: suppressing $q\la 0.3$, $\chi_{\rm
eff}\la 0.3$

Since the gas disk provides a preferential plane for binary formation,
aligning the spin of the heavier BH with the orbital angular momentum of the
binary ($\vec{L}_{\rm b}$) can be physically motivated by any mechanisms that
preferentially align the spin of the primary with the orbital angular momentum
of the disk ($\vec{L}_{\rm d}$, and hence the orbital angular momentum of the
binary, $\vec{L}_{\rm b}$). A few percent ($\sim 1-10\%$) mass accretion is
sufficient to torque the spin of a mis-aligned BH into alignment with a disk
(Bogdanović et al., 2007). So any mechanism that leads to an angular mass
segregation, such that more massive BH spend more time in the gas disk (or
form in the disk, or arrive there substantially before the less massive BH)
will provide appropriate conditions for the embedded BH population (i.e.
$\vec{\chi}_{1}$ is parallel to $\vec{L}_{\rm b}$).

Fortunately, if we begin with a spheroidal distribution of individual BH which
has no preferred orbital plane and relatively small natal spins that are
randomly aligned, the introduction of a gas disk will provide a torque that
will align the orbits of BH with the gas disk. The dominant aligning torque is
Bondi-Hoyle-Lyttleton drag, which depends on $M^{2}$, resulting in faster
alignment of more massive BH orbits with the AGN disk (Fabj et al., 2020).
Such BH have more time to accrete, and thus they have more time to have their
(presumably random initial) spins aligned with the disk compared to their less
massive partners. This mechanism requires the inner AGN disk to be dense
enough that orbital capture is efficient within AGN lifetimes (Fabj et al.,
2020; MacLeod & Lin, 2020). Another possibility is that there is a
preferential plane for AGN activity in galactic nuclei, such that angular mass
segregation occurs due to leftover mergers from prior AGN episodes in the same
plane (e.g. Berti & Volonteri, 2008; Fanidakis et al., 2011); or that mutual
interactions even in the absence of a gas disk provide a preferential plane
that results in disky assemblages of heavier BH at all times (e.g. Alexander
et al., 2007), or that compact objects in the disk form from a top heavy IMF
(e.g. Levin & Beloborodov, 2003), or are driven to higher masses through
accretion (Cantiello et al., 2021). Such mechanisms would enhance angular mass
segregation and drive the preferential spin alignment of the heavier partners
in BBH mergers.

##  4 Modifications to Monte Carlo simulations

Based on the phenomenological considerations above, we made some very simple
modifications to the Monte Carlo simulations of McKernan et al. (2020b) as
described below. We emphasize again that these initial choices of parameters,
while physically motivated, are only to illustrate that a population with a
($q,\chi_{\rm eff}$) anti-correlation can be produced by an AGN channel.
Detailed parameter ranges will be studied in future work.

For initial mass distributions, we introduced 1) a disk mass function
($M^{-1}$) to correspond to a heavier disk BH population ([$20,50M_{\odot}$]).
We also introduced 2) a spheroid BH mass function ($M^{-2}$) from which we
randomly draw the lighter captured BH population ($[5,30]M_{\odot}$). We
normalized both distributions so that $\sim 40$ BH from 1) start in the disk
and $\sim 200$ BH from 2) are captured by the disk in $1\rm{Myr}$, for each
realization of the AGN disk.

For initial spin distributions, we assumed that the population in 2) had
random spin orientation but a narrow range of dimensionless spin magnitudes
parameterized by $a=s(1-s)$ where $s=[-0.98,0.98]$ is a random draw from a
uniform distribution. This is not a physical model, but limits initial spin
magnitudes to a draw between $a=[-0.25,0.25]$ among the captured BH population
and corresponds to an assumption of small, but non-zero initial BH spin
magnitude at birth. For the population in 1), we assumed that by starting off
embedded within the AGN disk, BH had been slightly spun-up at $t=0$. We
parameterized this by drawing initial spin magnitudes from a distribution
given by $s(1-s)+0.15$ where $s=[-0.98,0.98]$. Again, this is not a physical
model, but has the effect of biasing initial BH dimensionless spin parameters
in the disk to a draw between $a=[-0.1,0.35]$, corresponding to an assumption
that BH have been marginally spun up, either from direct accretion from the
disk, or possibly from previous AGN episodes in the same plane. Our average
spin bias of $s(1-s)+0.15$ selected here corresponds to around a Myr of gas
accretion at the Eddington rate, or equivalently, shorter times at super-
Eddington rates. The corresponding spin alignments in population 2) are
random. A detailed study of the allowed range of spin-up by disks will be
carried out in future work.

For accretion over the course of the simulation, we assumed that BH that are
initially embedded in the disk accrete at slightly super-Eddington ($x2$)
rates. This has the effect of driving a faster torquing of spin orientation
into alignment with the AGN disk (Bogdanović et al., 2007) over 1Myr, with a
modest change in mass.

For migration, we assumed that feedback from accretion modifies the gas
torques on embedded objects. We used the parameterization from (Hankla et al.,
2020), which approximates the ratio of the heating torque ($\Gamma_{\rm
heat}$) to the migration torque ($\Gamma_{\rm mig}$) is

$\frac{\Gamma_{\rm heat}}{\Gamma_{\rm mig}}\approx 0.07\left(\frac{c}{v_{\rm
orb}}\right)\epsilon\tau^{-1}\alpha^{-3/2}$ (2)

where $v_{\rm orb}$ is the orbital velocity of the migrator, $\epsilon$ is the
Eddington ratio of accretion, $\tau$ is the disk optical depth and $\alpha$ is
the disk viscosity parameter. The result is a net outward heating torque on
embedded objects, which gets co-added to the (generally) inward migration
torque and has the effect of slowing down the general drift inwards of
population 1) BH. The heating torque depends strongly on the disk optical
depth, so a non-uniform quality now appears in the rate of migration of
embedded objects, depending on their radial location in the disk.

For disk size, we assumed that the dense disk cuts off at $10^{4}r_{g}$ (or
$0.05~{}$pc for $M_{SMBH}=10^{8}~{}M_{\odot}$). This has the effect of
inhibiting the number of embedded objects from population 1) and therefore
also reduces the relative number of mergers of within population 1). This size
scale is larger than the O($10^{3}r_{g}$)inferred from reverberation mapping
or microlensing (Chartas et al., 2016), so could be revised in a parameter
study, but is roughly where star formation is expected to consume the unstable
outer disk (e.g. Sirko & Goodman, 2003).

For dynamical interactions, we simply assumed that all binaries which would
have formed retrograde binaries, fail to form, and the individual components
continue their migration undisturbed. We assumed that there is a migration
trap in this disk model. The efficiency of the migration trap depends on the
surface density change in the inner disk (Bellovary et al., 2016) as well as
the effect of a Bondi accretion headwind (Pan & Yang, 2021). The latter
depends on the details of gas flow around embedded objects, which remains
poorly understood. If an AGN disk exhbits large changes in surface density,
particularly in the inner disk, where the pressure is high and disk mass is
small, we should expect regions of the disk where in-migration can slow or
stall and a pile-up can occur. Indeed in McKernan et al. (2020b) we found that
when we excluded the migration trap in a Thompson et al. (2005) disk model,
following the treatment of Dittmann & Miller (2020), the change in disk
surface density remained sufficient to substantially decrease migration
torques. This allowed a large number of interactions and mergers in what we
termed a migration ‘swamp’. Importantly, migration stalling strengthens the
symmetry-breaking, binary ionization dynamic shown in Fig. 3, since most
encounters for a stalled binary are due to ‘catch-up’ from larger disk radii.

Fig. 2b shows the results of implementing all of the above modifications in
our simulations. The points in Fig. 2b correspond to all mergers in $100$
different realizations of a migration trap (red points) and bulk disk (blue
points) in a Sirko & Goodman (2003) model AGN disk around a $10^{8}M_{\odot}$
SMBH. The migration trap is located at $700r_{g}$ and the disk has a lifetime
of 1Myr. This most nearly matches run R7 in McKernan et al. (2020b), though
the number of BH, their capture rate, their mass functions, and the existence
of retrograde binaries have all been altered as discussed above. A clear trend
in ($q,\chi_{\rm eff}$) is apparent for the trap mergers (red points),
although the overlap with the 1$\sigma$ distribution from (Callister et al.,
2021) is not particularly good. There is significant scatter about the trend,
particularly at $q\ga 0.5$, but smaller scatter in $\chi_{\rm eff}$ at $q\la
0.5$ in trap mergers. We emphasize that our goal here is to establish what
sort of phenomenological conditions are required in order to generate a
possible ($q,\chi_{\rm eff}$) anti-correlation in Monte Carlo simulations of
BBH mergers in AGN disks. While our initial assumptions here yield plausible
results relatively similar to (Callister et al., 2021), a N-body study is
required to test ionization dynamics and a detailed parameter space MC study
is required to establish what conditions in AGN disks are ruled out by
Callister et al. (2021). In particular, since the trap mergers yield a
stronger anti-correlation, it may be that we need to reduce the size of the
bulk disk to $<10^{4}r_{g}$ to minimize the bulk disk contribution.
Nevertheless, at a qualitative level, our results strongly suggest that while
gas torques may play an important role in forming BBH in AGN disks, dynamical
interactions are necessary to reach a stage of GW-driven merger. This, in
turn, suggests a modest stalling radius for gas-driven binary hardening
torques, and the common existence of migration traps in AGN disks.

Figure 3:  Breaking the symmetry of hardening/softening encounters: The hard-
soft boundary is different for pro- and retro-grade binaries for encounters at
larger and smaller disk radii. This symmetry breaking effect is strongest at
migration traps (or other stalling orbits), since encounters with binaries at
a trap will mostly be with migrators from larger disk radii.

##  5 Discussion

The possibility of an anti-correlation in ($q,\chi_{\rm eff}$) parameter space
of BBH mergers is exciting, since it will require distinct constraints on all
of the proposed merger channels (Callister et al., 2021). Here we point out
that the AGN channel has multiple conditions that may contribute to symmetry
breaking in spherical mass distribution, binary formation and binary hardening
in the disk.

The key details underpinning a phenomenological anti-correlation between
$(q,\chi_{\rm eff})$ in the AGN channel are symmetry-breaking in both the
angular mass distribution of BH and in the hardening of binaries to merger
(either dynamical or gas-driven). An AGN disk can be thought of as breaking
the spherical symmetry of a gas-free galactic nucleus and promoting angularly
dependent mass segregation. More massive BH are more quickly captured by the
AGN disk (Fabj et al., 2020) and will therefore spend more overall time in the
disk, spinning up and experiencing torquing towards alignment with the disk.
However, without the assumption of an offset for the spin magnitudes of the
initial, disk-embedded, heavy BH population, we are unable to obtain such a
clear anti-correlation among trap mergers. Thus, a separate possibility is
that mass segregation in gas-free galactic nuclei leads to disky distributions
of massive BH. Such disky distributions would provide a preferential plane for
AGN disk formation, exchanging angular momentum with incoming gas. Or, a long-
term preferential plane could be established if there are multiple AGN
episodes delivered from the same fuel reservoir. Likewise, a disky
distribution of BH implies the interval between fuelling episodes (duty cycle)
is less than the relaxation timescale. The fuel reservoir (torus) can be
volume filling and must therefore contain a fraction of the nuclear star
cluster. Objects embedded in the torus will accrete, spin up and torque over
time and would naturally end up embedded in inflowing gas. Yet another
possibility is a top-heavy IMF for star formation in AGN disks which naturally
promotes a more massive disky component (Levin & Beloborodov, 2003). A
combination of these mechanisms could account for a more massive disk
component.

Fig. 4 is a cartoon summarizing (in blue) what we can naturally generate in
AGN models, and (in red), _predictions_ for BBH mergers with negative
$\chi_{\rm eff}$ from this channel. If retrograde binary ionization were
irrelevant to BBH mergers in AGN disks then we should broadly expect a mirror
image distribution of Fig. 2(b) at $\chi_{\rm eff}<0$ (i.e. mapping through
$\chi_{\rm eff}=0$). However, if retrograde binary ionization is important for
this channel, then there are two clear expectations: First, the absolute
number of mergers with $\chi_{\rm eff}<0$ must be less than the corresponding
number of mergers with $\chi_{\rm eff}>0$. Second, those retrograde binaries
that survive to merger may be massive enough compared to dynamical encounters
that they can survive. Thus, we should expect the slope of the ($q,\chi_{\rm
eff}$) distribution in AGN to be shallower at $\chi_{\rm eff}<0$ than at
$\chi_{\rm eff}>0$.

We show that an anti-correlation in ($q,\chi_{\rm eff}$) parameter space can
arise naturally in the AGN channel provided that a few phenomenological
assumptions hold: 1) Heavier BH live in the AGN disk and tend to spin up into
alignment with the disk. 2) Lighter BH live in the spheroid surrounding the
AGN disk and have random spin alignments. 3) The inner AGN disk is dense, in
order that BH from 2) are captured, but 4) not large radially, to limit the
number of BH from population 1). 5) Migration must not be smooth, either due
to a turbulent disk and/or feedback from the embedded objects, to cap the
number of mergers between BH in 1). 6) Dynamical encounters within the disk
must be common enough that binaries which orbit retrograde about their center
of mass (anti-aligned with the disk angular momentum) are preferentially
ionized. The ionization symmetry-breaking effect is strongest at stalled
migration orbits, e.g. migration traps. Our understanding is that at present,
no other BBH merger model channel can naturally produce such an anti-
correlation. In particular, isolated binary evolution tends to lead to BH
spins which are both aligned with the orbital angular momentum of the binary,
and also tends towards $q\sim 1$ mergers; similarly, globular cluster dynamics
tends to sort the masses in binaries towards $q\sim 1$ (through exchange
interactions), and robustly produces symmetric distributions of $\chi_{\rm
eff}$.

Figure 4: Cartoon representation of populations expected from this model. In
blue are the $\chi_{\rm eff}>0$ populations as discussed in the text. In red
are predictions for a (small) population of mergers due to formation of some
BBH with $\chi_{\rm eff}<0$. If the hard-soft boundary symmetry breaking for
retro- and pro-grade binaries is important, then we expect the slope of the
$\chi_{\rm eff}<0$ merger population to be shallower than that for the
$\chi_{\rm eff}>0$ population (see text).

## Acknowledgements.

We would like to thank Imre Bartos for helpful feedback and comments. BM &
KESF are supported by NSF AST-1831415 and Simons Foundation Grant 533845.

## Data Availability

Any data used in this analysis are available on reasonable request from the
first author (BM).

## References

  * Aasi et al. (2015) Aasi J., Abbott B. P., Abbott R., Abbott T., Abernathy M. R., et al., 2015, Class. Quant. Grav., 32, 074001
  * Abbott et al. (2020a) Abbott R., et al., 2020a, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2020b) Abbott R., et al., 2020b, ApJ, 896, L44
  * Abbott et al. (2020c) Abbott R., et al., 2020c, ApJ, 900, L13
  * Abbott et al. (2021a) Abbott R., et al., 2021a, Phys. Rev. X, 11, 021053
  * Abbott et al. (2021b) Abbott R., et al., 2021b, ApJ, 913, L7
  * Acernese et al. (2015) Acernese F., Agathos M., Agatsuma K., Aisa D., Allemandou N., et al., 2015, Class. Quant. Grav., 32, 024001
  * Alexander et al. (2007) Alexander R. D., Begelman M. C., Armitage P. J., 2007, ApJ, 654, 907
  * Ashton et al. (2020) Ashton G., Ackley K., Hernandez I. M., Piotrzkowski B., 2020
  * Bartos et al. (2017) Bartos I., Kocsis B., Haiman Z., Márka S., 2017, ApJ, 835, 165
  * Baruteau et al. (2011) Baruteau C., Cuadra J., Lin D. N. C., 2011, ApJ, 726, 28
  * Bellovary et al. (2016) Bellovary J. M., Mac Low M.-M., McKernan B., Ford K. E. S., 2016, ApJ, 819, L17
  * Berti & Volonteri (2008) Berti E., Volonteri M., 2008, ApJ, 684, 822
  * Bogdanović et al. (2007) Bogdanović T., Reynolds C. S., Miller M. C., 2007, ApJ, 661, L147
  * Callister et al. (2021) Callister T. A., Haster C.-J., Ng K. K. Y., Vitale S., Farr W. M., 2021, arXiv e-prints, p. arXiv:2106.00521
  * Cantiello et al. (2021) Cantiello M., Jermyn A. S., Lin D. N. C., 2021, ApJ, 910, 94
  * Chartas et al. (2016) Chartas G., et al., 2016, Astronomische Nachrichten, 337, 356
  * Derdzinski et al. (2021) Derdzinski A., D’Orazio D., Duffell P., Haiman Z., MacFadyen A., 2021, MNRAS, 501, 3540
  * Disk Dynamics Collaboration et al. (2020) Disk Dynamics Collaboration et al., 2020, arXiv e-prints, p. arXiv:2009.04345
  * Dittmann & Miller (2020) Dittmann A. J., Miller M. C., 2020, MNRAS, 493, 3732
  * Fabj et al. (2020) Fabj G., Nasim S. S., Caban F., Ford K. E. S., McKernan B., Bellovary J. M., 2020, MNRAS, 499, 2608
  * Fanidakis et al. (2011) Fanidakis N., Baugh C. M., Benson A. J., Bower R. G., Cole S., Done C., Frenk C. S., 2011, MNRAS, 410, 53
  * Gayathri et al. (2021) Gayathri V., Yang Y., Tagawa H., Haiman Z., Bartos I., 2021, arXiv e-prints, p. arXiv:2104.10253
  * Generozov et al. (2018) Generozov A., Stone N. C., Metzger B. D., Ostriker J. P., 2018, MNRAS, 478, 4030
  * Gerosa & Berti (2019) Gerosa D., Berti E., 2019, Phys. Rev. D, 100, 041301
  * Graham et al. (2020) Graham M. J., et al., 2020, Phys. Rev. Lett., 124, 251102
  * Gröbner et al. (2020) Gröbner M., Ishibashi W., Tiwari S., Haney M., Jetzer P., 2020, A&A, 638, A119
  * Hailey et al. (2018) Hailey C. J., Mori K., Bauer F. E., Berkowitz M. E., Hong J., Hord B. J., 2018, Nature, 556, 70
  * Hankla et al. (2020) Hankla A. M., Jiang Y.-F., Armitage P. J., 2020, ApJ, 902, 50
  * Jermyn et al. (2021) Jermyn A. S., Dittmann A. J., Cantiello M., Perna R., 2021, arXiv e-prints, p. arXiv:2102.13114
  * Leigh et al. (2018) Leigh N. W. C., et al., 2018, MNRAS, 474, 5672
  * Levin & Beloborodov (2003) Levin Y., Beloborodov A. M., 2003, ApJ, 590, L33
  * Li et al. (2021) Li Y.-P., Dempsey A. M., Li S., Li H., Li J., 2021, ApJ, 911, 124
  * MacLeod & Lin (2020) MacLeod M., Lin D. N. C., 2020, ApJ, 889, 94
  * McKernan et al. (2012) McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MNRAS, 425, 460
  * McKernan et al. (2014) McKernan B., Ford K. E. S., Kocsis B., Lyra W., Winter L. M., 2014, MNRAS, 441, 900
  * McKernan et al. (2018) McKernan B., Ford K. E. S., Bellovary J., et al., 2018, ApJ, 866, 66
  * McKernan et al. (2019) McKernan B., Ford K. E. S., Bartos I., et al., 2019, ApJ, 884, L50
  * McKernan et al. (2020a) McKernan B., Ford K. E. S., O’Shaugnessy R., Wysocki D., 2020a, MNRAS, 494, 1203
  * McKernan et al. (2020b) McKernan B., Ford K. E. S., O’Shaughnessy R., 2020b, MNRAS, 498, 4088
  * Miralda-Escudé & Gould (2000) Miralda-Escudé J., Gould A., 2000, ApJ, 545, 847
  * Morris (1993) Morris M., 1993, ApJ, 408, 496
  * Nitz et al. (2021) Nitz A. H., Capano C. D., Kumar S., Wang Y.-F., Kastha S., Schäfer M., Dhurkunde R., Cabero M., 2021, arXiv e-prints, p. arXiv:2105.09151
  * Palmese et al. (2021) Palmese A., Fishbach M., Burke C. J., Annis J. T., Liu X., 2021, arXiv e-prints, p. arXiv:2103.16069
  * Pan & Yang (2021) Pan Z., Yang H., 2021, arXiv e-prints, p. arXiv:2101.09146
  * Samsing et al. (2020) Samsing J., et al., 2020, arXiv e-prints, p. arXiv:2010.09765
  * Schawinski et al. (2015) Schawinski K., Koss M., Berney S., Sartori L. F., 2015, MNRAS, 451, 2517
  * Secunda et al. (2019) Secunda A., Bellovary J., Mac Low M.-M., et al., 2019, ApJ, 878, 85
  * Secunda et al. (2020a) Secunda A., Hernandez B., Goodman J., Leigh N. W. C., McKernan B., Ford K. E. S., Adorno J. I., 2020a, arXiv e-prints, p. arXiv:2009.03910
  * Secunda et al. (2020b) Secunda A., et al., 2020b, ApJ, 903, 133
  * Sirko & Goodman (2003) Sirko E., Goodman J., 2003, MNRAS, 341, 501
  * Stone et al. (2017) Stone N. C., Metzger B. D., Haiman Z., 2017, MNRAS, 464, 946
  * Tagawa et al. (2019) Tagawa H., Haiman Z., Kocsis B., 2019, arXiv e-prints, p. arXiv:1912.08218
  * Tagawa et al. (2020a) Tagawa H., Haiman Z., Kocsis B., 2020a, ApJ, 898, 25
  * Tagawa et al. (2020b) Tagawa H., Haiman Z., Bartos I., Kocsis B., 2020b, ApJ, 899, 26
  * Tagawa et al. (2021) Tagawa H., Haiman Z., Bartos I., Kocsis B., Omukai K., 2021, arXiv e-prints, p. arXiv:2104.09510
  * Thompson et al. (2005) Thompson T. A., Quataert E., Murray N., 2005, ApJ, 630, 167
  * Venumadhav et al. (2020) Venumadhav T., Zackay B., Roulet J., Dai L., Zaldarriaga M., 2020, Phys. Rev. D, 101, 083030
  * Yang et al. (2019) Yang Y., Bartos I., Gayathri V., et al., 2019, Phys. Rev. Lett., 123, 181101
  * Yang et al. (2020) Yang Y., Gayathri V., Bartos I., Haiman Z., Safarzadeh M., Tagawa H., 2020, ApJ, 901, L34
  * Zhu et al. (2012) Zhu Z., Stone J. M., Rafikov R. R., 2012, ApJ, 758, L42