# Ferromagnetic domains in the large-$U$ Hubbard model with a few holes: an
FCIQMC study

Sujun Yun<EMAIL_ADDRESS>Max Planck Institute for Solid State Research,
Heisenbergstr. 1, 70569 Stuttgart, Germany  School of electronic engineering,
Nanjing XiaoZhuang University, Hongjing Road, Nanjing 211171, China  Werner
Dobrautz<EMAIL_ADDRESS>Department of Chemistry and Chemical
Engineering, Chalmers University of Technology, 41296 Gothenburg, Sweden
Hongjun Luo<EMAIL_ADDRESS>Max Planck Institute for Solid State Research,
Heisenbergstr. 1, 70569 Stuttgart, Germany  Vamshi Katukuri  Max Planck
Institute for Solid State Research, Heisenbergstr. 1, 70569 Stuttgart, Germany
Niklas Liebermann  Max Planck Institute for Solid State Research,
Heisenbergstr. 1, 70569 Stuttgart, Germany  Ali Alavi<EMAIL_ADDRESS>Max
Planck Institute for Solid State Research, Heisenbergstr. 1, 70569 Stuttgart,
Germany  Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom

###### Abstract

Two-dimensional Hubbard lattices with two or three holes are investigated as a
function of $U$ in the large-$U$ limit. In the so-called Nagaoka limit (one-
hole system at infinite $U$), it is known that the Hubbard model exhibits a
ferromagnetic ground-state. Here, by means of exact FCIQMC simulations applied
to periodic lattices up to 24 sites, we compute spin-spin correlation
functions as a function of increasing $U$. The correlation functions clearly
demonstrate the onset of ferromagnetic domains, centred on individual holes.
The overall total spin of the wavefunctions remain the lowest possible (0, or
$\frac{1}{2}$, depending on the number of holes). The ferromagnetic domains
appear at interaction strengths comparable to the critical interaction
strengths of the Nagaoka transition in finite systems with strictly one hole.
The existence of such ferromagnetic domains is the signature of Nagaoka
physics in Hubbard systems with a small (but greater than one) number of
holes.

##  I Introduction

The Hubbard model is a simple yet important model in the study of correlated
electrons, as it captures complex correlations between electrons on a lattice
with a fairly simple Hamiltonian [1]. Exact results of two-dimensional Hubbard
model are helpful for understanding a plethora of phenomena in strongly
correlated systems, including pairing mechanisms in unconventional
superconductors [2], the Mott metal-insulator transition [3], optical
conductivity [4, 5], and itinerant magnetism [6, 7]. For the single-band 2D
Hubbard model on a square lattice, Nagaoka [8] proved analytically that the
limit of infinitely strong interactions, in the presence of a single hole on
top of a Mott-insulating state with one electron per site, results in a
ferromagnetic ground state. Intuitively, the Nagaoka ferromagnetism can be
understood as resulting from an interference effect between the different
paths that the hole can take across the lattice. When the spins are aligned,
these paths interfere constructively and lower the kinetic energy of the hole
[9, 10, 11, 12].

While Nagaoka ferromagnetism has been analytically proven under extreme
conditions, and has also been observed in a quantum dot plaquette [13], the
stability of the Nagaoka ferromagnetic state at finite interaction strengths
on finite lattices has also been actively studied [14, 15, 16, 17, 18, 19, 20,
21]. However, open questions still exist, especially concerning the
thermodynamic stability of the ferromagnetic state for systems with more than
one hole. Extrapolations from the results on finite lattices have been used to
study properties in the thermodynamic limit. Thus it is important to obtain
exact results in systems with two-, three-, and perhaps more-hole systems, on
lattices as large as possible.

In the two-hole system, the total spin of the ground-state is zero ($S=0$).
This has been numerically studied by exact diagonalization (ED) [22], the
spin-adapted full configuration interaction quantum Monte Carlo (FCIQMC)
method [21] and analytical studies of arbitrarily large systems [19]. However,
the specific type of magnetism is unknown, because anti-ferromagnetism, para-
magnetism, as well as low-spin-coupled ferromagnetic domains, all correspond
to $S=0$. For the three-hole system, on the other hand, ED results of an
effective Hamiltonian [22] show that the total spin of ground-states on the 8-
and 16-site lattices are $\frac{3}{2}$ and $\frac{7}{2}$ respectively. Here
the effective Hamiltonian is constructed to exclude double occupation and thus
should be related to the Hubbard model in the large $U$ limit. It is
interesting to see whether the partial magnetization will still be observed on
larger lattices. In order to answer the above open questions, in this work we
investigate the two- and three-hole systems with the full configuration
interaction quantum Monte Carlo (FCIQMC) method.

FCIQMC is based on stochastic simulations of the dynamic evolution of the wave
function in imaginary time. Different to the diffusion quantum Monte Carlo
(DMC) [23] with the fixed node approximation and the auxiliary-field quantum
Monte Carlo (AFQMC) [24] with the phaseless approximation, no systematic
approximation is made in FCIQMC [25, 26, 27], and it thus serves as a highly
accurate method to approach the ground state wave functions. The annihilation
procedure of the algorithm enables it to overcome the fermionic sign problem
exactly, as long as enough walkers are used. In practice, for Hubbard systems
at large $U$, this means, with the currently available hardware, systems up to
26 sites can be studied [21]. (A 26-site lattice represents a useful increase
in size compared to exact-diagonalisation, for which 20-sites is the largest
lattice size so far reported [4]). In the present paper, we extend this to the
study of a system with a few holes, as well as report on spin-spin correlation
functions for the obtained exact ground-states.

The software NECI, a state-of-the-art implementation of the FCIQMC algorithm,
utilizes a very powerful parallelization and scales efficiently to more than
24000 central processing unit cores [27]. The FCIQMC method in a Slater
determinant (SD) basis has been extended to calculate ground and excited state
energies, spectral and Green’s functions for ab initio and model systems, as
well as properties via the one-, two-, three- and four-body reduced density
matrices(RDMs). To study magnetism, we need to use the replica-sampled 2-RDMs
[28, 29, 30] to obtain the spatial spin distribution. The replica-sampling
technique removes the systematic error in the RDM, at the expense of requiring
a second walker distribution. The premise is to ensure that these two walker
distributions are entirely independent and propagated in parallel, sampling
the same (in this instance ground-state) distribution. This ensures an
unbiased sampling of the desired RDM, by ensuring that each RDM contribution
is derived from the product of an uncorrelated amplitude from each replica
walker distribution. By using replica-sampled 2-RDMs the spin-spin correlation
function, $\langle\hat{\bf S}_{i}\cdot\hat{\bf S}_{j}\rangle$, can be
calculated, where $i$ and $j$ are lattice site indices. This spin-spin
correlation function can then be used to identify the specific type of
magnetism of the ground states.

FCIQMC in a spin-adapted basis is also used to study the partial polarization
in three-hole systems. Spin-adapted FCIQMC uses $SU(2)$ symmetry (arising from
the vanishing commutator $[\hat{H},\hat{{\bf S}}^{2}]=0$) conservation.
$SU(2)$ symmetry is imposed via the graphical unitary group approach (GUGA)
[31, 32, 33] which dynamically constrains the total spin $S$ of a multi-
configuration and highly open-shell wave function in an efficacious manner.
The spin-adapted version of the FCIQMC algorithm based on GUGA has been
developed in our group [34, 35], with – among others – applications to _ab
initio_ system [36, 37] and Nagaoka ferromagnetism in one-hole system [21].
With the spin-adapted method, the magnetisation of the ground state can be
determined in a reliable way, especially for systems with small spin gaps. The
results of spin-adapted FCIQMC show the partial spin-polarization only appears
in small, three-hole system (less than 18 sites) [22], which is the second
important result of this work.

The rest of this paper is organized as follows: In Sec. II, we briefly
describe the methods, where we mainly provide some more details on the
measurements of the spin-spin correlation function, $\langle\hat{\bf
S}_{i}\cdot\hat{\bf S}_{j}\rangle$ from the replica-sampled 2-RDMs in FCIQMC.
In Sec. III, results about the spatial spin distribution and partial spin-
polarization are discussed. Finally, we conclude in Sec. IV.

##  II methods

The Hamiltonian of the Hubbard model in real space takes the form

$\hat{H}=-t\sum_{\langle
ij\rangle\sigma}a^{\dagger}_{i,\sigma}a_{j,\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}$
(1)

where $a^{\dagger}_{i\sigma}$ ($a_{i\sigma}$) creates (annihilates) an
electron with spin $\sigma$ on site $i$, and
$n_{i\sigma}=a^{\dagger}_{i\sigma}a_{i\sigma}$ is the particle number
operator. $U$ refers to the Coulomb interaction strength. We consider only
nearest neighbour hopping terms, where $t$ is positive and is used as the unit
of the energy. When $U$ is infinitely large, there will be no double occupancy
and the system can be treated with an effective Hamiltonian with constrained
hopping terms [22]

