# Isovector spin susceptibility: Isotopic evolution of collectivity in spin
response

Kenichi Yoshida<EMAIL_ADDRESS>Department of Physics,
Kyoto University, Kyoto, 606-8502, Japan

###### Abstract

Background

Response to spin-dependent operators has been investigated in the magnetic
dipole and Gamow–Teller transitions, that provides magnetic properties of a
nuclear system.

Purpose

I investigate an isotopic dependence of the collectivity generated by the
spin-dependent interactions in the Ca and Ni isotopes through the isovector
(IV) spin-flip excitations. The responses in the neutral $(t_{z})$ and charge-
exchange $(t_{\pm})$ channels are considered in a unified way.

Method

A nuclear energy-density functional approach is employed for calculating the
response functions based on the Skyrme–Kohn–Sham–Bogoliubov method and the
quasiparticle-random-phase approximation (QRPA). I adopt the like-particle
QRPA and the proton–neutron QRPA for the neutral and charge-exchange channels,
respectively. I consider the fluctuation of the proton–neutron pair fields.

Results

The collective shift due to RPA correlations for the response in the neutral
channel is explained by the occupation probability of neutrons in the
$j_{>}=\ell+1/2$ orbital. Many particle–hole or two-quasiparticle excitations
have a coherent contribution to form a giant resonance in neutron-rich nuclei
for the charge-exchange channel. The IV spin susceptibility displays the
isotopic evolution of the collectivity and the underlying shell structure.

Conclusions

A repulsive character of the residual interaction in the spin–isospin channel
diminishes the IV spin susceptibility due to the collectivity, while the
dynamic ${}^{3}S$ pairing appearing in the charge-exchange channel opposes the
reduction.

††preprint: KUNS-2866

##  I Introduction

Susceptibility to an applied magnetic field is a fundamental property
characterizing the matter. A detailed investigation unveils the single-
particle and collective behavior of the spin degree of freedom. An interplay
of superconductivity and magnetism has attracted much interest and has been a
central issue in strongly-correlated many-body systems [1]. The fermionic
superfluidity or superconductivity is well understood by the ${}^{1}S_{0}$
pair correlation mostly [2, 3]. The onset of the spin susceptibility below the
critical temperature is a manifestation of the spin-triplet $(S=1)$
superfluidity of 3He [4, 5]. The $S=1$ pairing in electronic systems has been
investigated not only in the superfluid 3He [6] but in other systems recently
[7, 8]. In contrast, there has been a continuing discussion on the $S=1$
pairing in nuclear systems: the ${}^{3}S_{1}$ correlation of isoscalar (IS)
proton–neutron (pn) pairs in $N\sim Z$ nuclei [9] and the ${}^{3}P_{2}$
pairing of neutrons in the neutron-star matter [10].

Response of a nucleus to an external field displays the correlations among
constituent nucleons. The nuclear response is characterized by the transferred
angular momentum $\Delta L$, spin $\Delta S$, isospin $\Delta T$, and nucleon
number $\Delta N$ [11]. The response categorized as $\Delta L=0,\Delta
S=0,\Delta T=1,\Delta N=2$ pins down the ${}^{1}S_{0}$ pairing [3], and the
response $\Delta L=0,\Delta S=1,\Delta T=0,\Delta N=2$ has been studied
recently to uncover the ${}^{3}S_{1}$ pairing [12, 13, 14]. The response to
the spin-dependent electromagnetic and weak probes reveals the magnetic
properties of nuclear matter [15, 16, 17]. The IV spin magnetic-dipole (M1)
response with $\Delta L=0,\Delta S=1,\Delta T=1,\Delta N=0$ has been
investigated extensibly because they impact diverse quantities such as the
neutral-current neutrino-nucleus cross sections relevant to supernova physics
and neutrino physics [18] and the neutron-capture cross sections relevant to
the nucleosynthesis [19, 20, 21]. Furthermore, the IV spin-M1 response can be
seen in a wider perspective when it is considered as a single component
$\Delta T_{z}=0$ of the IV modes [22]. The additional components $\Delta
T_{z}=\pm 1$ represent the charge-exchange modes, namely, the Gamow–Teller
(GT) response.

The GT response can be a probe to disclose the ${}^{3}S_{1}$ IS pairing. The
dynamic IS pairing lowers the GT states in energy and thus shortens the
$\beta$-decay half-lives of neutron-rich nuclei [23] including deformed nuclei
[24]. The effects of the IS pairing on the GT transition strengths have been
studied in $N=Z$ odd-odd nuclei with a three-body model of two nucleons around
a spherical core [25]. A remarkable feature found in $N=Z$ odd-odd nuclei with
an $LS$-closed core (4He, 16O, 40Ca) is the appearance of the low-energy state
with a strong GT strength. A similar trait is also found by employing a
microscopic nuclear energy-density functional (EDF) method [26]. The low-
energy GT states have been indeed identified experimentally in the transitions
of ${}^{18}\text{O}\to{}^{18}\text{F}$ [27] and
${}^{42}\text{Ca}\to{}^{42}\text{Sc}$ [28, 29].

In this article, I am going to investigate the collectivity generated by the
spin response in view of the pair correlations. An isotopic dependence is
studied to see the shell effects. To this end, the EDF method is employed: a
theoretical model capable of handling nuclides with an arbitrary mass number
in a single framework [30, 31]. I introduce and evaluate the IV spin
susceptibility to discuss the neutron number dependence of the magnetic
property quantitatively.

This paper is organized in the following way: the theoretical framework for
describing the spin responses is given in Sec. II and details of the numerical
calculation are also given; Sec. III is devoted to the numerical results and
discussion based on the model calculation; the discussion on the response of
the Ca isotopes in the neutral and charge-exchange channels is given in Sec.
III.1 and Sec. III.2, respectively; the discussion for the Ni isotopes is
given in Sec. III.3; then, a summary is given in Sec. IV.

##  II Framework

###  II.1 KSB and QRPA for neutron-rich nuclei

Since the details of the formalism can be found in Refs. [32, 33, 24], here I
briefly recapitulate the basic equations relevant to the present study. In the
framework of the nuclear EDF method I employ, the ground state of a mother
(target) nucleus is described by solving the Kohn–Sham–Bogoliubov (KSB)
equation [34]:

$\displaystyle\begin{bmatrix}h^{q}(\bm{r}s)-\lambda^{q}&\tilde{h}^{q}(\bm{r}s)\\\
\tilde{h}^{q}(\bm{r}s)&-h^{q}(\bm{r}s)+\lambda^{q}\end{bmatrix}\begin{bmatrix}\varphi^{q}_{1,\alpha}(\bm{r}s)\\\
\varphi^{q}_{2,\alpha}(\bm{r}s)\end{bmatrix}=E_{\alpha}\begin{bmatrix}\varphi^{q}_{1,\alpha}(\bm{r}s)\\\
\varphi^{q}_{2,\alpha}(\bm{r}s)\end{bmatrix},$ (1)

where the KS potentials $h$ and $\tilde{h}$ are given by the functional
derivative of the EDF with respect to the particle density and the pair
density, respectively. An explicit expression of the potentials is found in
the Appendix of Ref. [35]. The superscript $q$ denotes $\nu$ (neutron,
$t_{z}=1/2$) or $\pi$ (proton, $t_{z}=-1/2$). The average particle number is
fixed at the desired value by adjusting the chemical potential $\lambda^{q}$.
When neutrons and/or protons are unpaired, the chemical potential is adapted
to lie in the center of the highest occupied and the lowest unoccupied
orbitals. Assuming the system is axially symmetric, the KSB equation (1) is
block diagonalized according to the quantum number $\Omega$, the $z$-component
of the angular momentum.

The excited states $|i\rangle$ are described as one-phonon excitations built
on the ground state $|0\rangle$ of the mother nucleus as

$\displaystyle|i\rangle$ $\displaystyle=\hat{\Gamma}^{\dagger}_{i}|0\rangle,$
(2) $\displaystyle\hat{\Gamma}^{\dagger}_{i}$
$\displaystyle=\sum_{\alpha\beta}\left\\{X_{\alpha\beta}^{i}\hat{a}^{\dagger}_{\alpha}\hat{a}^{\dagger}_{\beta}-Y_{\alpha\beta}^{i}\hat{a}_{\bar{\beta}}\hat{a}_{\bar{\alpha}}\right\\},$
(3)

where $\hat{a}^{\dagger}$ and $\hat{a}$ are the quasiparticle (qp) creation
and annihilation operators that are defined in terms of the solutions of the
KSB equation (1) with the Bogoliubov transformation. The phonon states, the
amplitudes $X^{i},Y^{i}$ and the vibrational frequency $\omega_{i}$, are
obtained in the quasiparticle-random-phase approximation (QRPA). The EDF gives
the residual interactions entering into the QRPA equation. For the axially
symmetric nuclei, the QRPA equation is block diagonalized according to the
quantum number $K=\Omega_{\alpha}+\Omega_{\beta}$. The deformed QRPA developed
in Refs. [32, 33] is extended to describe the response to the spin-dependent
operators.

###  II.2 Numerical procedures

To describe the developed neutron skin and the neutrons pair correlation
coupled with the continuum states that emerge uniquely in neutron-rich nuclei,
I solve the KSB equation in the coordinate space using cylindrical coordinates
$\bm{r}=(\rho,z,\phi)$ with a mesh size of $\Delta\rho=\Delta z=0.6$ fm and a
box boundary condition at $(\rho_{\mathrm{max}},z_{\mathrm{max}})=(14.7,14.4)$
fm. Since I assume further the reflection symmetry, only the region of $z\geq
0$ is considered. The qp states are truncated according to the qp energy
cutoff at 70 MeV, and the qp states up to the magnetic quantum number
$\Omega=23/2$ with positive and negative parities are included. I introduce
the truncation for the two-quasiparticle (2qp) configurations in the QRPA
calculations, in terms of the 2qp-energy as 60 MeV.

For the normal (particle–hole, p–h) part of the EDF, I employ the SGII
functional [36]. For the pairing energy, I adopt the one in Ref. [37] that
depends on both the IS and IV densities, in addition to the pair density, with
the parameters given in Table III of Ref. [37]. The same pairing EDF is
employed for the ${}^{1}S$ pn-pairing in the proton–neutron QRPA (pnQRPA)
calculation, while the linear term in the IV density is dropped. The effect of
the dynamic IS ($T=0$) paring is investigated by changing the strength of the
interaction. Note that the pnQRPA calculations including the dynamic IS
pairing with more or less the same strength as the ${}^{1}S$ pairing describe
well the characteristic low-lying GT strength distributions in the light
$N\simeq Z$ nuclei [28, 29, 27].

##  III Results and discussion

I consider the response of the Ca isotopes with the mass number
$A=42\text{--}70$ and the Ni isotopes with $A=52\text{--}78$ to the IV spin-
flip operator defined by

$\vec{F}_{K}=\mu\dfrac{1}{\sqrt{2}}\sum_{ss^{\prime}}\sum_{tt^{\prime}}\int\mathrm{d}\bm{r}\psi^{\dagger}(\bm{r}s^{\prime}t^{\prime})\psi(\bm{r}st)\langle
s^{\prime}|\sigma_{K}|s\rangle\langle t^{\prime}|\vec{\tau}|t\rangle,$ (4)

where $\mu=\mu_{\rm N}(\mathrm{g}_{\nu}-\mathrm{g}_{\pi})$, with $\mu_{\rm
N}=e\hbar/2mc$ being the nuclear magneton, $\mathrm{g}_{\nu}$,
$\mathrm{g}_{\pi}$ the spin g factor of a neutron and a proton, and
$\sigma_{K}$ and $\vec{\tau}=(\tau_{+1},\tau_{0},\tau_{-1})$ denote the
spherical components of the Pauli matrix of spin and isospin. With this
operator, the spin modes of the neutral (isospin unchanging: $\tau_{0}$) and
isospin changing ($\tau_{\pm 1}$) excitations are investigated in a unified
way. Note that I use here a probing operator that is $\mu\sqrt{2\pi}$ times
the one used in Ref. [22]. The IV spin-M1 and GT operators are indeed related
to the IV spin-flip operator (4) as
$F^{\textrm{IV,s}}_{K}(\textrm{M1})=\sqrt{\frac{3}{8\pi}}F_{K,0}$ and
$\vec{F}_{K}(\textrm{GT})=\frac{g_{A}}{\mu}\vec{F}_{K}$ when the contribution
of the meson-exchange and isobar currents is discarded [38]. Here $g_{A}$ is
the axial-vector coupling constant.