$H_{\textnormal{eff}}=-t\sum_{\langle
ij\rangle\sigma}\tilde{a}^{\dagger}_{i,\sigma}\tilde{a}_{j,\sigma},$ (2)

with $\tilde{a}^{\dagger}_{i,\sigma}=a_{i,\sigma}(1-n_{i,\sigma})$. In our
current work, we want to study the magnetic properties for finite $U$ and thus
will stay with the original Hamiltonian (1). Tough, we find that our results
for the three-hole systems in the large $U$ limit (see Sec. III(B)), coincide
with the result of Riera _et al._ [22] for the effective Hamiltonian, Eq. (2).
In our investigation we apply two different FCIQMC methods, which are based on
full CI expansions in terms of SDs and in terms of spin eigenfunctions (spin-
adapted basis states) respectively.

FCIQMC is a projector QMC method for obtaining the ground state wave function
$\ket{\Psi_{0}}$. By Monte Carlo simulation of the imaginary-time evolution of
the wave function

$\ket{\Psi(\tau)}=e^{-\tau(\hat{H}-E_{0})}\ket{\Psi(0)},$ (3)

the ground state wave function is approached in the long time limit
$\ket{\Psi(\tau\rightarrow\infty)}\propto\ket{\Psi_{0}}$.

In a previous work [21], we have investigated the magnetism for one hole and
two holes systems by using the spin-adapted ($SU(2)$ conserving) FCIQMC
method. We extend these investigations to three-hole systems in this work.
With the spin-adapted method, the magnetisation of the ground state can be
determined in a reliable way, especially for systems with small spin gaps. For
details of the spin-adapted FCIQMC method, we refer to previous work [21].

We also calculate the spin-spin correlation function to get a knowledge of the
spatial spin distribution. For these calculations, we find that the FCIQMC
method based on SDs is more efficient, in particular when it is combined with
the replica-sampling techniques for the 2-RDMs. With this technique, two
independent FCIQMC simulations are performed in parallel and the wave function
expansions

$\ket{\Psi(\tau)^{I/II}}=\sum_{\mu}c_{\mu}^{I/II}(\tau)|D_{\mu}\rangle,$ (4)

are sampled simultaneously for these two replicas. The sampled coefficients
$\\{c_{\mu}^{I}\\}$ and $\\{c_{\mu}^{II}\\}$ are then used to calculate the
2-RDM elements

$\displaystyle\Gamma_{pq}^{rs}$ $\displaystyle=$
$\displaystyle\langle\Psi|a^{\dagger}_{p}a^{\dagger}_{q}a_{s}a_{r}|\Psi\rangle$
(5) $\displaystyle=$
$\displaystyle\sum_{\mu,\nu}c^{I}_{\mu}c^{II}_{\nu}\langle
D_{\mu}|a^{\dagger}_{p}a^{\dagger}_{q}a_{s}a_{r}|D_{\nu}\rangle,$

where $p$, $q$, $r$, and $s$ are spin-orbital indices. Because the two
simulations are uncorrelated, the two-body reduced density matrix in Eq. (5)
becomes an unbiased one.

To calculate the spin-spin correlation, we need only the following 2-RDM
elements

$\displaystyle\Gamma_{i\uparrow(\downarrow)j\uparrow(\downarrow)}^{i\uparrow(\downarrow)j\uparrow(\downarrow)}$
$\displaystyle=$
$\displaystyle\langle\Psi|a^{\dagger}_{i,\uparrow(\downarrow)}a^{\dagger}_{j,\uparrow(\downarrow)}a_{j,\uparrow(\downarrow)}a_{i,\uparrow(\downarrow)}|\Psi\rangle$
$\displaystyle=$
$\displaystyle\langle\Psi|n_{i,\uparrow(\downarrow)}n_{j,\uparrow(\downarrow)}|\Psi\rangle$
$\displaystyle\Gamma_{i\uparrow(\downarrow)j\downarrow(\uparrow)}^{i\uparrow(\downarrow)j\downarrow(\uparrow)}$
$\displaystyle=$
$\displaystyle\langle\Psi|a^{\dagger}_{i,\uparrow(\downarrow)}a^{\dagger}_{j,\downarrow(\uparrow)}a_{j,\downarrow(\uparrow)}a_{i,\uparrow(\downarrow)}|\Psi\rangle$
$\displaystyle=$
$\displaystyle\langle\Psi|n_{i,\uparrow(\downarrow)}n_{j,\downarrow(\uparrow)}|\Psi\rangle$
$\displaystyle\Gamma_{i\uparrow(\downarrow)j\downarrow(\uparrow)}^{i\downarrow(\uparrow)j\uparrow(\downarrow)}$
$\displaystyle=$
$\displaystyle\langle\Psi|a^{\dagger}_{i,\uparrow(\downarrow)}a^{\dagger}_{j,\downarrow(\uparrow)}a_{j,\uparrow(\downarrow)}a_{i,\downarrow(\uparrow)}|\Psi\rangle.$
(6)

By using the following expressions of the local spin operators

$\displaystyle S_{i}^{x}$
$\displaystyle=\frac{1}{2}\left(a_{i\uparrow}^{\dagger}a_{i\downarrow}+a_{i\downarrow}^{\dagger}a_{i\uparrow}\right)$
$\displaystyle S_{i}^{y}$
$\displaystyle=\frac{i}{2}\left(a_{i\downarrow}^{\dagger}a_{i\uparrow}-a_{i\uparrow}^{\dagger}a_{i\downarrow}\right)$
$\displaystyle S_{i}^{z}$
$\displaystyle=\frac{1}{2}\left(n_{i\uparrow}-n_{i\downarrow}\right),$ (7)

the $\hat{S}^{z}$-spin correlation
$\langle\hat{S}_{i}^{z}\cdot\hat{S}_{j}^{z}\rangle$ can be evaluated as

$\displaystyle\langle$
$\displaystyle\hat{S}_{i}^{z}\cdot\hat{S}_{j}^{z}\rangle=\frac{1}{4}(\langle\Psi|n_{i,\uparrow}n_{j,\uparrow}|\Psi\rangle-\langle\Psi|n_{i,\uparrow}n_{j,\downarrow}|\Psi\rangle$
$\displaystyle-\langle\Psi|n_{i,\downarrow}n_{j,\uparrow}|\Psi\rangle+\langle\Psi|n_{i,\downarrow}n_{j,\downarrow}|\Psi\rangle)$
$\displaystyle=\frac{1}{4}(\Gamma_{i\uparrow j\uparrow}^{i\uparrow
j\uparrow}+\Gamma_{i\downarrow j\downarrow}^{i\downarrow
j\downarrow}-\Gamma_{i\uparrow j\downarrow}^{i\uparrow
j\downarrow}-\Gamma_{i\downarrow j\uparrow}^{i\downarrow j\uparrow}),$ (8)

and the total spin correlation function as

$\displaystyle\langle\hat{\bf S}_{i}\cdot\hat{\bf S}_{j}\rangle$
$\displaystyle=\frac{1}{2}(\Gamma_{i\uparrow j\downarrow}^{i\downarrow
j\uparrow}+\Gamma_{i\downarrow j\uparrow}^{i\uparrow j\downarrow})$
$\displaystyle+$ $\displaystyle\frac{1}{4}(\Gamma_{i\uparrow
j\uparrow}^{i\uparrow j\uparrow}+\Gamma_{i\downarrow j\downarrow}^{i\downarrow
j\downarrow}-\Gamma_{i\uparrow j\downarrow}^{i\uparrow
j\downarrow}-\Gamma_{i\downarrow j\uparrow}^{i\downarrow j\uparrow})$ (9)

For the special case when $\langle\hat{S}_{z}\rangle=0$, the above expressions
can be simplified to [38]

$\displaystyle\langle\hat{S}^{z}_{i}\cdot\hat{S}^{z}_{j}\rangle=\frac{1}{2}\left(\Gamma_{i\uparrow
j\uparrow}^{i\uparrow j\uparrow}-\Gamma_{i\uparrow j\downarrow}^{i\uparrow
j\downarrow}\right)$ (10)

and

$\displaystyle\langle\hat{\bf S}_{i}\cdot\hat{\bf
S}_{j}\rangle=\Gamma_{i\uparrow j\downarrow}^{i\downarrow
j\uparrow}+\frac{1}{2}(\Gamma_{i\uparrow j\uparrow}^{i\uparrow
j\uparrow}-\Gamma_{i\uparrow j\downarrow}^{i\uparrow j\downarrow})$ (11)