The Ca and Ni isotopes are all spherical in the ground state within the
present model. The neutrons are paired in all the isotopes except 70Ca, 56Ni,
and 78Ni while the protons are unpaired except in 52Ni.

###  III.1 Neutral channel: IV spin-M1 states in Ca isotopes

According to the shell model, in a spin-saturated closed-shell nucleus, such
as 40Ca, the M1 excitations are forbidden to occur. Thus the spin-flip
excitations in the Ca isotopes that I investigate here are by virtue of the
excitation of neutrons. It is noticed, however, that the M1 strengths in 40Ca
were observed around 10 MeV [39] and the beyond-RPA correlation is
indispensable to describe such a forbidden state [40].

Figure 1:  Calculated distributions (shifted) of the IV spin-flip transition
strengths in the neutral channel as functions of the excitation energy by
employing the SGII functional. The smearing parameter $\gamma=1$ MeV is used.
The results obtained without the RPA correlations (unp.) are depicted by the
dotted lines.  Figure 2:  Moments $m_{k}$ ($k=-1,0$ and $1$) of the transition
strengths in the Ca isotopes as functions of the neutron number. The results
obtained without the RPA correlations (unp.) are depicted by the dotted lines
with open symbols.

Figure 1 shows the calculated transition strengths for the IV spin-flip
excitation in the neutral channel, corresponding to the IV spin-M1 excitation:

$\displaystyle S(E)$
$\displaystyle=\sum_{K=-1}^{1}\dfrac{\mathrm{d}B(E,F_{K})}{\mathrm{d}E},$ (5)
$\displaystyle\dfrac{\mathrm{d}B(E,F_{K})}{\mathrm{d}E}$
$\displaystyle=\dfrac{2E\gamma}{\pi}\sum_{i}\dfrac{\tilde{E}_{i}|\langle
i|F_{K,0}|0\rangle|^{2}}{(E^{2}-\tilde{E}_{i}^{2})^{2}+E^{2}\gamma^{2}},$ (6)

where $\tilde{E}_{i}^{2}=(\hbar\omega_{i})^{2}+\gamma^{2}/4$ [41]. The
smearing width $\gamma$ is set to 1 MeV, which is supposed to simulate the
spreading effect, $\Gamma^{\downarrow}$, missing in the QRPA.

In the ${}^{42\text{--}48}$Ca isotopes, the spin-M1 resonant state is
generated by the $\nu f_{7/2}\to\nu f_{5/2}$ excitation. The unperturbed state
is shifted up in energy due to the residual interaction of the
$\sigma_{1}\cdot\sigma_{2}\tau_{1}\cdot\tau_{2}$ type. To discuss the isotopic
dependence of the strength distributions quantitatively, I show in Fig. 2 the
moment of the strengths defined by

$m_{k}=\sum_{K}\sum_{i}(\hbar\omega_{i})^{k}|\langle i|F_{K,0}|0\rangle|^{2}$
(7)

for $k=-1,0$, and 1. A linear rise seen in $m_{0}$ from 42Ca to 48Ca is
understood by the occupation of neutrons in the $f_{7/2}$ orbital. A p–h-type
transition strength is proportional to $u^{2}v^{2}$ of the
Bardeen–Cooper–Schrieffer (BCS) amplitude. Since the $f_{5/2}$ orbital is
completely empty $u_{f_{5/2}}=1$, the transition strengths develop
monotonically as $v_{f_{7/2}}^{2}$ increases. Here, the occupation probability
$v^{2}_{\alpha}$ can be evaluated by
$\int\mathrm{d}\bm{r}\sum_{s}|\varphi_{2,\alpha}(\bm{r}s)|^{2}$ using the qp
wave functions. Comparing the results with and those without the RPA
correlations, the effect of RPA correlations is also enhanced from 42Ca to
48Ca. The mean excitation energy evaluated by $m_{1}/m_{0}$ with (without) the
RPA correlations is 13.8 (13.2), 12.5 (11.5), 11.6 (10.1), and 10.9 (8.7) MeV
for 42,44,46,48Ca, respectively. This is again due to a gradual occupation of
neutrons in the $f_{7/2}$ orbital. A two-body matrix element entering the QRPA
equation is approximately proportional to $uvuv$ of the BCS amplitude. In the
present case, only the $\nu f_{7/2}\to\nu f_{5/2}$ excitation accounts for the
appearance of spin-M1 resonance in essence. Thus, the energy shift is due to
the diagonal matrix element, which is proportional to
$u^{2}_{f_{5/2}}v^{2}_{f_{7/2}}$ and evidently the occupation number of
neutrons in the $f_{7/2}$ orbital.

One sees an onset of the low-energy state in ${}^{50\text{--}56}$Ca. This is
due to the $\nu p_{3/2}\to\nu p_{1/2}$ excitation. The increase in $m_{0}$
from 50Ca to 52Ca is understood by the occupation of neutrons in the $p_{3/2}$
orbital similarly in the case of the $f_{7/2}$ orbital in
${}^{42\text{--}48}$Ca. When the neutron number increases further, neutrons
start to occupy the $p_{1/2}$ orbital. Then, the occupation of neutrons in the
$p_{1/2}$ orbital diminishes the transition strength as
$u^{2}_{p_{1/2}}v^{2}_{p_{3/2}}=1-v^{2}_{p_{1/2}}$. A linear reduction in
$m_{0}$ and the energy shift due to the RPA correlations is seen up to 60Ca,
where the occupation of neutrons in the $f_{5/2}$ orbital is fulfilled. Note
that the $f_{5/2}$ orbital is bound though the separation energy
$E_{f_{5/2}}-\lambda^{\nu}$ is as low as 1.4 MeV. In 60Ca, the spin excitation
is forbidden. A tiny strength, however, appears in the calculation. The
appearance of the strengths is due to the pairing of neutrons. In the present
calculation, the occupation probability of the $\nu 1g_{9/2}$ orbital is 0.08.

Beyond 60Ca, a linear increase in $m_{0}$ and the energy shift due to the RPA
correlations is mainly due to a gradual occupation of neutrons in the
$1g_{9/2}$ orbital. A distinct feature in the neuron-rich nuclei beyond $N=40$
is that the $g_{7/2}$ orbital is no more a bound state. Furthermore, there can
be an appreciable contribution of 2qp excitations in the continuum because the
resonance shape is not simply a single peak as calculated below $N=40$.
Therefore, the spin-M1 state acquires a wider width due to a stronger coupling
to the continuum states. Since the continuum states are discretized in the
present calculation, a quantitative discussion on the resonance shape is not
easy. It is thus an interacting future work to extend the continuum
calculation like in Ref. [42] to heavier nuclei and discuss the neutron-number
dependence of the width of the spin-M1 resonance [43]. The evolution of the
spin-M1 strength distributions in the Ca isotopes is similar to the one
predicted by the relativistic EDF approach [44]. It is noticed that the
pairing governing the isotopic dependence of the RPA correlations in the
spin-M1 excitation is the static pairing effect. The dynamic pairing caused by
the residual pair interaction does not affect the collectivity of the
unnatural parity states in the present model, while it enhances the
collectivity of the natural parity states as discussed, e.g., in Refs. [45,
46, 47, 48].

The inverse energy-weighted sum $m_{-1}$ shown in Fig. 2 is related to the IV
spin susceptibility [49] as

$\chi_{z}\equiv 2\mu^{2}\sum_{i}\dfrac{|\langle
i|S_{z}T_{z}|0\rangle|^{2}}{\hbar\omega_{i}}=\dfrac{m_{-1}}{3},$ (8)

where
$S_{z}T_{z}=\int\mathrm{d}\bm{r}\sum\psi^{\dagger}(\bm{r}s^{\prime}t^{\prime})\psi(\bm{r}st)\langle
s^{\prime}|\frac{\sigma_{z}}{2}|s\rangle\langle t^{\prime}|\tau_{0}|t\rangle$.
This is of particular interest because it is not sensitive to the high-energy
side of the response and hence is less affected by the 2p2h excitations [15].
For nuclear matter, the IV magnetic susceptibility per nucleon is simply
expressed [50] as

$\bar{\chi}_{z}=\dfrac{1+\frac{1}{3}F_{1}}{1+G_{0}^{\prime}}\mu^{2}N_{0}$ (9)

with the Landau parameters $N_{0}$, $G_{0}^{\prime}$ and $F_{1}$, and
$\mu^{2}N_{0}$ is regarded as the Pauli paramagnetic susceptibility. The
susceptibility of the correlated systems is suppressed due to the repulsive
nature of the interaction $G_{0}^{\prime}>0$. Note that the unperturbed value
displayed in the top panel of Fig. 2 includes the effect of $F_{1}$, i.e., the
effective mass that reduces the susceptibility and the like-particle pairing
that also reduces the susceptibility as for the Belyaev inertia [51]. In the
present case for the Ca isotopes, since protons do not take part in generating
the collectivity, the spin-M1 state is simply described as a neutron p–h
excitation of the spin–orbit partners. Therefore, the reduction of the
susceptibility strongly depends on the shell structure and the pairing
correlation of neutrons. The effect of the residual interaction is seen more
clearly than in $m_{0}$.

###  III.2 Charge-exchange channel: GT states in Ca isotopes

Figure 3:  Similar to Fig. 1 but for the response to the $t^{-}$ channel of
the operator (4) in (a) ${}^{42\text{--}48}$Ca ($f_{7/2}$ shell), (b)
${}^{50\text{--}60}$Ca ($p\mathchar 45\relax f_{5/2}$ shell) and (c)
${}^{62\text{--}70}$Ca ($g_{9/2}$ shell). The results including the $T=0$
(spin-triplet) pairing are drawn by the thin lines.

Figure 3 shows the calculated transition strengths for the IV spin-flip
excitation in the charge-exchange channel ($t_{-}$), corresponding to the GT
excitation, as given similarly to Eqs. (5), (6) except by replacing the
excitation energy $E$ with $E_{T}=E-\lambda^{\nu}+\lambda^{\pi}$ so that the
excitation is with respect to the ground state of the mother (target) nucleus.
The transition in the $t_{+}$ channel is strongly suppressed as the excitation
of protons in the neutral channel.

Let me discuss the GT excitations in ${}^{42\text{--}48}$Ca shown in Fig.
3(a). One sees a two-hump structure while there is only a single peak in the
spin-M1 response. These states are mainly constructed by the $\nu
f_{7/2}\to\pi f_{7/2}$ and $\nu f_{7/2}\to\pi f_{5/2}$ configurations. Since
the $(\nu f_{7/2})^{2}$ configuration couples to $S=0$, this does not take
part in the formation of the spin-M1 state. With the growth of the neutron
number, the transition strength increases. This is understood by a gradual
increase in the occupation of neutrons in the $f_{7/2}$ orbital as discussed
for a linear rise in the transition strengths in the neutral channel. The
excitation energy is predominantly given by the unperturbed energy and the
diagonal matrix element of the
$\sigma_{1}\cdot\sigma_{2}\tau_{1}\cdot\tau_{2}$ interaction, and the two
peaks are shifted up in energy due to a repulsive character of the
interaction. The energy shift is also nearly proportional to the occupied
neutron number in the $f_{7/2}$ orbital.