by using the relations $\Gamma_{i\uparrow j\uparrow}^{i\uparrow
j\uparrow}=\Gamma_{i\downarrow j\downarrow}^{i\downarrow j\downarrow}$ and
$\Gamma_{i\uparrow j\downarrow}^{i\uparrow j\downarrow}=\Gamma_{i\downarrow
j\uparrow}^{i\downarrow j\uparrow}$, see the Supporting Information (SI) and
Ref. [38] for details.

##  III results

Calculations are performed on three square lattices with tilted periodic
boundary conditions, namely the 18-site lattice with lattice vectors
$(3,3),(3,-3)$, the 20-site lattice with lattice vectors $(4,-2),(-2,-4)$, and
the 24-site lattice with lattice vectors $(5,1),(-1,-5)$. The supercell of the
24-site lattice is not square, but the underlying lattice still is. This
24-site lattice has been discussed comprehensively by Betts [39] with regard
to its finite-size properties, which are indeed favourable, having a low
topological imperfection as defined by Betts. The above three tilted lattices
are presented in [21] (see Fig.1 therein).

###  III.1 The two-hole system

Figure 1:  Spatial spin distributions given by
$\langle\hat{S}_{z}(i)\hat{S}_{z}(j)\rangle$ on the 18 sites lattice for two-
hole (a-c) and half-filled (d) systems. The critical interaction strength is
$U_{c}=92$ and the site $i$ is set in the center.

Figure 2: Spatial spin distributions given by
$\langle\hat{S}_{z}(i)\hat{S}_{z}(j)\rangle$ on the 20 sites lattice for two-
hole (a-c) and half-filled (d) systems.The critical interaction strength is
$U_{c}=103$.

Figure 3:  Spatial spin distributions given by
$\langle\hat{S}_{z}(i)\hat{S}_{z}(j)\rangle$ on the 24 sites lattice for two-
hole (a-c) and half-filled (d) systems. The critical interaction strength is
$U_{c}=127$. Figure 4:  Binding energy of two holes vs. 1/U on 18 sites Figure
5:  Spin distributions of the Slater determinant with the most weight in the
ground state of the 20 sites lattice with two holes. The red (blue) arrows
represent spin up (down), and the green dots represent a hole. The lattice
cutout is composed of four supercells.

Figure 6:  Statistical distribution of the distance between two holes on
18-(a), 20-(b) and 24-site(c) lattices. The number of considered SDs in the
ground-state wave function considered for the analysis are $5\times 10^{5}$
(0.93% of total occupied SDs in the FCIQMC ground state estimate), $8\times
10^{5}$(0.85% of occupied SDs), and $10^{7}$(0.76% of occupied SDs),
respectively. Please note that the ground-state wave function is not
normalized and the weight shown in the figure is a relative value.

The two-hole system has been investigated in our previous work [21], where we
find the total spin of the ground state is always zero for all the different
lattices. This supports the early conclusions based on exact diagonalization
[22] and analytical studies of arbitrarily large systems [19].

Table 1:  The ground state $S=0$ spin-spin correlations $\langle{\bf S}(i)\cdot{\bf S}(j)\rangle$ obtained from FCIQMC calculations on the 18-site square lattice Hubbard model with 16 electrons for different values of $U$. For comparison the values from the exact diagonalization(ED) calculation taken from Ref. [40] are also shown. $i$ and $j$ are the lattice site indices, and $R$ is their distance.  U R $\langle{\bf S}(i)\cdot{\bf S}(j)\rangle$ | 1 | $\sqrt{2}$ | 2 | $\sqrt{5}$ | 3
---|---|---|---|---|---
| 2-RDM | ED | 2-RDM | ED | 2-RDM | ED | 2-RDM | ED | 2-RDM | ED
20 | $-0.1910(1)$ | $-0.19096$ | $0.0448(1)$ | $0.04473$ | $-0.0154(1)$ | $-0.01541$ | $-0.01401(1)$ | $-0.01398$ | $0.0493(1)$ | $0.04931$
40 | $-0.1820(1)$ | $-0.18197$ | $0.0341(1)$ | $0.03417$ | $-0.0274(1)$ | $-0.02739$ | $-0.0062(1)$ | $-0.0061$ | $0.0619(1)$ | $0.06188$

To check the performance of our algorithm, we first perform calculations on
the 18-site lattice, where Lanzcos-based ED results are available [40]. In
Table 1, the results of the spin-spin correlation, $\langle\hat{\bf
S}(i)\cdot\hat{\bf S}(j)\rangle$, on the 18-site lattice for $U=20$ and $40$
are presented in comparison with the ED results. In our simulation we use a
time step of $\tau=0.001$, $N_{w}=5\cdot 10^{7}$ number of walkers and treat
the $5\cdot 10^{4}$ most populated states deterministically. As shown in Table
1, the two results agree very well for different $R$ (the distance between two
different sites $i$ and $j$).

In order to get a visual impression of the spatial spin distribution, we plot
the pair correlation function $\langle\hat{S}_{z}(i)\hat{S}_{z}(j)\rangle$ for
different interaction strengths and number of holes in a lattice cutout
composed of four supercells, as being presented in Figures 1, 2 and 3,
respectively for the 18, 20 and 24-sites lattices. The lattice site $i$ is set
in the center of the cutout.

On the 18-site lattice, the spatial spin distribution in the half-filling case
shows a clear anti-ferromagnetic pattern as shown in Fig. 1(d). For two-hole
systems, we still find some anti-ferromagnetic features when $U$ is not too
large, as the cases showing in Fig. 1(a) and 1(b), where spins on neighboring
sites are still anti-parallel. However, when $U$ is larger than a certain
critical value ($U_{c}$), as being presented in Fig. 1(c), the neighboring
spins become parallel, and we find some kind of ferromagnetic domains
structure. It is interesting to find that the critical $U_{c}$ is very close
to that of Nagaoka ferromagnetism in one-hole systems. We therefore do not try
to determine these critical $U_{c}$ values, but rather use those $U_{c}$’s
determined for the Nagaoka ferromagnetisation as references. Such transitions
of the spatial spin distribution patterns are also found on the 20- and
24-site lattices, see Fig. 2 and Fig. 3, where we also find ferromagnetic
domains when $U$ is large enough. When the lattice size increases, the size of
the ferromagnetic domains becomes larger, and their shape changes.

The fact that critical $U_{c}$’s are close to those of Nagaoka ferromagnetism
in one-hole systems, provides a physical picture that holes tend to separate
from each other. This picture, originally suggested by Tian [41], is now
tested and verified based on the analysis of Slater determinants of the
ground-state. For the two-hole systems, we have also found that for large
$U$’s the binding energies are positive as shown in Fig. 4. The binding energy
of two holes is defined as
$\delta=(E_{2}-E_{0})-2(E_{1}-E_{0})=E_{2}-2E_{1}+E_{0}$, where $E_{n}$ is the
ground-state energy of the $n$-hole system. The positive binding energy
implies that two holes tend to separate, rather than to bind for large $U$’s.
In this case, it is reasonable to assume that each hole carries a halo of
ferromagnetic texture.

To further confirm the picture of the hole separations we have looked at the
spin distribution of the most populated SD in the ground-state of the 20-site
lattice. For $U=200$, the spin distribution of the SD with highest weight is
plotted in Fig. 5, where we find that the distance between two holes is the
farthest, $\sqrt{10}$, and spins around each hole are the same. The analysis
of the $8\times 10^{5}$ most populated SDs shows that the most common hole-
distance is not $\sqrt{10}$, but $\sqrt{5}$ instead (as shown in Fig. 6(b)).
Similarly on the 24-site lattice, $\sqrt{5}$ is the most common hole-distance
(Fig. 6(c)), and on the 18-site lattice, the distance $2$ is slightly more
common than $\sqrt{5}$ (Fig. 6(a)). Ferromagnetic domains appear on finite
lattices when holes are far away from each other. For systems with a fixed
number of holes, the possibility of the formation of ferromagnetic domains
increases with the system size, which means that the ferromagnetic domain
structure favors a low hole density.

###  III.2 The three-hole system

Figure 7:  $E(S)-E(S_{max})$ versus $S$ on 16- (a) and 18- (b) sites lattices
with three holes, and the width of spin spectrum $\Delta E$ as the function of
$1/U$ on 18- (c) and 20- (d) sites lattices, respectively.

Figure 8:  The spatial spin distribution of the 18-, 20- and 24-sites lattices
with 3 holes are provided in (a-c), $U_{c}$ are the same as those in two-hole
system.