As neutrons occupy the $2p$ orbitals, the $\nu p_{3/2,1/2}\to\pi p_{3/2,1/2}$
excitations show up in low energies as depicted by the dotted lines in Fig.
3(b). Thanks to the RPA correlations, most of the strengths are gathered into
the giant resonance though there remains a tiny strength in the low energy
region. These low-lying states are responsible for the $\beta$ decay.

The IS pairing affects the GT states [23, 24] and explains well the
enhancement of the strength for the low-energy GT state in 18O (42Ca) [28, 29,
27, 26, 52], where the particle–particle (p–p) type excitations of $\nu
p_{3/2}\to\pi p_{3/2}$ ($\nu f_{7/2}\to\pi f_{7/2}$) and $\nu p_{3/2}\to\pi
p_{1/2}$ ($\nu f_{7/2}\to\pi f_{5/2}$) cooperatively participate in generating
the low-energy GT state. The numerical results including the dynamic IS
pairing are depicted by the thin line in Fig. 3. The effect of the IS pairing
is strong in 42Ca as discussed in Ref. [26]. Here, the strength of the IS pair
interaction was set as the same as that for the IV pair interaction.

Beyond $N=34$, one sees an appreciable amount of the unperturbed strengths is
developed in low energies with an increase in the neutron number. This comes
from the $\nu f_{5/2}\to\pi f_{7/2}$ excitation. Due to the RPA correlations,
however, most of the strengths are brought together into the high-lying giant
resonance. The remaining low-lying state is expected to be affected by the IS
pairing because of the p–p nature of the configuration. In 54,56Ca, the
occupation probability of neutrons in the $f_{5/2}$ orbital is 0.2 and 0.4. In
the present cases, the low-lying state is lowered in energy and enhanced in
strengths only slightly , as shown by the thick line in Fig. 3(b). This is
because the $\nu f_{7/2}\to\pi f_{5/2}$ excitation is purely a p–h type and is
not affected by the IS pairing.

The IS pairing affects the giant resonance , as neutrons partially occupy the
$g_{9/2}$ orbital. Beyond $N=40$, the $\nu g_{9/2}\to\pi g_{9/2}$
configuration appears more clearly in the high energy region according to a
gradual increase in the occupation probability of neutrons in the $g_{9/2}$
orbital. Since this configuration is a p–p type excitation, the IS pairing is
active. As one can see in Fig. 3(c), the giant resonance is lowered in energy
due to the attractive nature of the IS pairing.

In neutron-rich isotopes, the number of the available p–h or 2qp excitations
increases because of the imbalanced Fermi level of neutrons and protons.
Therefore, the GT strengths are concentrated in the giant resonance in
${}^{62\text{--}70}$Ca as a collective effect notwithstanding that a
significant fraction of the GT strengths is found inside the $Q_{\beta}$
window without the RPA correlations. The super-allowed GT resonance predicted
in the light neutron drip-line nuclei [53] is thus unlikely to occur in the Ca
isotopes even near the drip line.

Figure 4:  Top: IV spin susceptibility $\chi_{z}$ of the Ca isotopes obtained
with and without the RPA correlations. Middle: IV spin susceptibility
$\chi_{\perp}$ in the charge-exchange channel. Bottom: Ratio of the IV-spin
susceptibilities $\chi_{\perp}/\chi_{z}$ obtained with and without the dynamic
$T=0$ (spin-triplet) pairing.

To investigate the magnetic property systematically, I introduce the static
susceptibility. The IV spin susceptibility in the charge-exchange channel is
given as

$\chi_{\perp}=\mu^{2}\left[\sum_{i}\dfrac{|\langle
i|S_{z}T_{-}|0\rangle|^{2}}{\hbar\omega_{i}}+\sum_{i}\dfrac{|\langle
i|S_{z}T_{+}|0\rangle|^{2}}{\hbar\omega_{i}}\right],$ (10)

where
$S_{z}T_{\pm}=\int\mathrm{d}\bm{r}\sum\psi^{\dagger}(\bm{r}s^{\prime}t^{\prime})\psi(\bm{r}st)\langle
s^{\prime}|\frac{\sigma_{z}}{2}|s\rangle\langle t^{\prime}|\tau_{\pm
1}|t\rangle$. For nuclear matter, one may expect $\chi_{\perp}$ coincides with
$\chi_{z}$. The IV spin susceptibility in the charge-exchange channel is
reduced due to the RPA correlations as in the neutral channel.

I show in the middle panel of Fig. 4 the calculated IV spin susceptibilities
$\chi_{\perp}$ obtained with and without the RPA correlations. As the neutron
number increases, the p–h or 2qp excitations possessing the non-vanishing GT
matrix-element appear in low energies. The susceptibility in the free system
thus increases. Since the operator (4) changes only the direction of spin and
isospin, and does not change the spatial structure, the p–h or 2qp
configurations in the $0\hbar\omega_{0}$ excitation are only possible to
appear in the GT response. Around $N=40$, the $-1\hbar\omega_{0}$ excitations
start to appear in low energies [54], and the $0\hbar\omega_{0}$ excitations
show up in a relatively higher energy region. Therefore, $\chi_{\perp}$ keeps
almost unchanged beyond $N=40$. The effect of the RPA correlations becomes
apparent as the collectivity becomes strong. With an increase in the neutron
number, the number of p–h or 2qp excitations in the $pf$-shell increases, and
they make a coherent contribution to the formation of the giant resonance. One
can thus see a systematic reduction of $\chi_{\perp}$ from the unperturbed one
in the range of $22\leq N\leq 40$.

The dynamic IS pairing acts on the spin response oppositely to the p–h
residual interaction, the former causes an attractive shift of the transitions
and enhances the susceptibility while the latter a repulsive shift and reduces
the susceptibility. For nuclear matter one may assume to express the ratio of
the IV spin susceptibilities as

$\dfrac{\chi_{\perp}}{\chi_{z}}=\dfrac{1+G_{0}^{\prime}}{1+G_{0}^{\prime}+V_{\textrm{pair}}}=1-\dfrac{V_{\textrm{pair}}}{1+G_{0}^{\prime}+V_{\textrm{pair}}}$
(11)

with the IS pair interaction $V_{\textrm{pair}}<0$. For asymmetric systems,
there should be a correction due to the neutron excess. A recently developed
technique [55] can be used to evaluate this, but it is beyond the scope of the
present work. It is also noticed that the dynamic IS pairing affects the p–h
type excitation only through the static IV pairing. Within the scope of
consideration, roles of the IS pairing can thus be investigated only at a
minimum. In plotting the ratio $\chi_{\perp}/\chi_{z}$, nevertheless, one can
expect to see the magnetic properties from different perspectives. Shown in
the bottom panel of Fig. 4 is the ratio for the Ca isotopes. Since protons do
not take part in the spin excitations as mentioned in the beginning, one sees
a high asymmetry between $\chi_{z}$ and $\chi_{\perp}$ even in the case of
$N\sim Z$. The enhancement in the ratio $\chi_{\perp}/\chi_{z}$ is governed by
the shell effect in $\chi_{z}$, and displayed is only a tiny contribution of
the dynamic IS pairing to the susceptibility.

Figure 5:  Same as Fig. 4 but for the Ni isotopes.

###  III.3 Ni isotopes

Figure 6:  Same as Fig. 1 but for the Ni isotopes.

An outline of the isotopic dependence of the magnetic properties can be
illustrated in the IV spin susceptibilities, $\chi_{z}$ and $\chi_{\perp}$.
Thus I briefly discuss the isotopic evolution of the collectivity in the spin
responses for the Ni isotopes by showing Fig. 5.

Let me discuss first $\chi_{z}$. In the present case, the $\pi f_{7/2}\to\pi
f_{5/2}$ excitation participates in the spin-M1 states in addition to the
excitation of neutrons. As in the case for the Ca isotopes, the isotopic
evolution of the collectivity is mainly explained by the occupation of
neutrons in the $j_{>}$ orbital. To fortify a conclusion, I show in Fig. 6 the
strength distributions of the IV spin-M1 excitation. A gradual increase in the
collectivity up to $N=28$, enhancement in the transition strength in the giant
resonance region and suppression of the susceptibility due to RPA
correlations, is because of the occupation of neutrons in the $f_{7/2}$
orbital. Taking a close look at the strength distribution, one sees that two
peaks appear in the giant resonance region in 52Ni, while two peaks merge into
a single peak in 56Ni. The lower-energy and higher-energy states are mainly
generated by the $\pi f_{7/2}\to\pi f_{5/2}$ excitation and the $\nu
f_{7/2}\to\nu f_{5/2}$ excitation, respectively. The energy difference for
these p–h excitations is about 2 MeV in 52Ni, while about 0.2 MeV in 56Ni.
Thus, the coherence between these excitations develops from 52Ni to 56Ni. When
neutrons occupy the $p_{3/2}$ orbital, a low-lying state shows up. The
appearance of the low-lying state has a contribution to the increase in
$\chi_{z}$ in ${}^{50\text{--}52}$Ca. In ${}^{58\text{--}60}$Ni, however, the
mixing between the $\pi f_{7/2}\to\pi f_{5/2}$ and $\nu f_{7/2}\to\nu f_{5/2}$
excitations becomes weak, which leads to the cancellation of $\chi_{z}$. As
neutrons occupy the $p_{1/2}$ and $f_{5/2}$ orbitals, the collectivity
decreases further. At $N=40$, the spin excitation is forbidden for neutrons.
Beyond $N=40$, the $\nu g_{9/2}\to\nu g_{7/2}$ excitation starts to appear in
the high energy. The lower-energy peak around 8–9 MeV and the higher-energy
peak around 10–11 MeV are predominantly generated by the $\pi f_{7/2}\to\pi
f_{5/2}$ excitation and the $\nu g_{9/2}\to\nu g_{7/2}$ excitation,
respectively. As in the case for the Ca isotopes, the transition strengths in
high energies develop with an increase in the occupation of neutrons in the
$g_{9/2}$ orbital.

Figure 7:  Similar to Fig. 3 but for (a) ${}^{52\text{--}56}$Ni ($f_{7/2}$
shell), (b) ${}^{58\text{--}68}$Ni ($p\mathchar 45\relax f_{5/2}$ shell) and
(c) ${}^{70\text{--}78}$Ni ($g_{9/2}$ shell).

I then discuss $\chi_{\perp}$. In ${}^{52\text{--}56}$Ni, the excitation in
the $t_{+}$ channel has an appreciable contribution to the IV spin
susceptibility as well as that in the $t_{-}$ channel. Therefore,
$\chi_{\perp}$ value is higher than in the Ca isotopes below $N=28$. To see
what is happening more clearly, I show in Fig. 7(a) the transition strengths
in the $t_{+}$ channel by the dashed line, where $E$ is replaced by
$E_{T}=E+\lambda^{\nu}-\lambda^{\pi}$ in Eq. 6; the transition strengths in
the $t_{-}$ channel are drawn by the solid line. The low-lying state in the
$t_{+}$ channel is mainly generated by the $\pi f_{7/2}\to\nu f_{7/2}$
excitation. In 52Ni, the neutron occupation probability in the $f_{7/2}$
orbital is 0.48. This configuration is then considered as a hole–hole type
excitation and is thus affected by the dynamic IS pairing, though not so
strongly. The GT state in the $t_{-}$ channel is constructed by the $\nu
f_{7/2}\to\pi f_{5/2}$ excitation, which can be considered as a p–p type
excitation at $N=24$. Since the $\nu f_{7/2}\to\pi f_{7/2}$ excitation is
forbidden, the collectivity and the influence by the IS pairing are weaker
than in the Ca isotopes.

Beyond $N=28$, one sees in the middle panel of Fig. 5 that the collectivity
develops with an increase in the neutron number as in the Ca cases. At $N>34$,
the collectivity does not increase in the Ni isotopes, while one has a gradual
increase in the collectivity toward $N=40$ in the Ca isotopes. This different
isotopic dependence results from the fact that the $\nu f_{5/2}\to\pi f_{7/2}$
excitation does not take part in generating the collectivity due to the Pauli
effect in the Ni isotopes. Beyond $N=40$, the IV spin susceptibility keeps
almost constant similarly in the Ca isotopes because the strength distribution
in the low energy is almost unchanged and the $\nu g_{9/2}\to\pi g_{9/2}$
excitation appears in a high energy region. Note that the negative parity
states show up instead in the low energy, leading to an interplay between the
allowed and first-forbidden $\beta$-decays [56]. Since the $\nu g_{9/2}\to\pi
g_{9/2}$ excitation is a p–p type excitation around $N=40$, the IS pairing
affects the giant resonance as shown in Figs. 7(b) and 7(c).

Finally, let me mention the isotopic dependence of the ratio of the IV spin
susceptibilities, $\chi_{\perp}/\chi_{z}$, showing in the bottom panel of Fig.
5. In 56Ni, its value is unity, where protons and neutrons behave
symmetrically. Furthermore, the IS pairing is ineffective for the spin
susceptibility because of the shell closure. Above $N=28$ the collectivity
becomes stronger and stays almost constant beyond $N=34$ in the charge-
exchange channel, while that in the neutral channel decreases toward $N=40$.
Then, one sees a gradual increase in the ratio. Beyond $N=40$, the
collectivity is enhanced in the neutral channel, leading to a decrease in the
ratio as in the Ca isotopes. The effect of the IS pairing is much weaker than
in the case of the Ca isotopes.

##  IV Summary

A comparative study on the isotopic dependence of the magnetic properties has
been performed by investigating the IV spin-flip excitations in the Ca and Ni
isotopes. The responses in the neutral and charge-exchange channels were
considered in a unified way. I made use of the nuclear EDF method for
calculating the response functions based on the Skyrme–KSB and the QRPA. The
like-particle QRPA and the proton–neutron QRPA were employed for the neutral
and charge-exchange channels, respectively. The collective shift due to RPA
correlations for the response in the neutral channel is mainly explained by
the occupation probability of neutrons in the $j_{>}$ orbital, where the IV
pairing of neutrons gives a dominant role. In the charge-exchange channel,
many p–h or 2qp excitations have a coherent contribution to form a giant
resonance in neutron-rich nuclei. The dynamic IS pairing lowers the low-lying
and giant resonance states sensitively to the shell structure. I have found
that the isotopic evolution of the collectivity and the shell structure are
nicely displayed by the IV spin susceptibility. A repulsive character of the
residual interaction in the spin–isospin channel diminishes the
susceptibility, while the IS pairing appearing in the charge-exchange channel
opposes the suppression.

###### Acknowledgements.

This work was supported by the JSPS KAKENHI (Grants No. JP19K03824 and No.
JP19K03872). The numerical calculations were performed on Yukawa-21 at the
Yukawa Institute for Theoretical Physics, Kyoto University.

## References

  * Bulaevskii  _et al._ [1985] L. N. Bulaevskii, A. I. Buzdin, M. L. Kulić, and S. V. Panjukov, Coexistence of superconductivity and magnetism theoretical predictions and experimental results, Adv. Phys. 34, 175 (1985).
  * Schrieffer [1964] J. Schrieffer,  _Theory of Superconductivity_ , Frontiers in physics (W.A. Benjamin, 1964).
  * Brink and Broglia [2005] D. M. Brink and R. A. Broglia,  _Nuclear Superfluidity: Pairing in Finite Systems_, Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology (Cambridge University Press, 2005).
  * Balian and Werthamer [1963] R. Balian and N. R. Werthamer, Superconductivity with pairs in a relative $p$ wave, Phys. Rev. 131, 1553 (1963).
  * Leggett [1965] A. J. Leggett, Spin susceptibility of a superfluid fermi liquid, Phys. Rev. Lett. 14, 536 (1965).
  * Vollhardt and Wolfle [1990] D. Vollhardt and P. Wolfle,  _The superfluid phases of helium 3_ (Taylor and Francis Inc., 1990).
  * Joynt and Taillefer [2002] R. Joynt and L. Taillefer, The superconducting phases of $\mathrm{UPt}_{3}$, Rev. Mod. Phys. 74, 235 (2002).
  * Mackenzie and Maeno [2003] A. P. Mackenzie and Y. Maeno, The superconductivity of $\mathrm{Sr}_{2}\mathrm{RuO}_{4}$ and the physics of spin-triplet pairing, Rev. Mod. Phys. 75, 657 (2003).
  * Frauendorf and Macchiavelli [2014] S. Frauendorf and A. O. Macchiavelli, Overview of neutron–proton pairing, Prog. Part. Nucl. Phys. 78, 24 (2014), arXiv:1405.1652 [nucl-th] .
  * Haensel  _et al._ [2007] P. Haensel, A. Potekhin, and D. Yakovlev,  _Neutron Stars 1: Equation of State and Structure_ , Astrophysics and Space Science Library (Springer New York, 2007).
  * Harakeh and Woude [2001] M. Harakeh and A. Woude,  _Giant Resonances: Fundamental High-frequency Modes of Nuclear Excitation_, Oxford science publications (Oxford University Press, 2001).
  * Yoshida [2014] K. Yoshida, Proton-neutron pairing vibrations in $N=Z$ nuclei: Precursory soft mode of isoscalar pairing condensation, Phys. Rev. C 90, 031303 (2014), arXiv:1409.4884 [nucl-th] .
  * Litvinova  _et al._ [2018] E. Litvinova, C. Robin, and I. A. Egorova, Soft modes in the proton–neutron pairing channel as precursors of deuteron condensate in $N=Z$ nuclei, Phys. Lett. B 776, 72 (2018), arXiv:1612.09182 [nucl-th] .
  * Chazono  _et al._ [2021] Y. Chazono, K. Yoshida, K. Yoshida, and K. Ogata, Proton induced deuteron knockout reaction as a probe of an isoscalar proton-neutron pair in nuclei, Phys. Rev. C 103, 024609 (2021), arXiv:2007.06771 [nucl-th] .
  * Ericson [1984] M. Ericson, Nuclear spin/isospin response functions, Prog. Part. Nucl. Phys. 11, 277 (1984).
  * Osterfeld [1992] F. Osterfeld, Nuclear spin and isospin excitations, Rev. Mod. Phys. 64, 491 (1992).
  * Heyde  _et al._ [2010] K. Heyde, P. von Neumann-Cosel, and A. Richter, Magnetic dipole excitations in nuclei: elementary modes of nucleonic motion, Rev. Mod. Phys. 82, 2365 (2010), arXiv:1004.3429 [nucl-ex] .
  * Nakamura  _et al._ [2017] S. X. Nakamura, H. Kamano, Y. Hayato, M. Hirai, W. Horiuchi, S. Kumano, T. Murata, K. Saito, M. Sakuda, T. Sato, and Y. Suzuki, Towards a unified model of neutrino-nucleus reactions for neutrino oscillation experiments, Rept. Prog. Phys. 80, 056301 (2017), arXiv:1610.01464 [nucl-th] .
  * Käppeler  _et al._ [2011] F. Käppeler, R. Gallino, S. Bisterzo, and W. Aoki, The $s$ process: Nuclear physics, stellar models, and observations, Rev. Mod. Phys. 83, 157 (2011).
  * Arnould  _et al._ [2007] M. Arnould, S. Goriely, and K. Takahashi, The r-process of stellar nucleosynthesis: Astrophysics and nuclear physics achievements and mysteries, Phys. Rept. 450, 97 (2007), arXiv:0705.4512 [astro-ph] .
  * Kajino  _et al._ [2019] T. Kajino, W. Aoki, A. B. Balantekin, R. Diehl, M. A. Famiano, and G. J. Mathews, Current status of r -process nucleosynthesis, Prog. Part. Nucl. Phys. 107, 109 (2019), arXiv:1906.05002 [astro-ph.HE] .
  * Auerbach and Klein [1984] N. Auerbach and A. Klein, Structure of isovector spin excitations in nuclei, Phys. Rev. C 30, 1032 (1984).
  * Engel  _et al._ [1999] J. Engel, M. Bender, J. Dobaczewski, W. Nazarewicz, and R. Surman, Beta decay of r process waiting point nuclei in a selfconsistent approach, Phys. Rev. C 60, 014302 (1999), arXiv:nucl-th/9902059 .
  * Yoshida [2013] K. Yoshida, Spin–isospin response of deformed neutron-rich nuclei in a self-consistent Skyrme energy-density-functional approach, Prog. Theor. Exp. Phys. 2013, 113D02 (2013), [Erratum: Prog. Theor. Exp. Phys., 2021, 019201 (2021)], arXiv:1308.0424 [nucl-th] .
  * Tanimura  _et al._ [2014] Y. Tanimura, H. Sagawa, and K. Hagino, Three-body model calculations for $N=Z$ odd–odd nuclei with $T=0$ and $T=1$ pairing correlations, Prog. Theor. Exp. Phys. 2014, 053D02 (2014), arXiv:1309.3394 [nucl-th] .
  * Bai  _et al._ [2014] C. L. Bai, H. Sagawa, G. Colò, Y. Fujita, H. Q. Zhang, X. Z. Zhang, and F. R. Xu, Low-energy collective Gamow-Teller states and isoscalar pairing interaction, Phys. Rev. C 90, 054335 (2014).
  * Fujita  _et al._ [2019] H. Fujita, Y. Fujita, Y. Utsuno, K. Yoshida, T. Adachi, A. Algora, M. Csatlós, J. M. Deaven, E. Estevez-Aguado, C. J. Guess, J. Gulyás, K. Hatanaka, K. Hirota, R. Hutton, D. Ishikawa, A. Krasznahorkay, H. Matsubara, F. Molina, H. Okamura, H. J. Ong, G. Perdikakis, B. Rubio, C. Scholl, Y. Shimbara, G. Süsoy, T. Suzuki, A. Tamii, J. H. Thies, R. G. T. Zegers, and J. Zenihiro, Experimental study of Gamow-Teller transitions via the high-energy-resolution ${}^{18}\mathrm{O}(^{3}\mathrm{He},t)^{18}\mathrm{F}$ reaction: Identification of the low-energy “super” -Gamow-Teller state, Phys. Rev. C 100, 034618 (2019).
  * Fujita  _et al._ [2014] Y. Fujita, H. Fujita, T. Adachi, C. L. Bai, A. Algora, G. P. A. Berg, P. von Brentano, G. Colò, M. Csatlós, J. M. Deaven, E. Estevez-Aguado, C. Fransen, D. De Frenne, K. Fujita, E. Ganioğlu, C. J. Guess, J. Gulyás, K. Hatanaka, K. Hirota, M. Honma, D. Ishikawa, E. Jacobs, A. Krasznahorkay, H. Matsubara, K. Matsuyanagi, R. Meharchand, F. Molina, K. Muto, K. Nakanishi, A. Negret, H. Okamura, H. J. Ong, T. Otsuka, N. Pietralla, G. Perdikakis, L. Popescu, B. Rubio, H. Sagawa, P. Sarriguren, C. Scholl, Y. Shimbara, Y. Shimizu, G. Susoy, T. Suzuki, Y. Tameshige, A. Tamii, J. H. Thies, M. Uchida, T. Wakasa, M. Yosoi, R. G. T. Zegers, K. O. Zell, and J. Zenihiro, Observation of Low- and High-Energy Gamow-Teller Phonon Excitations in Nuclei, Phys. Rev. Lett. 112, 112502 (2014).
  * Fujita  _et al._ [2015] Y. Fujita, H. Fujita, T. Adachi, G. Susoy, A. Algora, C. L. Bai, G. Colò, M. Csatlós, J. M. Deaven, E. Estevez-Aguado, C. J. Guess, J. Gulyás, K. Hatanaka, K. Hirota, M. Honma, D. Ishikawa, A. Krasznahorkay, H. Matsubara, R. Meharchand, F. Molina, H. Nakada, H. Okamura, H. J. Ong, T. Otsuka, G. Perdikakis, B. Rubio, H. Sagawa, P. Sarriguren, C. Scholl, Y. Shimbara, E. J. Stephenson, T. Suzuki, A. Tamii, J. H. Thies, K. Yoshida, R. G. T. Zegers, and J. Zenihiro, High-resolution study of Gamow-Teller excitations in the ${}^{42}\mathrm{Ca}(^{3}\mathrm{He},t)\phantom{\rule{2.79999pt}{0.0pt}}^{42}\mathrm{Sc}$ reaction and the observation of a “low-energy super-Gamow-Teller state”, Phys. Rev. C 91, 064316 (2015).
  * Bender  _et al._ [2003] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Self-consistent mean-field models for nuclear structure, Rev. Mod. Phys. 75, 121 (2003).
  * Nakatsukasa  _et al._ [2016] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Yabana, Time-dependent density-functional description of nuclear dynamics, Rev. Mod. Phys. 88, 045004 (2016), arXiv:1606.04717 .
  * Yoshida and Van Giai [2008] K. Yoshida and N. Van Giai, Deformed quasiparticle-random-phase approximation for neutron-rich nuclei using the Skyrme energy density functional, Phys. Rev. C 78, 064316 (2008), arXiv:0809.0169 [nucl-th] .
  * Yoshida and Nakatsukasa [2013] K. Yoshida and T. Nakatsukasa, Shape evolution of giant resonances in Nd and Sm isotopes, Phys. Rev. C 88, 034309 (2013), arXiv:1305.6437 [nucl-th] .
  * Dobaczewski  _et al._ [1984] J. Dobaczewski, H. Flocard, and J. Treiner, Hartree-fock-bogolyubov description of nuclei near the neutron-drip line, Nucl. Phys. A 422, 103 (1984).
  * Kasuya and Yoshida [2021] H. Kasuya and K. Yoshida, Hartree–Fock–Bogoliubov theory for odd-mass nuclei with a time-odd constraint and application to deformed halo nuclei, Prog. Theor. Exp. Phys. 2021, 013D01 (2021), arXiv:2005.03276 [nucl-th] .
  * van Giai and Sagawa [1981] N. van Giai and H. Sagawa, Spin-isospin and pairing properties of modified Skyrme interactions, Phys. Lett. B 106, 379 (1981).
  * Yamagami  _et al._ [2009] M. Yamagami, Y. R. Shimizu, and T. Nakatsukasa, Optimal pair density functional for description of nuclei with large neutron excess, Phys. Rev. C 80, 064301 (2009), arXiv:0812.3197 [nucl-th] .
  * Richter  _et al._ [1990] A. Richter, A. Weiss, O. Haüsser, and B. A. Brown, New evidence for meson-exchange-current enhancement of isovector M1 strength, Phys. Rev. Lett. 65, 2519 (1990).
  * Pringle  _et al._ [1982] D. Pringle, E. Garman, S. Chew, K. Snover, W. Catford, S. Hesmondhalgh, and K. Allen, Decay of the lowest $T=2$ state in ${}^{40}\mathrm{Ca}$, Phys. Lett. B 115, 291 (1982).
  * Kamerdzhiev and Tkachev [1989] S. P. Kamerdzhiev and V. N. Tkachev, A microscopic model taking into account 2p2h configurations in magic nuclei. Calculations of M1 excitations, Z. Phys. A 334, 19 (1989).
  * Bohr and Mottelson [1969] A. Bohr and B. Mottelson,  _Nuclear Structure: Volume II, Nuclear Deformations_ (Benjamin, 1969).
  * Hamamoto and Sagawa [1999] I. Hamamoto and H. Sagawa, Response of light drip line nuclei to spin dependent operators, Phys. Rev. C 60, 064314 (1999).
  * Kamerdzhev  _et al._ [1993] S. Kamerdzhev, J. Speth, G. Tertychny, and J. Wambach, M1 resonances in unstable magic nuclei, Z. Phys. A 346, 253 (1993), arXiv:nucl-th/9303001 .
  * Oishi  _et al._ [2020] T. Oishi, G. Kružić, and N. Paar, Role of residual interaction in the relativistic description of M1 excitation, J. Phys. G 47, 115106 (2020), arXiv:1912.11147 [nucl-th] .
  * Matsuo [2001] M. Matsuo, Continuum linear response in coordinate space Hartree-Fock-Bogolyubov formalism for collective excitations in dripline nuclei, Nucl. Phys. A 696, 371 (2001), arXiv:nucl-th/0105025 .
  * Matsuo  _et al._ [2005] M. Matsuo, K. Mizuyama, and Y. Serizawa, Di-neutron correlation and soft dipole excitation in medium mass neutron-rich nuclei near drip-line, Phys. Rev. C 71, 064326 (2005), arXiv:nucl-th/0408052 .
  * Yoshida  _et al._ [2006] K. Yoshida, M. Yamagami, and K. Matsuyanagi, Pairing and continuum effects on low-frequency quadrupole vibrations in deformed Mg isotopes close to the neutron drip line, Nucl. Phys. A 779, 99 (2006), arXiv:nucl-th/0605073 .
  * Yoshida and Yamagami [2008] K. Yoshida and M. Yamagami, Low-frequency $K^{\pi}=0^{+}$ modes in deformed neutron-rich nuclei: Pairing- and $\beta$-vibrational modes of neutron, Phys. Rev. C 77, 044312 (2008), arXiv:0802.2341 [nucl-th] .
  * Lipparini and Stringari [1989] E. Lipparini and S. Stringari, Sum rules and giant resonances in nuclei, Phys. Rep. 175, 103 (1989).
  * Alberico  _et al._ [1982] W. Alberico, A. Molinari, R. Cenni, and M. B. Johnson, On the collective modes of infinite nuclear matter, Ann. Phys. 138, 178 (1982).
  * Belyaev [1965] S. Belyaev, Time-dependent self-consistent field and collective nuclear hamiltonian, Nucl. Phys. 64, 17 (1965).
  * Sagawa  _et al._ [2016] H. Sagawa, C. L. Bai, and G. Colò, Isovector spin-singlet ($T=1,S=0$) and isoscalar spin-triplet ($T=0,S=1$) pairing interactions and spin-isospin response, Phys. Scr. 91, 083011 (2016), arXiv:1512.07589 [nucl-th] .
  * Sagawa  _et al._ [1993] H. Sagawa, I. Hamamoto, and M. Ishihara, Gamow-Teller states in light nuclei near the neutron drip line, Phys. Lett. B 303, 215 (1993).
  * Yoshida [2017] K. Yoshida, Charge-exchange dipole excitations in neutron-rich nuclei: $-1\hbar\omega_{0}$ , anti-analog pygmy and anti-analog giant resonances, Phys. Rev. C 96, 051302 (2017), arXiv:1709.10272 [nucl-th] .
  * Davesne  _et al._ [2019] D. Davesne, A. Pastore, and J. Navarro, Linear response theory in asymmetric nuclear matter for Skyrme functionals including spin-orbit and tensor terms II: Charge Exchange, Phys. Rev. C 100, 064301 (2019), arXiv:1905.12049 [nucl-th] .
  * Yoshida [2019] K. Yoshida, Suddenly shortened half-lives beyond 78Ni: $N=50$ magic number and high-energy non-unique first-forbidden transitions, Phys. Rev. C 100, 024316 (2019), arXiv:1903.03310 [nucl-th] .