For three-hole systems, we first study the possibility of partial spin-
polarization on lattices larger than 16 sites. As been mentioned in the
introduction, results of exact diagonalization (ED) based on an effective
Hamiltonian [22] show that the total spin of the ground-state on 8- and
16-sites lattices are $\frac{3}{2}$ and $\frac{7}{2}$ respectively. If this
partial spin-polarization also exists for larger lattices, then for large
enough $U$ there would be a phase of ferrimagnetism (a mixed type of anti-
ferromagnetism and ferromagnetism). To investigate this problem, we have
mainly calculated ground state energy as a function of total spin, $E(S)$, on
16- and 18-sites lattices.

In Fig. 7(a), the ground state energy as a function of the total spin for the
16-site lattice is shown. The results are presented for different $U$’s in
comparison with the ED result for $U=\infty$ [22], according to Eq. (2). With
increasing $U$, the $E(S)$ curves clearly converges to the ED result, and in
fact the result at $U=20000$ already coincides with the ED result. These
results show a partial spin polarization at $S=\frac{7}{2}$ in the large $U$
regime. However, on the 18-sites lattice, such a partial spin-polarization
does not exist for any $U$. In Fig. 7(b), the results of $E(S)$ are plotted
for different $U$’s for the 18-site lattice. The coincidence of the two curves
at $U=20000$ and $U=50000$ indicates the convergence in the large $U$ limit.
For any $U$ the total spin of the ground-state always take the smallest value
$S=0.5$ and thus there is no partial spin polarization. Similar results are
also obtained on the 20-site lattice, where we have only performed
calculations for $U=20000$. We therefore believe that the partial spin
polarization result obtained in Ref. [22] is only a finite size effect.

To investigate the possibility of ferromagnetic domains in three-hole systems
we have calculated the width of the spin spectrum $\Delta E=E_{max}-E_{min}$,
where $E_{max}=\max(E(S)),E_{min}=\min(E(S))$ are the maximal and minimal
values of the energy over all spin states for a given $U$. In our previous
work [21], we find that for two holes systems a clear change in the slope of
$\Delta E(1/U)$ usually implies the formation of ferromagnetic domains. The
width of spin spectrum $\Delta E$ for three holes systems on the 18- and
20-sites lattice are plotted in Fig. 7(c) and (d) respectively as functions of
$1/U$, in comparison with those for two holes systems. For three-hole systems,
the curve of $\Delta E(1/U)$ for the 18-sites lattice is simply a straight
line, while for the 20-sites lattice we see a small but clear change of slope.
This indicates that on large enough lattices the three holes systems can also
have ferromagnetic domains in the large $U$ regime. To support this, we have
also investigated spatial spin distributions of three holes systems by
calculating the spin-spin correlation functions. The spatial spin
distributions are presented in Fig. 8, where ferromagnetic domains are found
on 20- and 24- site lattices. Comparing with the ferromagnetic domains in two-
hole systems, the extent of ferromagnetic domains is smaller in this case.

##  IV Conclusion

The instability of Nagaoka ferromagnetism in the Hubbard model is studied for
two and three holes systems on finite lattices. It is shown that the total
spin of the ground-state takes the minimal value $S=S_{min}$, and there exists
no partial spin polarization on lattices larger than 16 sites both on two- and
three-hole systems. In the large $U$ regime ($U\geq U_{c}$, where $U_{c}$ is
the critical $U$ of Nagaoka ferromagnetism in the one-hole system), we find
the formation of ferromagnetic domains. Based on the analysis of binding
energy and the statistical distance between two holes, we find a general
feature of the ferromagnetic domain structure, that holes tend to be far away
from each other.

Acknowledgments

Sujun Yun would like to extend her sincere gratitude to Nikolay Bogdanov and
Giovanni Li Manni of the Max Planck Institute for Solid State Research, Youjin
Deng in the University of Science and Technology of China and Qianghua Wang in
Nanjing University for participating in discussion and providing support. The
authors gratefully acknowledge funding from the Max Planck Society. Sujun Yun
is supported by the national natural science foundation of China (Grant No.
11805103 and No. 11447204) and Jiangsu natural science foundation (Grant No.
BK20190137).

## References

  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, Annual Review of Condensed Matter Physics 13, 275 (2022).
  * Scalapino [2012] D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * Tohyama  _et al._ [2005] T. Tohyama, Y. Inoue, K. Tsutsui, and S. Maekawa, Phys. Rev. B 72, 045113 (2005).
  * Ohki  _et al._ [2018] D. Ohki, G. Matsuno, Y. Omori, and A. Kobayashi, Crystals 8, 137 (2018).
  * Dai  _et al._ [2022] H. Dai, M. Cai, Q. Hao, Q. Liu, Y. Xing, H. Chen, X. Chen, X. Wang, H.-H. Fu, and J. Han, ACS NANO 16, 12437 (2022).
  * Holzmann and Moroni [2020] M. Holzmann and S. Moroni, Phys. Rev. Lett. 124, 206404 (2020).
  * Nagaoka [1966] Y. Nagaoka, Phys. Rev 147, 392 (1966).
  * Thouless [1965] D. J. Thouless, Proceedings of the Physical Society 86(5), 893 (1965).
  * Lieb [1989] E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  * Tasaki [1989] H. Tasaki, Phys. Rev. B 40, 9192 (1989).
  * Tasaki [1998] H. Tasaki, J. Phys. Condens. Matter 10, 4353 (1998).
  * Ivantsov  _et al._ [2017] I. Ivantsov, A. Ferraz, and E. Kochetov, Phys. Rev. B 95, 155115 (2017).
  * Barbieri  _et al._ [1990] A. Barbieri, J. Riera, and A. Young, Phys. Rev. B 41, 11697 (1990).
  * Becca and Sorella [2001] F. Becca and S. Sorella, Phys. Rev. Lett. 86, 3396 (2001).
  * Liu  _et al._ [2012] L. Liu, H. Yao, E. Berg, S. R. White, and S. A. Kivelson, Phys. Rev. Lett. 108, 126406 (2012).
  * Müller-Hartmann  _et al._ [1993] E. Müller-Hartmann, T. Hanisch, and R. Hirsch, Physica B: Condensed Matter 186-188, 834 (1993).
  * Wurth  _et al._ [1996] P. Wurth, G. Uhrig, and E. Müller-Hartmann, Annalen der Physik 508, 148 (1996).
  * Doucot and Wen [1989] B. Doucot and X. G. Wen, Phys. Rev. B 40, 2719(R) (1989).
  * Fang  _et al._ [1989] Y. Fang, A. E. Ruckenstein, E. Dagotto, and S. Schmitt-Rink, Phys. Rev. B 40, 7406(R) (1989).
  * Yun  _et al._ [2021] S. Yun, W. Dobrautz, H. Luo, and A. Alavi, Phys. Rev. B 104, 235102 (2021).
  * Riera and Young [1989] J. A. Riera and A. P. Young, Phys. Rev. B 40, 5285(R) (1989).
  * Foulkes  _et al._ [2001] W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Shi and Zhang [2021] H. Shi and S. Zhang, J. Chem. Phys. 154, 024107 (2021).
  * Booth  _et al._ [2009] G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem. Phys. 131, 054106 (2009).
  * Cleland  _et al._ [2010] D. Cleland, G. H. Booth, and A. Alavi, The Journal of Chemical Physics 132, 041103 (2010).
  * Guther  _et al._ [2020] K. Guther, R. J. Anderson, N. S. Blunt, N. A. Bogdanov, D. Cleland, N. Dattani, W. Dobrautz, K. Ghanem, P. Jeszenszki, N. Liebermann, G. L. Manni, A. Y. Lozovoi, H. Luo, D. Ma, F. Merz, C. Overy, M. Rampp, P. K. Samanta, L. R. Schwarz, J. J. Shepherd, S. D. Smart, E. Vitale, O. Weser, G. H. Booth, and A. Alavi, The Journal of Chemical Physics 153, 034107 (2020).
  * Overy  _et al._ [2014] C. Overy, G. H. Booth, N. S. Blunt, J. S. Spencer, D. Cleland, and A. Alavi, J. Chem. Phys. 141, 244117 (2014).
  * E. Thomas  _et al._ [2015] R. E. Thomas, G. H. Booth, and A. Alavi, Phys. Rev. Lett. 114, 033001 (2015).
  * Blunt  _et al._ [2017] N. S. Blunt, G. H. Booth, and A. Alavi, J. Chem. Phys. 146, 244105 (2017).
  * Paldus [1976] J. Paldus, Phys. Rev. A 14, 1620 (1976).
  * Shavitt [1977] I. Shavitt, Int. J. Quantum Chem. 11, 131 (1977).
  * Shavitt [1978] I. Shavitt, Int. J. Quantum Chem. 12, 5 (1978).
  * Dobrautz  _et al._ [2019] W. Dobrautz, S. D. Smart, and A. Alavi, J. Chem. Phys. 151, 094104 (2019).
  * Dobrautz  _et al._ [2021] W. Dobrautz, O. Weser, N. A. Bogdanov, A. Alavi, and G. Li Manni, Journal of Chemical Theory and Computation 17, 5684 (2021), pMID: 34469685\.
  * Manni  _et al._ [2020] G. L. Manni, W. Dobrautz, and A. Alavi, Journal of Chemical Theory and Computation 16, 2202 (2020).
  * Li Manni  _et al._ [2021] G. Li Manni, W. Dobrautz, N. A.Bogdanov, K. Guther, and A. Alavi, J. Phys. Chem. A 125, 4727 (2021).
  * Kutzelnigg  _et al._ [2010] W. Kutzelnigg, K. Shamasundar, and D. Mukherjee, Mol. Phys. 108, 433 (2010).
  * Betts  _et al._ [1999] D. D. Betts, H. Q. Lin, and J. S. Flynn, Can.J.Phys. 77, 353 (1999).
  * Becca  _et al._ [2000] F. Becca, A. Parola, and S. Sorella, Phys. Rev. B 61, R16287 (2000).
  * Tian [1991] G. Tian, Phys. Rev. B 44, 4444 (1991).


# Black hole mass estimate in OJ 287 based on the bulk-motion comptonization
model

Sergey Kuznetsov1 & Lev Titarchuk2
1Prosto Science AI Lab, Montreal, QC, H2X 3X2, Canada
2Dipartimentodi Fisica, Università di Ferrara, via Saragat 1, 44122 Ferrara,
Italy E-mail<EMAIL_ADDRESS>(SK)

(Last updated 2024 February 5; in original form 2023 January 24)

###### Abstract

The multi-wavelength outburst activity in the BL Lacertae source OJ 287 has
sparked a lot of controversy about whether the source contains one or two
black holes (BHs) and what characteristics of this black hole binary would be.
In this article we present the results of analysis of the X-ray flaring
activity of OJ 287 using the data of Swift/XRT observations. We discovered
that the energy spectra in all spectral states can be adequately fit with the
XSPEC BMC model (the Comptonization one). As a result we found that the X-ray
photon index of the BMC model, $\Gamma$ correlates with the mass accretion
rate, $\dot{M}$. We found the photon index $\Gamma$ to increase monotonically
with accretion rate $\dot{M}$ from $\Gamma\sim 2$ in the intermediate state
(IS) to $\Gamma\sim 2.5$ the high/soft state (HSS) with subsequent saturation
at $\Gamma\sim$ 2.6 level at higher luminosities. This type of behavior of the
spectral index is remarkably similar to the pattern observed in a number of
established stellar-mass black hole candidates. Assuming the universality of
the observed pattern of the correlation between the photon index and mass
accretion rate, we estimate the BH mass in OJ 287 to be around $2\times
10^{8}$ solar masses, using the well studied BH binaries GX 339–4 and XTE
J1859–226 to calibrate the model.

###### keywords:

black hole physics—BL Lacertae objects:individual:OJ 287–galaxies:active

††pubyear: 2023††pagerange: Black hole mass estimate in OJ 287 based on the
bulk-motion comptonization model–References

##  1 Introduction

The primary goal of our paper is to examine the value of the secondary black
hole (BH) mass in OJ 287 applying the scaling method of BH mass determination
[Shaposhnikov & Titarchuk (2009), hereafter ST09]. OJ 287 demonstrated a
noticeable 12 yr cycle [Valtonen et al. (2006)] that can be explained by the
presence of the accreting black hole binary system (BHBS) with systematic
disturbance of the accretion disk (AD) around the primary by the secondary.

Valtonen et al. (2012), hereafter V12, presented results of the two X-ray
observations of OJ 287 by XMM-Newton on April 12 and November 3-4, 2005. In
Fig. 2 they plotted the observations of the 2005 campaigns. V12 claimed that
for the April 12, 2005 the spectral energy distribution (SED), the spectrum
from radio to X-rays followed nicely a synchrotron self-Compton (SSC) model.
V12 further stated that the BH mass value of the secondary determined from the
orbit solution $M_{sec}\sim 1.4\times 10^{8}$ $M_{\odot}$ agrees with the
measured properties of the optical outburst. V12 also formulated a question
regarding the primary BH mass necessary to guarantee the stability of the
primary accretion disc. They found that the minimum value of the primary mass
$1.8\times 10^{10}M_{\odot}$ is quite close to the BH mass determined from the
orbit, $1.84\times 10^{10}M_{\odot}$.

To find evidence for the emission of the secondary, V12 needed to look at the
short time-scale variability in OJ 287. It has been found to be variable from
15 min time-scale upwards on many occasions [see for example, Gupta et al.
(2012)] and on one occasion the light curve has shown sinusoidal variations
with a period of 228 min (Sagar et al., 2004). If this variation is associated
with the last stable orbit of a maximally rotating black hole, the mass of the
black hole would be $1.46\times 10^{8}M_{\odot}$ (Gupta et al., 2012), i.e.
identical to the mass obtained from the orbit solution (Valtonen et al.,
2010).

Komossa et al. (2021), hereafter KOM21 made a detailed analysis of XMM-Newton
spectra of OJ 287 spread over 15 years. KOM21 also presented data from Swift
UVOT and XRT observations of OJ 287, which began in 2015, along with all
public Swift information after 2005. During this period, OJ 287 was found in
the low/hard extreme state and outburst soft states. In addition, they
established that the OJ 287 X-ray spectrum was highly variable, passing all
states seen in blazars from the low/hard (flat) state to exceptionally soft
steep state (ST). KOM21 found that these spectra can consist of two parts: IC
radiation which is prevailing at the low-state, and very soft radiation
component that becomes extremely powerful when OJ 287 becomes more luminous.
KOM21 have concluded that their 2018 XMM-Newton spectra almost identical with
the EHT examination of OJ 287, and were well characterized by a model with a
hard IC component of the photon index $\Gamma\sim 1.5$ and a soft component.

We reanalyzed the data of the X-ray telescope (XRT) onboard of Swift
Observatory for OJ 287 and applied the ST09 scaling method to this data in
order to make a BH mass estimate. Using this method requires an accurate
knowledge of the distance to the source OJ 287. The method was proposed back
in 2007 by Shaposhnikov & Titarchuk (2007), hereafter ST07 and also by ST09.
It is worth noting that there are two scaling methods: the first method based
on the correlation between the photon index $\Gamma$ and the quasi-periodic
oscillation frequency (QPO) $\nu_{L}$; and the second method is based on the
correlation between $\Gamma$ and normalization of the spectrum proportional to
$\dot{M}$. If the first method ($\Gamma-\nu_{L}$) is used to estimate the mass
of the BH, the distance to the source is not required (ST07), but for the
second method ($\Gamma-\dot{M}$) the source distance and the inclination of
the accretion disk relative to the Earth observer are needed (ST09).

For both methods, it is important that the source shows a change in spectral
states during the outburst and a characteristic behavior of $\Gamma$, namely,
a monotonic increase of $\Gamma$ with $\nu_{L}$ or $\dot{M}$ in the
LHS$\to$IS$\to$HSS transition and reaching a constant level (saturating) at
high values of $\nu_{L}$ or $\dot{M}$, followed by a monotonic decrease in
$\Gamma$ during the HSS$\to$IS$\to$LHS transition during the outburst decay.

The saturation of $\Gamma$ (so-called “$\Gamma$-saturation phase”) during the
outburst is a specific signature that this particular object contains a BH
(Titarchuk & Zannias, 1998). Indeed, the $\Gamma$-saturation phase can be
caused only by an accretion flow converging to the event horizon of a BH [see
the Monte-Carlo simulation results in Laurent & Titarchuk (1999, 2011)].
Because of this, it makes sense to compare BH sources that have the same
$\Gamma$-saturation levels. In the second method ($\Gamma-\dot{M}$), it is
assumed that BH luminosity is directly proportional to $\dot{M}$ (and,
consequently, to the mass of the central BH), and inversely proportional to
the squared distance to the source. Thus, we can determine the BH mass by
comparing the corresponding track $\Gamma-\dot{M}$ for a pair of sources with
BHs, in which all parameters are known for one source, and for the other
source parameters except a BH mass are known (for more details on the scaling
method, see ST09).

The scaling method has a number of advantages in determination of a BH mass
compared to other methods, namely, the calculation of the X-ray spectrum
originating in the innermost part of disk based on first-principle
(fundamental) physical models, taking into account the Comptonization of the
soft disk photons by hot electrons of the internal disk part and in a
converging flow to BH.