# Galaxy Deblending using Residual Dense Neural networks

Hong Wang  Stony Brook University, Stony Brook, NY 11794  Computational
Science Initiative, Brookhaven National Laboratory, Upton, NY 11973
Sreevarsha Sreejith  Physics Department, Brookhaven National Laboratory,
Upton, NY 11973  Yuewei Lin  Computational Science Initiative, Brookhaven
National Laboratory, Upton, NY 11973  Anže Slosar  Physics Department,
Brookhaven National Laboratory, Upton, NY 11973  Shinjae Yoo  Computational
Science Initiative, Brookhaven National Laboratory, Upton, NY 11973

###### Abstract

We present a new neural network approach for deblending galaxy images in
astronomical data using Residual Dense Neural network (RDN) architecture. We
train the network on synthetic galaxy images similar to the typical
arrangements of field galaxies with a finite point spread function (PSF) and
realistic noise levels. The main novelty of our approach is the usage of two
distinct neural networks: i) a _deblending network_ which isolates a single
galaxy postage stamp from the composite and ii) a _classifier network_ which
counts the remaining number of galaxies. The deblending proceeds by
iteratively peeling one galaxy at a time from the composite until the image
contains no further objects as determined by classifier, or by other stopping
criteria. By looking at the consistency in the outputs of the two networks, we
can assess the quality of the deblending. We characterize the flux and shape
reconstructions in different quality bins and compare our deblender with the
industry standard, SExtractor. We also discuss possible future extensions for
the project with variable PSFs and noise levels.

††preprint: arXiv:xxxx.xxxx

##  I Introduction

Astronomical survey science is entering a decade which will witness a surfeit
of new optical imaging data that will transform astronomical and cosmological
research. The new datasets obtained from upcoming ground-based and space-based
observing facilities will image vast areas of the sky at unprecedented depths.
Among the most important of these will be the Legacy Survey of Space and Time
(LSST) from the Vera Rubin Observatory [1], which will complement not only the
already existing datasets from Dark Energy Survey (DES, [2]) and Hyper-Suprime
Cam (HSC,[3]), but also the data from upcoming space-based surveys like Euclid
[4]. Due to the large field of view, most of these surveys will not be able to
employ adaptive optics resulting in an arcsec sized point spread function
(PSF). At the depths of these surveys, this will result in blending, affecting
over half of the galaxies in the survey (for example, in [5], the authors
estimate that approximately 60% of the objects in the HSC are affected by this
problem).

Blending, especially if undetected, can introduce serious systematic errors in
the survey analysis. These include potentially catastrophically wrong
inference of fluxes, leading to biased photometric redshifts, biased estimates
of shear and local object density dependent systematics that can interact with
galaxy and cluster detections. The main approach to solving the blending
problem is two-pronged. On one hand, there is the realization that source
separation (deblending) will never be perfect and we should focus on
understanding its properties and its effects on the analysis. This is done
through simulations and artificial source injections into real data. On the
other hand, the better the deblender we have, the smaller the corrections we
need to do to make our analysis unbiased. We therefore need both good
deblenders and good schemes for understanding their imperfections.

Deblending is in principle a well defined problem. The basic model is that the
images of individual galaxies are combined on the projected plane, assuming
perfect transparency (i.e. intensities add), and then observed through the
atmosphere and telescope with known PSF and noise properties. The most
sophisticated deblenders on the market combine machine-learning approaches for
setting priors on galaxy shapes, and physical modelling for things that we can
model explicitly, like the noise and the PSF [6, 7].

Several machine-learning based methods have been proposed recently, in order
to grapple with the galaxy deblending problem. [8] designed a branched
deblender with generative adversarial networks, which can work on the
deblending of images with two overlapped galaxies. [9] developed a framework
for measuring the photometry of blended galaxies as well as perform
segmentation with a standard convolutional neural network (CNN) and a U-Net.
In 2020, [5] introduced an algorithm where a Variational Auto
Encoder(VAE)-like neural network was used for galaxy deblending. Most of the
current generation neural-network based galaxy deblenders require the galaxy
to be located at the center of the image in order to give the network a sign
as to which galaxy to recover, which is often impractical and prevents the
model from being used iteratively, i.e., it will degrade the model’s
performance on the image after the first galaxy is removed. Besides, the
number of galaxies in the blended image for the previous methods is always
fixed, and to the best of our knowledge, there is no neural network based
galaxy deblending framework that can work on images with an arbitrary and
unknown number of galaxies. In this paper, we propose an innovative deblending
framework with a galaxy deblender, and a classifier that can deblend galaxies
from a blended image iteratively without any prior information about the
number of galaxies. In addition, since the deblender recovers galaxies based
on their luminosity, our framework has no constraint on the position of the
galaxies. Nevertheless, like other machine-learning approaches to deblending,
our work remains exploratory as we aim to better understand the applicability
of neural networks to astronomical image analysis.