It is necessary to emphasize that our approach is only justifiable if the
observed X-ray emission in OJ 287 is indeed produced by BMC in the accretion
flow. In reality, it may be dominated by emission produced in relativistic
jets, which is a popular explanation of the multi-wavelength emission
properties of blazars in general. Nevertheless, we apply the scaling method to
OJ 287 because, as will be shown below, it yields results (specifically, the
spectral slope vs. normalization dependence) that indicate that at least a
significant fraction of the observed X-ray emission in this blazar is produced
by the BMC mechanism.

In this paper, based on Swift data analysis, we estimate a BH mass in OJ 287
by applying the scaling technique. In §2 we provide details of our data
analysis while in §3 we present a description of the spectral models used for
fitting of this data. In §4 we focus on observational results and their
interpretation. In §5 we discuss the main results of the paper. In §6 we
present our final conclusions.

##  2 Observation and Data Reduction

We analyzed observational data of OJ 287 carried out by XRT instrument aboard
Swift observatory from April 2016 to June 2017 (MJD 57501-57918). We processed
$\sim$140 ks of pointed observations of the source. We retrieved the data from
the UK Swift Science Data Centre Archive111https://www.swift.ac.uk/archive and
used XRT data both in Windowed Timing (WT) and Photon Counting (PC) modes.

The lowest source count rate was detected at the very beginning and at the end
of this observational cycle (April–June 2016, May–June 2017). It corresponds
approximately to XRT observations in PC mode with the flux below $\sim$1
mCrab. This covers 31 PC observations of $\sim$34 ks in total. The rest of the
observations were carried out during an epoch of strong flaring activity of OJ
287 in the middle of this period from October 2016 to April 2017 and they
consisted of $\sim$ 105 ks of 109 WT observations. We excluded only 3 very
short observations from our analysis with accumulation time less than 100 s
per ObsID and an additional 6 ones when the source position on the detector
was at the edge of the field of view (FOV).

For our comprehensive analysis we reprocessed original Level 1 XRT data using
standard xrtpipeline screening settings from HEASOFT package v6.30 and
calibration files version 20210915. We also checked PC data for the pile-up
effect and found its possible presence in 2 observations with contamination
radius of 2 detector pixels.

The source spectra were extracted from every individual PC and WT observations
for our preliminary analysis. For the PC mode we used a 30-pixel circular
region for the source spectrum and a 60-pixel circular region taken from
source-free area for background correction. Because of the lack of 2D-image in
WT mode, we used a rectangular 20x40 pixels region co-aligned to the axis of
1D-image for both source and background spectra.

A typical pointed XRT observation of OJ 287 consists of $\sim$1 ks or less. As
such, a short duration did not allow us to obtain statistically significant
best-fit spectral model parameters. At the same time, the estimation of the
source flux by the applicable model to the spectrum can be a tool for the
further averaging of individual source spectra over the flux. In order to
implement this approach we used the tbabs*bmc (Bulk motion Comptonization)
model with fixed $N_{H}=2.7\times 10^{20}$ and log$(A)=2.0$. The applicable
above absorption toward OJ 287 was adopted from Willingale et al.
(2013)222https://www.swift.ac.uk/analysis/nhtot and constant log$(A)$ was used
for simplification and reduction of free model parameters.

The preliminary spectral analysis demonstrated the 0.3–10 keV source flux
variations in the following ranges: 0.5$-$2.6$\times 10^{-11}$ and
0.5$-$4.1$\times 10^{-11}$ erg s-1 cm-2 for PC and WT observations. Since our
goal was the averaging of the source spectra with similar flux, we grouped 31
PC and 100 WT observations into 7 and 18 flux bins correspondingly.

The procedure of retrieving averaged spectrum is trivial for PC observations
and includes the following steps: reading of re-processed event files,
application of region filters for source and background spectra extraction
with xselect; building the accumulated exposure map with ximage; and
generation of Ancillary Response Files (ARFs) with xrtmkarf. Similar to the
preliminary spectral analysis described above per individual ObsID, we applied
the same 30/60-pixel circular regions for the source extraction and background
correction.

The retrieving of averaged source spectrum over WT observations is not trivial
if the orientation of the detector does not match among all observations plus
if the source positions on the detector were shifted. As a result, the
averaged 2D-image will be a combination of 1D-images of the same length under
different angles and with different shifts centered to the source position.
The hardest task in this procedure is the proper background extraction and
source/background area correction using event lists and builded exposure maps.

In order to extract the source spectrum from the set of accumulated WT
observations, we used a 20-pixel circular region as a filter with xselect. In
this case we extracted the source counts using the same 40-pixel width/length
region as we did before for individual WT observations (see above). For
background extraction we used the annulus centered to the source position with
inner/outer radii of 80/120 pixels. Since the detector size is only 200 pixels
in width, we are then getting the same 40 linear pixels of 1D-image on the
detector despite the individual 1D-images orientations and shifts. The last
step in getting of WT averaged spectrum is the manual correction of the
background scaling factor stored in BACKSCAL keyword of the spectral file.
Since the obtained background effective surface is actually the length and it
has the same width of 40-pixels as the source, we then manually matched this
value of BACKSCAL keyword in the header of the background file to the source
spectrum one.

##  3 Analysis and Results

For our comprehensive spectral analysis we fitted 7 PC and 18 WT averaged over
the flux spectra by tbabs*bmc model with fixed $N_{H}=2.7\times 10^{20}$ (see
above). The typical OJ 287 spectra with $\Gamma=2.0$ and $\Gamma=2.6$ are
shown in Fig. 1.

For our further analysis we plotted the spectral index $\Gamma$ vs $N$ in Fig.
2. The correlation between $\Gamma$ and $N$ is noticeable for $N$ lower than
$\sim 2.5\times 10^{-4}$. Taking into account that $\Gamma$ does not fluctuate
significantly at higher $N$ values, we separated our data points shown in Fig.
2 in two parts and fitted them by linear function and a constant at lower and
higher $N$ correspondingly. The result of this best-fit approximation is shown
in Fig. 2 by a solid line.

Additionally, we fitted the same data points using the analytic function

$F(N)=A-B\ln\left[\exp\left(1.0-(N/N_{tr})^{\beta}\right)+1\right]$ (1)

widely used as scalable pattern between reference and target sources in black
hole mass estimation (see ST09). The results of this approximation are shown
in Table 1. GX 339-4 best-fit parameters corresponding to observed LHS to IS
source transition in 2007 were taken from ST09 Fig. 3 and Table 2 (28 data
points, online extended table). XTE J1859-226 data were adopted from ST09 (see
Fig. 9 there).

Taking into account $F(N)$ tends asymptotically to a constant at infinity, we
projected the saturation value of $N$ we found for OJ 287 to $F(N)$. In Fig. 2
it corresponds to the projection of the break point of the solid curve to the
dotted curve. We then repeated the same approach in scaling $F(N)$ for GX
339-4 and XTE J1859-226. For the reference sources $N$ should be the same way
equidistant from the upper limit of $F(N)$ as $N$ for the target source. In
other words, the relative height of $N$ position between the asymptotes should
be the same.

The curves corresponding to analytical function $F(N)$ of the reference and
target sources from Table 1 and their $N_{t}$ and $N_{r}$ are shown in Fig. 3
in details.

Table 1: Parameters of scaling patterns of Eq. (1) for reference and target sources Source | $A$ | $B$ | $N_{tr}$ | $\beta$
---|---|---|---|---
XTE J1859-226 | 2.55 | 1.25 | 7.13E-2 | 1.00
GX 339-4 | 2.46 | 0.955 | 1.28E-1 | 1.60
OJ 287 | 2.54 | 0.42 | 1.33E-4 | 2.28
Figure 1:  X-ray observed spectra of OJ 287, using the tbabs*bmc fit model for
the high/soft (dotted line ) and intermediate (dashed line) states. Spectra
are presented as $E^{2}f(E)$-diagram where $f(E)$ is a photon spectrum. It is
evident using this plot that the high/soft luminosity is at least four times
higher than that in the intermediate state.

###  3.1 X-ray spectra

To fit the energy spectra of this source we used a XSPEC model consisting of
the Comptonization (bulk-motion Comptonization, hereafter BMC) component [see
Titarchuk & Zannias (1998); Laurent & Titarchuk (1999)]. We also used a
multiplicative tbabs model (Wilms et al., 2000) which takes into account
absorption by neutral material. We assume that accretion onto a BH is
described by two main zones [see, for example, Fig. 1 in Titarchuk, & Seifina
(2021)]: a geometrically thin accretion disk [e.g. the standard Shakura-
Sunyaev disk, see Shakura & Sunyaev (1973) and a transition layer (TL), which
is an intermediate link between the accretion disk, and a converging (bulk)
region (see Titarchuk & Fiorito (2004)], that is assumed to exist, at least,
below 3 Schwarzschild radii, $3R_{S}=6GM_{\rm BH}/c^{2}$. The spectral model
parameters are the equivalent hydrogen absorption column density $N_{H}$; the
photon index $\Gamma$; $\log(A)$ is related to the Comptonized factor $f$
[$={A}/{(1+A)}$]; the color temperature and normalization of the seed photon
blackbody component, $kT_{s}$ and $N$, respectively.

Similarly to the ordinary bbody XSPEC model, normalization is a ratio of the
source (disk) luminosity $L$ to the square of the distance $d$ (ST09, see Eq.
1 there):

$N=\biggl{(}\frac{L}{10^{39}\mathrm{erg/s}}\biggr{)}\biggl{(}\frac{10\,\mathrm{kpc}}{d}\biggr{)}^{2}$
(2)

This encompasses an important property of our model. Namely, using this model
one can correctly evaluate normalization of the original “seed” component,
which is presumably a correct $\dot{M}$ indicator (Seifina & Titarchuk, 2011).
In its turn

Figure 2:  $\Gamma$ vs $N$ (normalization, proportional to the mass accretion
rate) correlations those for OJ 287 (target source). We use a specific
function for $\Gamma=1.80+3.190\times N$ for $N<2.3\times 10^{-4}$ and
$\Gamma=2.53$ for $N>2.3\times 10^{-4}$ to fit the data. Dotted curve
represents best-fit model shape of Eq 1.

$L=\frac{GM_{BH}\dot{M}}{R_{*}}=\eta(r_{*})\dot{m}L_{Ed}$ (3)

Here $R_{*}=r_{*}R_{S}$ is an effective radius where the main energy release
takes place in the disk, $R_{S}=2GM/c^{2}$ is the Schwarzschild radius,
$\eta=1/(2r_{*})$, $\dot{m}=\dot{M}/\dot{M}_{crit}$ is the dimensionless
$\dot{M}$ in units of the critical mass accretion rate
$\dot{M}_{crit}=L_{Ed}/c^{2}$, and $L_{Ed}$ is the Eddington luminosity.

For the formulation of the Comptonization problem, one can see Titarchuk et
al. (1997); Titarchuk & Zannias (1998); Laurent & Titarchuk (1999); Borozdin
et al. (1999); Shaposhnikov & Titarchuk (2009).

Spectral analysis of the Swift/XRT data of OJ 287, in principle, provides a
general picture of the spectral evolution. We can trace the change in the
spectrum shape during the IS–HSS transition in Fig. 1, which demonstrates two
representative $E*F_{E}$ spectral diagrams for different states of OJ 287. The
emergent spectra (see Fig. 1) can be described as a sum of the low energy
blackbody and its fraction convolved with the Comptonization Green function
(CGF) [see Eqs. (16) and (B5) in Sunyaev & Titarchuk (1980)]. The HSS and IS
spectra are characterized by a strong soft blackbody component (presumably,
related to the accretion disk) and a power law (as the hard tail of the CGF),
while in the LHS the Comptonization component is dominant and the blackbody
component is barely seen because the innermost part of the disk is fully
covered by the scattering media which has a Thomson optical depth more than 2.

Analysis of the Swift/XRT data fits (see Fig. 1) showed that $\Gamma$
monotonically increases from 2 to 2.6, when normalization of the spectral
component (or $\dot{M}$) increases by a factor of about 5 (see points, in Fig.
2) at the outburst rise phase (IS–HSS).

Table 2: BH masses and distances

Reference sources |  $m_{r}^{(a)}$ (M${}_{\odot})$ |  $i_{r}^{(a)}$ (deg) |  $d_{r}^{(b)}$ (kpc) | References
---|---|---|---|---
XTE 1859–226 | 5.4(a) | $\sim 70$ | $\sim 4.2^{(b)}$ |  Corral-Santana et al. (2011), Hynes et al. (2002), Shaposhnikov & Titarchuk (2009)
GX 339-4 | 12(a) | $\sim 70$ | 5.75$\pm 0.8^{(a)}$ |  Shaposhnikov & Titarchuk (2009), Hynes et al. (2003), Heida et al. (2017)
Target source |  $m_{t}$ (M⊙) |  $i_{t}^{(a)}$ (deg) |  $d_{t}^{(b)}$ (kpc) |
OJ 287 | 2$\times 10^{8}$ | 50 | 1.037 Gpc | using XTE J1859–226 as a reference source
OJ 287 | 2$\times 10^{8}$ | 50 | 1.037 Gpc | using GX 339-4 as a reference source

(a) Dynamically and using the scaling method values of BH masses and system
inclination; (b) source distance found in literature.

Figure 3:  $\Gamma$ vs. N (normalization, proportional mass accretion rate and
BH mass) for OJ 287 (target source) and GX 339-4 and XTE 1859-226, the
reference sources. The transition from the power line (linear shape) to the
asymptotic horizontal asymptotic line occurs at $N_{t}=0.23\times 10^{-4}$ for
OJ 287 and at $N_{r}=0.32,~{}0.31$ for GX 339-4 and XTE 1859-226,
respectively.

##  4 A BH mass estimate

Now, we use the scaling method for a BH mass estimate ($M_{oj287}$) in OJ 287.
The BH mass scaling method using the $\Gamma-N$ correlation is described in
detail in ST09. This method consists of (i) searching for such a pair of BHs
for which the $\Gamma$ correlates with increasing normalization $N$ (which is
proportional to mass accretion rate $\dot{M}$ and a BH mass $M$, see ST09, Eq.
(7) there) and the saturation level, $\Gamma_{sat}$ are the same and (ii)
calculating the scaling coefficient $s_{N}$ which allows to determine a BH
mass of the target object. It is worthwhile to emphasize that one needs a
ratio of distances for the target and reference sources in order to estimate a
BH mass using the following equation for the scaling coefficient

$s_{N}=\frac{N_{r}}{N_{t}}=\frac{m_{r}}{m_{t}}\frac{d_{t}^{2}}{d_{r}^{2}}{f_{G}}$
(4)

where $N_{r}$, $N_{t}$ are normalizations of the spectra,
$m_{t}=M_{t}/M_{\odot}$, $m_{r}=M_{r}/M_{\odot}$ are the dimensionless BH
masses with respect to solar, $d_{t}$ and $d_{r}$ are distances to the target
and reference sources, correspondingly. A geometry factor, $f_{G}=\cos
i_{r}/\cos i_{t}$ where $i_{r}$ and $i_{r}$ are the disk inclinations for the
reference and target sources, respectively [see ST09, Eq. (7)].

We found that XTE 1859–226 and GX 339-4 can be used as the reference sources
because these sources met all aforementioned requirements to estimate a BH
mass of the target source OJ 287 [see items (i) and (ii) above].

In Fig. 1 we show two types of the spectral states observed in OJ 287. The
first one (in the bottom there) is a typical intermediate spectrum observed
when the resulting luminosity is relatively low. We presented all emergent
spectra of this source in terms of $EF_{(}E)$ diagrams in order to estimate
the resulting luminosity. In Fig. 2 we demonstrate the $\Gamma$ versus $N$
where $N$ is presented in the units of $L_{39}/d^{2}_{10}$ ($L_{39}$ is the
source luminosity in units of $10^{39}$ erg/s and $d_{10}$ is the distance to
the source in units of 10 kpc). As one can see, the correlations of these
sources are characterized by similar shapes and saturation levels
$\Gamma_{sat}>2.5$. In order to implement the scaling method we use an
analytical approximation $F(N)$ for the $\Gamma-N$ correlation,

$\displaystyle\Gamma(N)$ $\displaystyle=$ $\displaystyle
1.84+3081*N~{}~{}~{}~{}~{}~{}~{}{\rm for}~{}~{}~{}~{}~{}~{}~{}N<2.3\times
10^{-4}$ $\displaystyle\Gamma(N)$ $\displaystyle=$ $\displaystyle
2.53~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm
for}~{}~{}~{}~{}~{}~{}~{}N>2.3\times 10^{-4}$ (5)

We estimated a BH mass for OJ 287 using Fig 2 of this manuscript and Figs. 3
and 9 in ST09, for GX 339-4 and XTE1859-226 respectively. Thus, we found that
the $\Gamma$ vs $N$ correlations are self-similar, at least for $\Gamma>2$,
see Fig. 3 and from these correlations we could estimate $N_{t}$, $N_{r}$ for
OJ 287 and GX 339-4 and XTE1859-226, respectively. They are $N_{t}=2.3\times
10^{-4}$ and $N_{r,GX}=0.27$, $N_{r,1859}=0.24$. taken at the beginning of the
$\Gamma$-saturation part (see solid points).