This paper is organized as follows: in Section II, we introduce the
architecture of our framework including the deblender and the classifier, and
the experimental settings used to train the model. In Section III, we present
the experimental results and the comparison with the industry standard
deblending method - Source Extractor (SExtractor, [10]). The discussion and
conclusion follow, in the last Section IV.

##  II Neural Network Architecture and Training

###  II.1 Proposed Framework

Our goal is to deblend galaxies from astronomical images with an arbitrary and
unknown number of overlapped galaxies. The proposed framework consists of two
components: a _deblender_ , which isolates the image corresponding to a single
galaxy from an astronomical image, and a _classifier_ , which counts how many
galaxies remain in the image. The _deblender_ and _classifier_ are then used
iteratively to separate the scene into its constituent galaxy images. This is
illustrated in the Figure 1 and represented in the following meta-code:

    while True:
        num_galaxies = classifier(image)
        if num_galaxies == 0:
            exit loop
        deblended_galaxy = deblender(image)
        image -= deblended_galaxy

Figure 1: Framework.

Specifically, given a noisy blended image with multiple galaxies, the
deblender will take it as the input and output a noiseless image with a single
galaxy as it is trained. This single-galaxy image is then subtracted from the
input of the deblender to get the residual image which is another noisy
blended image. The classifier is used to determine if there are further
galaxies in the residual image, and if there are, this process is repeated
until there is no galaxy left in the residual image.

In an ideal case, the classifier will detect one fewer galaxy at each step and
the process will stop when there are no more galaxies left. We call such
deblends _High quality_ deblends. In a non-ideal case, one of the two
following scenarios often play out. The number of galaxies in the image does
not decrease by one galaxy per iteration, but the process still eventually
comes to a halt with zero galaxies detected in the final image. When this
happens, we refer to the result as _Medium quality_ deblends. The third case
is if the process gets stuck in an infinite loop where the _classifier_
maintains that there are more than zero galaxies in the image, but the
deblender fails to locate them. Then, if the classifier predicts the same non-
zero number of galaxies in the residual image for three consecutive iterations
we terminate the iteration and call the results _Low quality_ deblends.

The next subsection describes the architectures of our two main components.

Figure 2: The architecture of RDN.

###  II.2 Deblender

Given an astronomical image with multiple overlapped galaxies, the deblender
aims to deblend one galaxy from it. A Residual Dense Network (RDN) [11] is
trained as the deblender in this framework. RDN shows superior performance in
image super-resolution [11] and image restoration [12]. For the deblending
task, the RDN will take noisy blended images with multiple galaxies as input
and the output will be a noiseless image containing the brightest galaxy.

Figure 2 shows the architecture of the RDN, starting with a shallow feature
extraction net (SFTNet), which includes two convolutional layers that extract
features from the input. The SFTNet is followed by the main component of the
RDN, namely the residual dense blocks (RDBs). Each RDB contains $C$
convolutional layers with a ReLU activation function. The layers in the RDB
are densely connected, which means the feature maps from all the preceding
layers of the current RDB as well as the output from the preceding RDB will be
the input of the current convolutional layer. The feature maps from the
convolutional layers in the current RDB and the preceding RDB are also
concatenated as a local feature fusion, followed by a $1\times 1$
convolutional layer. At the end of an RDB, there is a local residual step by
the addition of the output from the previous RDB. Following $D$ RDBs, where
$D$ denotes the number of RDBs, is a global feature fusion where, similar to
the local feature fusion, the features from all the RDBs will be concatenated
and then passed to a $1\times 1$ convolutional layer. The last step before the
up-sampling net is global residual learning, referring to an addition between
the current features and the features from the first convolutional layer in
SFTNet. The efficient sub-pixel convolutional neural network (ESPCN) followed
by a convolutional layer forms the up-sampling net. Since the RDN in our
framework is used for deblending, the output image will have the same size as
the input image. During training, the RDN will take noisy blended images as
input and predict noiseless images with a single galaxy. $\ell_{1}$-norm is
used to compute the difference between the output of the RDN and the ground
truth, thus the loss function for RDN is written:

$l^{\theta}_{\rm
RDN}(I,I_{gt})=\frac{1}{WH}\sum_{x=1}^{W}\sum_{y=1}^{H}\left|f(I,\theta)_{x,y}-(I_{gt})_{x,y}\right|$
(1)

where $I$ is the input noisy blended image and $f_{\theta}(I)$ is the output
of the RDN with $\theta$ referring to the parameters in the model. $I_{GT}$
denotes the ground truth image, which is a noiseless image with the brightest
galaxy.

In Figure 2 we show a schematic representation of the _deblender_ network.

###  II.3 Classifier

For the classifier portion of the deblending network, we use a VGG-16 network
[13]. VGG-16 contains 13 convolutional layers and 3 fully connected layers. In
our framework, it is modified to have four classes as output, i.e., the
classifier can tell if the images have 0, 1, 2 or 3 galaxies.

In the first phase of training, the classifier is trained on images with
galaxies ranging from 0 to 3. During the second phase in the training process
where the deblender and the classifier are trained jointly, the deblended
images from the deblender will also be passed through the classifier to see if
the RDN has deblended only one galaxy from the blended image. This loss will
also be used to update the RDN’s parameters. Thus, the training set of the
classifier contains both noiseless and noisy images with 0 to 3 galaxies. The
cross-entropy (CE) loss is utilized to train the classifier parameterized by
$\phi$ as Eq. 2 where $y$ represents the one-hot ground truth label.

$l^{\phi}_{CE}(I)=-y^{\rm T}\cdot\log(g(I,\phi))$ (2)

###  II.4 Training & Test Data

In this work, we use BlendingToolKit
111https://github.com/LSSTDESC/BlendingToolKit in order to simulate the galaxy
images and to generate the blended images. BlendingToolkit is a complete
framework whose functionalities include generation of images with blended
objects as well as providing measurements of the deblending and detection
performance. It relies on the industry standard galsim [15] packages to make
the actual renderings of the galaxies and is built on top of the
WeakLensingDeblending 222https://github.com/LSSTDESC/WeakLensingDeblending
package.

We first generate noiseless images of single galaxies, which can be considered
as the pre-blended ground truth images. We employ the LSST DM galaxy catalog
with a span of 1 square deg (OneDegSq.fits supplied with the BlendingToolKit).
We use the default options, which impose a magnitude cut on the i-band,
$i<25.3$. Then, the noiseless blended images are generated by a pixelwise
summation between these single galaxy images. Additional random Poisson noise
is added to the image to get the noisy blended images corresponding to the
LSST 10-year depth. We selected $g$, $r$ and $i$ bands from the resultant
images and converted them into RGB images. The original dimension of the
generated images is (120, 120, 3). We crop each image to (80, 80, 3) from the
center. The pixel values of the blended images will be normalized to $[0,1]$
before entering the framework and the pre-blended images will be scaled to
$[-1,1]$ following [8] .

Unlike the other machine learning based deblending methods where one of the
galaxies has to be located at the center of the image [8, 9], there is no
constraint on the positions of the galaxies for our deblender, because we
train it to output the object with the highest luminosity. This is also the
reason why our framework can work iteratively. The use of the summation when
creating the blended images also makes the iterative process possible. During
training, the noisy blended images are treated as the input and the noiseless
image with the brightest galaxy as the pre-blended ground truth for the
deblender. Our goal is to train a deblender that can deblend and denoise at
the same time and will always recover the brightest galaxy from the blended
images. The classifier is trained with $0-3$ galaxies and for each of the four
classes, half of the training data contain noiseless images while the other
half contain noisy images. This is essential in order to train a classifier
that can tell the number of galaxies in both noiseless and noisy images at the
same time.

We generate 50,000 2-galaxy blended images with BlendingToolKit as the
training set for the deblender. The test set contains 1000 blended images with
2 galaxies in each. There is also a second test set with 1000 3-galaxy blended
images, to demonstrate the generalization of our framework. For the
classifier, we generated blended images containing 0, 1, 2 and 3 galaxies,
where each of the classes has 100,000 images (making it 400,000 in total),
half of which are noiseless images and the other half noisy. The test set for
the classifier contains 2000 images for each class and 8000 in total. The test
set for the classifier also has the same ratio of noiseless and noisy images.
Therefore, the classifier gives 4 classes as output indicating the number of
galaxies in the input image.

###  II.5 Training procedure

The training process can be divided into two phases. In the first phase, the
deblender and the classifier are pre-trained separately until they perform
reasonably. Then they are fine-tuned jointly to boost the performance of both,
while working in the iterative setting. Specifically, we expect that the
classifier will force the deblender to output only one galaxy, while the
deblender will provide more intermediate images to motivate the classifier to
have better discrimination ability.

During the pre-training of the RDN, the noisy two-galaxy blend images are used
as the input and the noiseless images with the brighter single galaxy as the
ground truth. Eq. 1 is the objective function during this process. A model
trained in this way will have the ability to not only deblend the brighter
galaxy, but also to denoise simultaneously. Meanwhile, in the first training
phase the classifier is trained on the dataset containing both noisy and
noiseless blended images with 0, 1, 2 and 3 galaxies. This pre-training phase
makes the trained classifier able to classify both noisy and noiseless images
and the trained RDN to determine the end of the whole workflow. They will be
used during together during the second phase of the training.

During the second phase, the deblender and the classifier will be trained
jointly. The framework will operate as it is designed, i.e., the deblender
will be run twice because the input contains two galaxies. The blended image,
denoted as $I_{{\rm blend}}$ will be passed through the deblender to get the
deblended $I_{{\rm deblend}-1}$, which is compared with the ground truth
$I_{{\rm gt}-1}$ to build the $\ell_{1}$-norm loss $l_{\rm RDN}(I_{{\rm
deblend}-1},I_{{\rm gt}-1})$. The first residual image $I_{{\rm res}-1}$ is
the subtraction between $I_{{\rm blend}}$ and $I_{{\rm deblend}-1}$.
Generally, the residual images are fainter than the blended images and
therefore they will be normalized by dividing with the maximum pixel values.
Then, the deblender will predict and scale back to recover the second galaxy
$I_{{\rm deblend}-2}$ from $I_{{\rm res}-1}$ and another residual image
$I_{{\rm res}-2}$ will be calculated. For simplicity, we use
$\mathbb{S}=\\{I_{{\rm blend}},I_{{\rm deblend}-1},I_{{\rm deblend}-2},I_{{\rm
res}-1},I_{{\rm res}-2},I_{{\rm gt}-1},I_{{\rm gt}-2}\\}$ to denote all the
images. In this phase, the classifier will be optimized by all the available
images in $\mathbb{S}$ as formulated in Eq. 3. Since $\mathbb{S}$ is class-
imbalanced with more 1-galaxy images than 0 and 2-galaxy images, $\alpha_{I}$
is used as a weight. For 1-galaxy images, $\alpha_{I}=0.2$, otherwise
$\alpha_{I}=1$. For the deblender, in addition to $l_{\rm RDN}(I_{{\rm
deblend}-1},I_{{\rm gt}-1})$, it is updated based on a loss from the
classifier’s prediction on the deblended images. The objective function for
the second phase can be formulated as Eq. 4. $\lambda$ here is a trade-off
coefficient.

$l^{\phi}_{phase-2}=\sum_{I\in\mathbb{S}}\alpha_{I}l_{CE}(I)$ (3)

$l^{\theta}_{phase-2}=l_{\rm RDN}(I_{{\rm deblend}-1},I_{{\rm gt}-1})+\lambda
l_{CE}(I_{{\rm deblend}})$ (4)

Figure 3: Iterative results for 2-galaxy blended images. We show two examples.
In each example, the top row shows the noisy image as it progresses through
the iterative process with the input image on the left and the remaining noise
image (after the galaxies were subtracted) on the right. The middle row show
the galaxy images that were isolated from the input image by the _deblender_.
The bottom row shows the ground truth, which can be visually compared with the
deblended images above. All the images were plotted with inverse fluxes to
improve the contrast (i.e. white corresponds to zero flux). The number at the
corner of each image represents the number of remaining galaxies as determined
by the classifier.

Figure 4: As Fig. 3, but for medium and low quality samples.

Figure 5: As Fig. 3, but for 3 galaxy blends.

##  III Results

###  III.1 Experimental settings

####  III.1.1 Implementation details

The RDN used in our framework contains $D=16$ RDBs with $C=8$ convolutional
layers in each block. During the first phase of the training process, the
learning rate for the RDN is $10^{-4}$, the batch size is 128, and it is
trained for 150 epochs. For the classifier, the initial learning rate is
$0.1$, the batch size is 200, and it is trained for 200 epochs. In the second
phase, the deblender and the classifier are updated jointly, where the
learning rate decays to $10^{-5}$ for the former and $10^{-6}$ for the latter.
After some experimentation, we chose the trade-off coefficient
$\lambda=10^{-4}$. The batch size is 8 due to GPU limit and the framework is
trained for 10 epochs in this phase.

####  III.1.2 Evaluation

We trained the framework on 2-galaxy blended images and tested the trained
model on both 2-galaxy and 3-galaxy blended images. Both test sets contain
1000 noisy blended images.

We start by showing some example results in Figures 3 and 5. These results
show the typical output from a high-quality deblend. The network is using both
the morphological as well as the color information to isolate images of
individual galaxies which is helpful in recovering images even in difficult
cases.

Table 1:  Iterative test on 2-galaxy and 3-galaxy blended images. | 2-galaxy | 3-galaxy
---|---|---
High-quality | 77.3% | 50.4%
Medium-quality | 15.3% | 36.9%
Low-quality | 7.4% | 12.7%

In Table 1 we present the fraction of deblends sorted by the quality category
(see Section II.1) for 2-galaxy and 3-galaxy problems. Note that this only
indicates how well the scheme thinks it is doing, rather than how well it is
actually doing. With this caveat, the process shows relatively good results
with a large fraction of nominally high-quality deblends.

| 2-galaxy problem | 3-galaxy problem
---|---|---
| PSRN | SSIM | PSRN | SSIM
| Mean | Median | Mean | Median | Mean | Median | Mean | Median
Deblended 1 | 56.51 | 57.93 | 0.9973 | 0.9997 | 54.66 | 57.19 | 0.9954 | 0.9997
Deblended 2 | 58.18 | 59.02 | 0.9967 | 0.9997 | 55.96 | 58.77 | 0.9929 | 0.9996
Deblended 3 |  |  |  |  | 57.71 | 59.29 | 0.9956 | 0.9995
Table 2:  PSNR(dB) and SSIM for 2-galaxy and 3-galaxy deblending problems.

To quantify the quality of the deblends, we start with some standard image
analysis metrics. We applied the peak signal-to-noise ratio (PSNR) and the
structural similarity index (SSIM) as metrics to evaluate the quality of
deblended images from the deblender when compared with the ground truth.

PSNR represents a peak error in the measurement of reconstruction quality in
image compression. It computes the logarithm of the ratio between the maximum
pixel value (MAX) of the ground truth and the mean squared error (MSE) between
the test image and the ground truth in decibels. In our experiment, the ground
truth consists of noiseless images with a single galaxy and the test images
are the corresponding deblended images from the RDN. PSNR is formulated as Eq.
5:

${\rm PSNR(dB)}=20\cdot\log_{10}({\rm MAX})-10\cdot\log_{10}({\rm MSE})$ (5)

SSIM [17], known as the structural similarity index, is commonly used to
evaluate the similarity between two images using the means $\mu_{x}$ and
$\mu_{y}$, the variances $\sigma_{x}$ and $\sigma_{y}$ and the covariance
$\sigma_{xy}$. In Eq. 6, $c_{1}=(k_{1}L)^{2}$ and $c_{2}=(k_{2}L)^{2}$ are two
small constants to avoid the instability with a weak denominator. We use
$k_{1}=0.01$, $k_{2}=0.03$ by default and $L$ is the dynamic range of pixel
values.

${\rm
SSIM}=\frac{(2\mu_{x}\mu_{y}+c_{1})(2\sigma_{xy}+c_{2})}{(\mu_{x}^{2}+\mu_{y}^{2}+c_{1})(\sigma_{x}^{2}+\sigma_{y}^{2}+c_{2})}$
(6)

Results are summarized in Table 2. We do observe a few trends here. First, the
PSNR of the second deblended galaxy is higher than that of the first. This
somewhat counter-intuitive result comes from the normalization of PSNR, since
the second deblend is fainter than the first one. In other words, the PSNR is
not telling us that the 2nd deblended galaxy is better deblended, only that
its quality is less degraded than expected, given how much fainter it is. This
interpretation is confirmed by the SSIM values, which are higher for the first
deblend, but only marginally. In both cases we see that the 3-galaxy problem
is more difficult than the 2-galaxy problem. We also see that the median
values are systematically above the mean values, telling us that the mean is
pulled down by a small number of catastrophic outliers. To study these results
in language that is more relevant to the astronomy community, we turn to the
recovery of the fluxes and moments, described in the next section.

###  III.2 Flux and image moment recovery using RDN

Figure 6: Plot showing the comparison of the fluxes recovered from the RDN and
the true fluxes for the x galaxies. The red line is the linear regression fit
to the data and the dashed black line is an arbitrary $y=x$ fit. The shading
in green denotes a $95\%$ confidence interval to the linear fit.

An important measure of the effectiveness of any algorithm such as ours is the
accurate recovery of physical parameters from the input image. To that end, we
have opted to compare the fluxes and the second order image moments (which
represent the equivalent ellipse of the image and by extension, the shape and
orientation) retrieved by the RDN against that of the truth images. For ease
of reference, the brighter galaxies in each field are referred to as ‘x’
galaxies and the fainter ones as ‘y’ galaxies.

Figure 6 shows a comparison of the fluxes of the x galaxies as recovered by
the RDN against the truth images for the test set. The larger panel on the
left shows the true fluxes on the X-axis and the RDN fluxes on the y-axis. The
solid red line is a linear regression fit to the data, and the dashed black
line is an arbitrary $y=x$ fit. The light green shading represents a $95\%$
confidence interval to the linear fit. The three smaller panels on the right
side show the same fits to the high, medium and low quality blends which are
detailed in Section II.1. The overlapping histograms in the bottom left panel
show the density distributions of the RDN and true fluxes. Figure 7 is a
similar representation for the fainter, y galaxies.

Figure 7: As Fig. 6, but for the y galaxies.

As is evident from Figure 6, the RDN is very efficient in recovering the
fluxes for the brighter galaxies in our fields, more or less uniformly so for
the high, medium and low quality images. The effectiveness of the RDN for this
group of images can be further observed from the overlapping histograms of the
true (red) and recovered (green) fluxes in the bottom left panel. Figure 7
shows the same plots for the fainter, y galaxies. As expected, the lower
quality assignments are associated with a higher scatter between real and
measured fluxes, although for the brighter galaxy, there is no measurable
difference between the medium and low quality marks (given the size of our
test catalog). For the fainter galaxies, however, the low quality deblends are
essentially uncorrelated with the true values. We find that except for
secondary objects in the low-quality bin, our RDN based estimator is
essentially unbiased in recovering flux values.

Similarly to the fluxes, the shapes and orientations of the galaxies,
represented by their second order image moments, are crucial characteristics
that we aimed to recover using the RDN, in particular with the weak lensing
application in mind. Figure 8 shows a comparison of the moment $\mu_{11}$ for
the truth images and the RDN recovered images. The bigger panels show the true
moment on the x-axis and the RDN recovered moment on the y-axis for the x and
y galaxies. The blue line is an arbitrary $y=x$ fit. The smaller panels on the
bottom show the same plot for the high, medium and low quality images. Results
for the other component of the second order moment, $\mu_{20}-\mu_{02}$ looks
qualitatively the same. In an ideal deblender, the results would lie on the
$y=x$ line, but we notice the appearance of a curious cross, which we will
comment on below.

Figure 8: Plot showing the comparison of true $\mu_{11}$ and RDN recovered
$\mu_{11}$ for the brighter (x) and the fainter (y) galaxies on the large top
plots. The smaller second row plots show the same split by quality bin. The
blue line is a $y=x$ fit.  Figure 9: Plot showing the second order image
moments for x and y ground truth images in green, with the RDN recovered
moments superimposed in blue. The true and recovered moments exhibit the same
general distribution while some scatter is observed at the edges of the plot.
Figure 10: Plot showing the histogram of the absolute value of the second
moment $|\mu_{2}|$ for the true images. The true simulated galaxies belong to
two categories: they are either round $|\mu_{2}|=0$ or not. The histogram is
truncated with a small number of galaxies with considerably higher
$|\mu_{s}|$s.

Figure 9 shows a comparison of the second order moments ($\mu_{11}$,
$\mu_{20}-\mu_{02}$) plotted against each other for both the ground truth x
and y galaxies (green) with the moments recovered by the RDN superimposed
(blue). We see no evidence of any preferred axis in the recovered distribution
of these quantities.

This leads to a curious result. The ensemble properties of the RDN recovered
fluxes are correctly distributed. However, the per-galaxy results have
structure that goes beyond the normal scatter around the truth result. In
particular, the cross seen in Figure 8 implies that galaxies with a
significant $\mu_{11}$ sometimes end up having zero recovered $\mu_{11}$, but
also vice-versa, that galaxies with zero true $\mu_{11}$ occasionally end up
having a significant recovered $\mu_{11}$.

One plausible explanation for this curious result is as follows. In Figure 10
we plot the histogram of the truth value of
$|\mu_{2}|=\sqrt{(\mu_{20}-\mu_{02})^{2}+(2\mu_{11})^{2}}$, i.e. the
rotationally invariant measure of the non-circularity of the image. This
histogram shows a distinct peak at around $\mu=0$ region. The simulated
galaxies fall into two classes, those that are approximately round and those
which are not. The RDN network seems to have learned this and tries to
classify the galaxy into one of the two categories before reconstructing its
image. The arms of the crosses correspond to these failures: one for the non-
circular galaxies reconstructed as circular ones and the other for the
circular galaxies reconstructed as the non-circular ones.

###  III.3 Comparison with SExtractor

Figure 11: Plot showing the comparison of flux recovery from the RDN and
SExtractor.The top left and right panels show the true fluxes on the X-axis
and the RDN and SExtractor fluxes on the y-axis respectively for x galaxies.
The red line is a linear regression fit with a $95\%$ confidence interval
represented by light green shading. The dashed black line is the $y=x$ line.
The bottom left and right panels are the same representations for y galaxies.

In this section, we compare our deblending strategy with what is widely
considered as the industry standard, Source Extractor (SExtractor, [10]),
which has been the baseline detection, deblending and image extraction
software in astronomy for over two decades. It returns a set of user specified
parameters from the more extensive default parameter file, by following the
specific configurations defined by the user. In this work, we have used the
Python formulation of SExtractor, sep [18] with the settings employed for
DES333https://github.com/esheldon/sxdes.

Figure 12: Plot showing the second order image moments for x and y ground
truth images in green, with the SExtractor recovered moments superimposed in
blue. While the distribution is similar to that in Figure 9, the SExtractor
recovered moments exhibit much larger scatter than the RDN recovered moments.

We note that this is a fundamentally unfair comparison for several reasons.
First, the RDN has been trained for this particular set of galaxies and thus
internally employs a correctly tuned prior for the distribution of
morphologies, fluxes and ellipticities, while SExtractor employs a general
algorithm and is thus intrinsically more robust. Moreover, following DES
configuration, SExtractor is run on the co-added image (R+G+B) and is thus
missing the color information. Nevertheless, it is an appropriate sanity check
to measure what kind of improvements can be brought about by employing more
sophisticated methods.