Thus, for GX 339-4 we obtained $s_{N}=N_{r}/N_{t}\sim 0.1\times 10^{4}$
applying $N$-values: $N_{r}=0.27$ and $N_{t}=2.3\times 10^{-4}$ (see Fig. 3).
On the other hand, a value of a factor $f_{G}=\cos{i_{r}}/\cos{i_{t}}\sim
0.53$ for the target and reference sources can be obtained using inclinations
$i_{t}=50^{o}$ and $i_{r}=70^{o}$ (see Table 2 and Hynes et al. (2003)). As a
result of the estimated target mass (OJ 287), $m_{t}$ we find that

$m_{t}=f_{G}\frac{m_{r}}{s_{N}}\frac{d_{t}^{2}}{d_{r}^{2}}$ (6)

or

$m_{t}\sim 2\times
10^{8}\frac{f_{G}}{0.53}\frac{m_{r}}{12}\frac{(d_{t}/1.073\times
10^{6})^{2}}{(d_{r}/5.75\times 10^{3})^{2}}(s_{N}/0.1\times 10^{4})^{-1}$ (7)

where we use values of $d_{t}=1.073$ Gpc and $d_{r}=5.75$ kpc (see Table 2).

Also using the reference source XTE J1859-226 we estimate

$m_{t}\sim 1.7\times 10^{8}$ (8)

applying Eq. (6) if we substitute the values of $f_{G}=\cos i_{r}/\cos
i_{t}=0.53$, $s_{N}=0.1\times 10^{4}$, $d_{r}=4.2$ kpc and $m_{r}=5.4$ (see
Table 2).

##  5 Discussion

Valtonen et al. (2012) suggested that a BH mass of the secondary in OJ 287 is
about $\sim 1.4\times 10^{8}M_{\odot}$ using the orbital solution. Komossa et
al. (2021), hereafter KOM21 found two X-ray spectral states in OJ 287 (see
Fig.4 there). Their so called ’hard state’ dominated by the power law
component and ’soft state’ characterized by a soft component which became very
strong when the source was very bright similar to those presented in Fig. 1
here. KOM21 also estimated the secondary BH mass, $M_{sec}$ as about
$1.5\times 10^{8}$ if $L_{x}/L_{\rm Edd}<1.7\times 10^{-2}$. On the other hand
KOM21 claimed that the Eddington ratio $L_{x}/L_{\rm Edd}$, that they found,
surprisingly low, $5.6\times 10^{-4}$ with an assumption that $M_{BH}\sim
1.8\times 10^{10}$ solar masses. In fact, our result confirmed their guess
that $L_{x}/L_{\rm Edd}$ is on the of order of $2\times 10^{-2}$ and thus, a
BH mass is about $2\times 10^{8}$ solar masses and is related to the secondary
BH in the OJ 287.

Our analysis of the Fall 2016 – Spring 2017 Swift observations identified two
distinct spectral states of OJ 287: the intermediate and a soft one. Those
states correspond to the initial rise of the X-ray spectral index correlated
to increase of the mass accretion rate (intermediate state) followed by
saturation of the spectral index after reaching a certain level of accretion
rate. The results of spectral analysis of the intermediate and soft states
were used to estimate the value of the mass accretion rate corresponding to
the saturation of the spectral index and consequently constrain the BH mass.
The third spectral state of OJ 287 that can be characterized by harder source
spectrum and represents a different accretion regime, was not the subject of
our analysis.

##  6 Conclusions

The multi-wavelength outburst activity in OJ 287 with X-ray telescope (XRT) on
the Swift call into question whether the source contains one or two black
holes (BHs). It is very important to reveal the characteristics of this
binary. In the presented work we demonstrated that the OJ 287 X-ray spectra
underwent the state transition from the intermediate state (IS) to the
high/soft state (HSS) (see Fig. 1). We found that energy spectra in all
spectral states can be modeled using a product of the tbabs and a
Comptonization component, BMC (see XSPEC). Moreover, in OJ 287 we discovered
the correlation of the photon index $\Gamma$ with normalization, N
(proportional to the disk mass accretion rate $\dot{M}$, see Fig. 2) similar
to those established in BH Galactic sources by ST09. We found that $\Gamma$
increases monotonically with $\dot{M}$ from the intermediate state (IS) to the
high-soft state HSS, and then saturates at $\Gamma\sim$ 2.6. This can be
considered as observational evidence of the presence of a BH in OJ 287. Based
on this correlation, we apply the scaling method (ST09) to estimate a BH mass
$\sim 2\times 10^{8}$ solar masses, using the well-studied X-ray BH binaries,
GX 339–4 and XTE J1859–226 as reference sources (see Fig. 3).

## Acknowledgements

This work made use of data supplied by the UK Swift Science Data Centre at the
University of Leicester. The authors are grateful to Dr. Kim Page of the UK
Swift Science Data Centre for the help with Swift/XRT data processing and
analysis.

The authors would like to thank Sergey Sazonov, Alexey Vikhlinin and Sergey
Trudolyubov for careful reading of the manuscript and valuable suggestions. We
are especially grateful to Colm Elliott, Marla Kouri-Towe and Michele Delmaire
for the proof reading and useful suggestions about corrections to the text of
the manuscript.

## Data availability

The Swift data is available to download through the UK Swift Data Science
website https://www.swift.ac.uk/archive.

## References

  * Borozdin et al. (1999) Borozdin, K., Revnivtsev, M., Trudolyubov, S. et al., 1999, ApJ, 517, 367
  * Corral-Santana et al. (2011) Corral-Santana, J.M. et al. 2011 MNRAS, 413, L15
  * Gupta et al. (2012) Gupta S. P., et al., 2012, New Astron., 17, 8
  * Heida et al. (2017) Heida, J., et al. 2017, ApJ, 846, 132
  * Hynes et al. (2002) Hynes, R.I. et al. 2002, MNRAS, 331, 189
  * Hynes et al. (2003) Hynes, R.I. et al. 2003, ApJ, 583, L95
  * Komossa et al. (2021) Komossa, S., et al., 2021, MNRAS, 504, 5575
  * Laurent & Titarchuk (2011) Laurent, P., & Titarchuk, L. 2011, ApJ, 727, 34 (LT11)
  * Laurent & Titarchuk (1999) Laurent, P., & Titarchuk, L. 1999, ApJ, 511, 289 (LT99)
  * Prince et al. (2021) Prince R., et al. 2021, MNRAS, 508, 315
  * Sagar et al. (2004) Sagar R., Stalin C. S., Gopal-Krishna, Wiita P. J., 2004, MNRAS, 348, 176
  * Seifina & Titarchuk (2011) Seifina, E., & Titarchuk, L. 2011, ApJ, 737, 128
  * Shakura & Sunyaev (1973) Shakura, N.I., & Sunyaev, R.A. 1973, A&A, 24, 337 (SS73)
  * Shaposhnikov & Titarchuk (2007) Shaposhnikov, N., & Titarchuk, L. 2007, ApJ, 663, 449 (ST07)
  * Shaposhnikov & Titarchuk (2009) Shaposhnikov, N., & Titarchuk, L. 2009, ApJ, 699, 453 (ST09)
  * Sunyaev & Titarchuk (1980) Sunyaev, R.A. & Titarchuk, L.G. 1980, A&A, 86, 121
  * Titarchuk et al. (1997) Titarchuk, L.G., Mastichiadis, A. & Kylafis, N.D. 1997, ApJ, 487, 834
  * Titarchuk & Zannias (1998) Titarchuk, L. & Zannias, T. 1998, ApJ, 499, 315
  * Titarchuk & Fiorito (2004) Titarchuk, L.G. & Fiorito, R. 2004, ApJ, 612, 988
  * Titarchuk, & Seifina (2021) Titarchuk, L. & Seifina, E. 2021, MNRAS, 501, 5659
  * Valtonen et al. (2006) Valtonen M. J. et al., 2006, ApJ, 646, 36
  * Valtonen et al. (2010) Valtonen M. J. et al., 2010, Celest. Mech. Dyn. Astron., 106, 235
  * Valtonen et al. (2012) Valtonen, M., Ciprini, S. & Pihajoki, P. 2012, MNRAS, 427, 77 (VAL12)
  * Willingale et al. (2013) Willingale, R., et al., 2013, MNRAS, 431, 394
  * Wilms et al. (2000) Wilms, J., Allen, A., McCray, R. 2000, ApJ, 542, 914