Figure 11 shows the comparison of the flux recovery by the RDN and SExtractor
for both the x and y galaxies. The top left and right panels feature the RDN
recovered and SExtractor recovered fluxes on the y-axis and the true fluxes on
the X-axis for the x galaxies. The bottom left and right panels are the same
representations for the y galaxies. The red line is a linear regression fit to
the data, with the light green shading representing a $95\%$ confidence
interval. The dashed black line shows $y=x$. As is evident from the plots, the
RDN does a better job of recovering the object fluxes for both the x and y
galaxies. SExtractor has a tendency to put a disproportionate amount of flux
into the first detected objects, and hence the x images are biased high, while
the y images are biased low.

Measuring moments for SExtractor images in considerably more difficult due to
the presence of noise in them and the results are overly noisy. Correctly
subtracting the noise contribution in a fair manner goes beyond the scope of
this paper.

###  III.4 Failure modes

In this section we investigate the most common failure modes in the two
approaches by cherry-picking cases, in which one method performs correctly
while the other fails.

We manually looked at all the cases where RDN detects just a single galaxy. In
a majority of those, the secondary object was very faint and overlapped almost
perfectly with the primary, leading to the deblending being suspended because
the _classifier_ determined that there are no further objects. However, for
$12$ out of the $1000$ objects in the test set, for the deblended x galaxies,
the result clearly shows two objects, examples of which are shown in Figure
13. These problems are trivial for the SExtractor approach. A visual
inspection of all these images was conducted and it was observed that, both
the objects in $7$ of these fields are very similar in shape and, in $5$
fields, they had similar brightnesses. It stands to reason that the RDN has
some difficulty in deblending objects with similar structural parameters, most
likely because it cannot decide which is the “brighter” object, especially in
the presence of noise. In fact, when presented with such examples in training,
it simply pays a large loss penalty whenever it starts with the “wrong”
object, without a clear rule about how to pick one of the two objects if the
noise hides the identity of the bright ones. We will discuss this further in
Section IV.

Figure 13: Instances where the RDN fails, while SExtractor succeeds in
detecting and deblending the objects in the field. On the left is the original
blend, in the centre is the output of the RDN and on the right is the blend
with the SExtractor detections superimposed in red.

When SExtractor is run on the test set, it returns photometric parameters for
$979$ x galaxies and $729$ y galaxies (the x counterparts of all y galaxies
are present in the x set). As seen in Figure 13, this includes $9$ instances
where the RDN fails to deblend the different galaxies in the field. However,
overall, there are more instances where the RDN successfully deblends the x
and y galaxies, while SExtractor either doesn’t detect either galaxy, or only
detects the x galaxy, but not the the y galaxy, a few examples of which are
illustrated by Figure 14. We note that these are traditional blend merges,
where the deblender merges two distinct objects into a single one. Since human
eye is pretty good at detecting these, we could possibly improve upon this by
fine-tuning the SExtractor settings, potentially at the cost of artificially
shredding objects with substructure. These are also examples where color
information is most helpful.

Figure 14: Instances where SExtractor fails to detect and deblend, while the
RDN succeeds. On the far left is the original blend, the RDN result is
represented in the two central panels for the x and y galaxies, and the blend
with the SExtractor detection superimposed in red is on the far right panel.

##  IV Summary and Conclusions

In this work we present a new approach to deblending using a Residual Dense
Network, which was trained on blended galaxy images with a realistic PSF and
noise levels, and which has performed decently in deblending and recovering
object fluxes and shapes. Compared to previous works, our set up does not
assume that the object to be deblended is located at the center of the image.
We also do not need to assume in advance, the number of objects to be
deblended. By using two neural networks, one to perform deblending and one to
determine the number of remaining objects in the field, we can classify the
quality of deblends. We have shown that deblends that are designed to be
higher quality have indeed, less noisy fluxes and shapes.

In the most current deblending networks found in literature, the network is
trained to deblend the central galaxy. Given that centering in the presence of
noise presents its own problems, we have trained our network to deblend the
brightest galaxy remaining in the image. This works very well, where the
choice of the brightest is clear, but leads to a particular failure mode where
the network cannot decide where to start, when two objects are approximately
equally bright (see Figure 13). This might be a generic feature of any kind of
deblender that proceeds by peeling off one galaxy at a time. For example, even
if we choose the galaxy closer to the center, there would likely be an
equivalent failure mode when two galaxies are equally close to the center. The
same thing could happen for other initial assumptions that attempt to
distinguish between the galaxies. There are several approaches to solve this.
One possibility would be to use a symmetrized loss function that would train
the network to deblend either galaxy. Another possibility would be to have a
more sophisticated network that deblends multiple galaxies at once, perhaps
again with an appropriately symmetrized loss function.

We have found that our method outperforms the industry standard SExtractor on
flux recovery, subject to strong caveats mentioned in the Section III.3.
However, we found that the RDN recovered shapes are either very good or
catastrophically wrong. The catastrophical failures result in the cross-effect
discussed in Section III.2. One possible explanation for this involves
bimodality in the shapes of the truth objects and the RDN internally assigning
the image to the wrong class. It also serves as a warning for use of a deeply
non-linear deblender based on NN for work requiring precision shapes, such as
weak lensing. At the same time, such work could lead to potentially
interesting new work on galaxy-star separation and morphological
classification of galaxies.

Similarly to the majority of other neural network approaches, our deblender
cannot deal optimally with a PSF that is variable and is different for the
different input channels (for e.g. for different bands), nor can it deal with
a pixel mask. This is a significant deficiency that is usually swept under the
rug under the assumption that the network can be retrained for a set of
different PSF. However, in practice, we have a different PSF for different
bands and with the typical astronomical observations in five or six bands, the
number of possible PSF combinations becomes unmanageable. The PSF shape in
each band thus needs to be part of the input and network training. We are
leaving these important aspects for future work.

One of the main advantages of most ML approaches to deblending is that they
also automatically denoise and deconvolve the image and can in principle fill
in the missing information in case of pixel masks. For a sufficiently
sophisticated neural network, this would trivialize many typical operations
such as flux and shear estimation on the resulting truth images. However, it
would also require us to have a more advanced way of propagation of both
statistical and systematic uncertainties.

This work has demonstrated the feasibility of using RDN in astronomical image
analysis and has highlighted some of the unresolved issues with the current
machine learning approaches to the deblending problem in the case of realistic
astronomical images.

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research through
the SciDAC grant “Accelerating HEP Science: Inference and Machine Learning at
Extreme Scales”. We thank our collaborators for useful input.

## References

  * LSST Science Collaboration  _et al._ [2009] LSST Science Collaboration, P. A. Abell, J. Allison, S. F. Anderson, J. R. Andrew, J. R. P. Angel, L. Armus, D. Arnett, S. J. Asztalos, T. S. Axelrod, and et al., LSST Science Book, Version 2.0, arXiv e-prints , arXiv:0912.0201 (2009), arXiv:0912.0201 [astro-ph.IM] .
  * Abbott  _et al._ [2018] T. M. C. Abbott, F. B. Abdalla, A. Alarcon, J. Aleksić, S. Allam, S. Allen, A. Amara, J. Annis, J. Asorey, S. Avila, and et al., Dark Energy Survey year 1 results: Cosmological constraints from galaxy clustering and weak lensing, Phys. Rev. D 98, 043526 (2018), arXiv:1708.01530 [astro-ph.CO] .
  * Hikage  _et al._ [2019] C. Hikage, M. Oguri, T. Hamana, S. More, R. Mandelbaum, M. Takada, F. Köhlinger, H. Miyatake, A. J. Nishizawa, H. Aihara, R. Armstrong, J. Bosch, J. Coupon, A. Ducout, P. Ho, B.-C. Hsieh, Y. Komiyama, F. Lanusse, A. Leauthaud, R. H. Lupton, E. Medezinski, S. Mineo, S. Miyama, S. Miyazaki, R. Murata, H. Murayama, M. Shirasaki, C. Sifón, M. Simet, J. Speagle, D. N. Spergel, M. A. Strauss, N. Sugiyama, M. Tanaka, Y. Utsumi, S.-Y. Wang, and Y. Yamada, Cosmology from cosmic shear power spectra with Subaru Hyper Suprime-Cam first-year data, PASJ 71, 43 (2019), arXiv:1809.09148 [astro-ph.CO] .
  * Refregier  _et al._ [2010] A. Refregier, A. Amara, T. D. Kitching, A. Rassat, R. Scaramella, and J. Weller, Euclid Imaging Consortium Science Book, arXiv e-prints , arXiv:1001.0061 (2010), arXiv:1001.0061 [astro-ph.IM] .
  * Arcelin  _et al._ [2021] B. Arcelin, C. Doux, E. Aubourg, C. Roucelle, and LSST Dark Energy Science Collaboration, Deblending galaxies with variational autoencoders: A joint multiband, multi-instrument approach, MNRAS 500, 531 (2021), arXiv:2005.12039 [astro-ph.IM] .
  * Melchior  _et al._ [2018] P. Melchior, F. Moolekamp, M. Jerdee, R. Armstrong, A. L. Sun, J. Bosch, and R. Lupton, SCARLET: Source separation in multi-band images by Constrained Matrix Factorization, Astronomy and Computing 24, 129 (2018), arXiv:1802.10157 [astro-ph.IM] .
  * Lanusse  _et al._ [2019] F. Lanusse, P. Melchior, and F. Moolekamp, Hybrid Physical-Deep Learning Model for Astronomical Inverse Problems, arXiv e-prints , arXiv:1912.03980 (2019), arXiv:1912.03980 [astro-ph.IM] .
  * Reiman and Göhre [2019] D. M. Reiman and B. E. Göhre, Deblending galaxy superpositions with branched generative adversarial networks, Monthly Notices of the Royal Astronomical Society 485 (2019).
  * Boucaud  _et al._ [2020] A. Boucaud, C. Heneka, E. E. Ishida, N. Sedaghat, R. S. de Souza, B. Moews, H. Dole, M. Castellano, E. Merlin, V. Roscani, _et al._ , Photometry of high-redshift blended galaxies using deep learning, Monthly Notices of the Royal Astronomical Society 491, 2481 (2020).
  * Bertin and Arnouts [1996] E. Bertin and S. Arnouts, SExtractor: Software for source extraction., A&AS 117, 393 (1996).
  * Zhang  _et al._ [2018] Y. Zhang, Y. Tian, Y. Kong, B. Zhong, and Y. Fu, Residual dense network for image super-resolution, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2018) pp. 2472–2481.
  * Zhang  _et al._ [2020] Y. Zhang, Y. Tian, Y. Kong, B. Zhong, and Y. Fu, Residual dense network for image restoration, IEEE Transactions on Pattern Analysis and Machine Intelligence (2020).
  * Simonyan and Zisserman [2014] K. Simonyan and A. Zisserman, Very deep convolutional networks for large-scale image recognition, arXiv preprint arXiv:1409.1556 (2014).
  * Note [1] https://github.com/LSSTDESC/BlendingToolKit.
  * Rowe  _et al._ [2015] B. T. P. Rowe, M. Jarvis, R. Mandelbaum, G. M. Bernstein, J. Bosch, M. Simet, J. E. Meyers, T. Kacprzak, R. Nakajima, J. Zuntz, H. Miyatake, J. P. Dietrich, R. Armstrong, P. Melchior, and M. S. S. Gill, GALSIM: The modular galaxy image simulation toolkit, Astronomy and Computing 10, 121 (2015), arXiv:1407.7676 [astro-ph.IM] .
  * Note [2] https://github.com/LSSTDESC/WeakLensingDeblending.
  * Wang  _et al._ [2004] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE transactions on image processing 13, 600 (2004).
  * Barbary  _et al._ [2015] K. Barbary, K. Boone, and C. Deil, Sep: V0.3.0 (2015).
  * Note [3] https://github.com/esheldon/sxdes.