# A New Galaxy Cluster Merger Capable of Probing Dark Matter: Abell 56

David Wittman Department of Physics and Astronomy, University of California,
Davis, CA 95616 USA  Rodrigo Stancioli Department of Physics and Astronomy,
University of California, Davis, CA 95616 USA  Kyle Finner IPAC, California
Institute of Technology, 1200 E California Blvd., Pasadena, CA 91125, USA
Faik Bouhrik California Northstate University, 2910 Prospect Park Dr, Rancho
Cordova, CA 95670 USA  Reinout van Weeren Leiden Observatory, Leiden
University, Niels Bohrweg 2, 2333 CA, Leiden, the Netherlands  Andrea Botteon
INAF - IRA, via P. Gobetti 101, I-40129 Bologna, Italy

###### Abstract

We report the discovery of a binary galaxy cluster merger via a search of the
redMaPPer optical cluster catalog, with a projected separation of 535 kpc
between the BCGs. Archival XMM-Newton spectro-imaging reveals a gas peak
between the BCGs, suggesting a recent pericenter passage. We conduct a galaxy
redshift survey to quantify the line-of-sight velocity difference ($153\pm
281$ km/s) between the two subclusters. We present weak lensing mass maps from
archival HST/ACS imaging, revealing masses of $M_{200}=4.5\pm 0.8\times
10^{14}$ and $2.8\pm 0.7\times 10^{14}$ M⊙ associated with the southern and
northern galaxy subclusters respectively. We also present deep GMRT 650 MHz
data revealing extended emission, 420 kpc long, which may be an AGN tail but
is potentially also a candidate radio relic. We draw from cosmological n-body
simulations to find analog systems, which imply that this system is observed
fairly soon (60-271 Myr) after pericenter, and that the subcluster separation
vector is within 22∘ of the plane of the sky, making it suitable for an
estimate of the dark matter scattering cross section. We find
$\frac{\sigma_{\rm DM}}{m}=1.1\pm 0.6$ cm2/g, suggesting that further study of
this system could support interestingly tight constraints.

Galaxy clusters (584); Dark matter (353); Galaxy spectroscopy (2171); Weak
gravitational lensing (1797); Hubble Space Telescope (761)

††facilities: Keck:II (Deimos), HST (ACS), GMRT, XMM††software: SAS (v19.0.0;
Gabriel et al. 2004), mc3gmm code (Golovich et al. 2019), FIATMAP code
(Wittman et al. 2006), SExtractor (Bertin & Arnouts 1996)

##  1 Introduction

A collision of two galaxy clusters dramatically reveals the contrasting
behaviors of gas, galaxies, and dark matter (DM). Seminal papers on the Bullet
Cluster provided a “direct empirical proof of dark matter” (Clowe et al.,
2006) as well as limits on the scattering cross-section of DM particles with
each other (Markevitch et al., 2004; Randall et al., 2008), aka DM “self-
interaction.” The Bullet constraint, $\frac{\sigma_{\rm DM}}{m}<0.7$ cm2/g, is
still quite large in particle physics terms—roughly at the level of neutron-
neutron scattering. In principle, ensembles of merging clusters enable tighter
constraints (Harvey et al., 2015; Wittman et al., 2018b) but these are
complicated by the fact that few systems have well-modeled dynamics.
Specifically, the time since pericenter, pericenter speed, and viewing angle
cannot be extracted from systems that have more than two merging subclusters.
Even with binary mergers, other factors may hinder the study of dark matter,
such as a merger axis closer to the line of sight rather than the plane of the
sky or bright stars that limit deep optical observations. Hence there is
interest in finding more “clean” binary systems with merger axis close to the
plane of the sky.

Figure 1: Abell 56: 0.4-1.25 keV XMM-Newton contours over SDSS multiband
(left) and HST/ACS F814W (right) images.

Historically, merging systems were discovered upon notice of disturbed X-ray
morphology, which typically happened serendipitously in pointed observations.
Meanwhile, modern optical sky surveys find tens of thousands of clusters and
promise to find many more as they get wider and deeper (Racca et al., 2016;
Ivezić et al., 2019). These surveys potentially contain new binary mergers, if
appropriate cuts can filter out tens of thousands of more ordinary clusters.
We have developed a new selection method based on the redMaPPer (Rykoff et
al., 2014) cluster catalog, which is in turn based on the ${\sim}$10,000 deg2
Sloan Digital Sky Survey (SDSS; York et al., 2000) imaging. We select clusters
not dominated by a single brightest cluster galaxy (BCG), in which there is
substantial angular separation between the top BCG candidates. These clusters
become candidate mergers, which are then checked against archival XMM-Newton
and Chandra data where available; if the X-ray peak is between the BCGs, the
candidate is worthy of additional followup. Because these two X-ray archives
consist of pointed observations rather than a uniform survey, our initial
candidates do not form a sample with well-defined selection criteria.
Nevertheless, a few initial candidates are worthy of immediate study in their
own right. In this paper we present the first such candidate, RM
J003353.1-075210.4, which we identify as Abell 56 as explained in §2.
Additional sections in this paper present a galaxy redshift survey of the
system (§3); a weak lensing analysis (§(1)); a search for analog systems in a
cosmological simulation (§5); radio observations in search of a radio relic
that could outline a merger shock (§6); and a constraint on the dark matter
scattering cross section $\frac{\sigma_{\rm DM}}{m}$(§7). We assume a flat
$\Lambda$CDM cosmology with $H_{0}=69.6$ km/s and $\Omega_{m}=0.286$.

##  2 Abell 56: initial overview

Nomenclature. The redMaPPer designation for this cluster is RM
J003353.1-075210.4. The original coordinates for Abell 56 (Abell et al., 1989
; hereafter ACO) are nearly $5^{\prime}$ north of the redMaPPer position.
Although ACO cite the positional uncertainty as $2.5^{\prime}$, inspection of
the SDSS imaging111https://skyserver.sdss.org/dr16/en/tools/chart/navi.aspx
reveals no other clusters in the area, suggesting that the ACO coordinates are
off by more than their nominal uncertainty. (The limited depth of the ACO
catalog is such that any real ACO cluster must be in the redMaPPer catalog.)
Indeed, the widely-used SIMBAD database222http://simbad.u-strasbg.fr/simbad/
resolves the name “Abell 56” to the redMaPPer position, as does the SDSS
Navigator noted above. Hence we identify this cluster as Abell 56. Note,
however, that the NASA/IPAC Extragalactic Database (NASA/IPAC Extragalactic
Database (2019), NED) resolves this name to the original ACO coordinates.

This cluster has also been detected by the Planck Sunyaev-Zel’dovich survey
(Planck Collaboration et al., 2016), with the designation PSZ2 G109.99-70.28.
This gas peak position is $1.3^{\prime}\pm 2.4^{\prime}$ from the redMaPPer
position and $4.0^{\prime}\pm 2.4^{\prime}$ from the ACO position. Hence
Planck Collaboration et al. (2016) adopted the redMaPPer position while
adopting “ACO 56” as the identifier in their union catalog. As a result, a
search for G109.99-70.28 on SIMBAD yields the redMaPPer position; however NED
yields the much more uncertain gas peak position.

BCGs and redshifts. Figure 1 presents two views of Abell 56: with XMM-Newton
contours (see below) over SDSS multiband imaging and over a single-band
(F814W) archival image from the Hubble Space Telescope Advanced Camera for
Surveys (HST/ACS) (see §(1)). At the redMaPPer photometric redshift of 0.30,
the physical scale is 4.5 kpc/arcsec. There are two galaxy subclusters
separated by close to $2^{\prime}$ (530 kpc), with the X-ray peak located
along the subcluster separation vector, about $32^{\prime\prime}$ (140 kpc)
from the southern subcluster. These numbers will be refined with further data
in later sections of this work.

The southern subcluster is dominated by a galaxy observed by the Baryon
Oscillation Spectroscopic Survey (BOSS; Dawson et al., 2013) to be at
$z=0.30231$; redMaPPer assigns this galaxy 84% probability of being the
overall BCG. The northern subcluster has two galaxies that appear nearly
equally salient in Figure 1; the eastern one (about 0.2 mag brighter) is
assigned 16% probability of being the overall BCG and BOSS places it at
$z=0.30475$. redMaPPer technically assigns some nonzero BCG probability to the
second-brightest galaxy in each subcluster, but in each case it is only
$0.017\%$, which we consider negligible.

Repp & Ebeling (2018) briefly considered this cluster as part of an 86-cluster
sample. It is classified in their Table 6 as being in the most disturbed of
their four optical morphology classes.

Merger basics. Taking the BCGs as tracers for a first calculation of the
merger geometry, we find a projected separation of 118′′ (535 kpc) and a line-
of-sight velocity difference of 565 km/s. For comparison, the projected
separation in the X-ray selected Bullet cluster is 720 kpc (Bradač et al.,
2006; Clowe et al., 2006) and those in Golovich et al. (2019) radio-selected
sample of merging clusters are generally closer to 1 Mpc, indicating that more
time since pericenter (TSP) has passed. This suggests the potential of optical
selection to find systems with smaller separations hence less TSP. Because a
complete understanding of the merging process will require snapshots of
systems spanning a range of TSP, optical selection may find its place within a
range of complementary selection methods.

Richness and related estimates. Rykoff et al. (2014) give the optical richness
$\lambda$ (a measure of how many galaxies are in the cluster, within a certain
luminosity range below the BCG) as 128\. Simet et al. (2017) calibrated the
relation between weak lensing mass and $\lambda$ (including its scatter), from
which we estimate the mass of Abell 56 to be $M_{200}=10.4^{+7.6}_{-4.6}\times
10^{14}\ h^{-1}$ M⊙. Sereno & Ettori (2017) implemented a system for mass
forecasting with proxies, taking into account various biases, and found
$M_{200}=11.49\pm 0.89\times 10^{14}$ M⊙ for this system based on its
redMapper richness. They also found $M_{500}=7.09\pm 0.77\times 10^{14}$ M⊙
using $Y_{\rm 500}$, a measure of the Sunyaev-Zel’dovich effect, as a proxy.
For comparison, Planck Collaboration et al. (2016) found
$M_{500}=5.62^{+0.54}_{-0.58}\times 10^{14}$ M⊙ from their scaling relation
based on the same $Y_{\rm 500}$ measurement.

From the scaling relations of Rozo & Rykoff (2014) one would expect the X-ray
temperature $T_{X}$ to be around 7 keV with up to 40% scatter at fixed
richness. Because this is a merging cluster, the X-ray properties may vary
from the scaling relations even more than usual.

X-ray properties from archival data. The cluster was observed with the XMM-
Newton European Photon Imaging Camera (EPIC) in 2010 (Obs.ID 0650380401, P.I.
Allen). The exposure times were 7121 s, 7127 s, and 5533 s for the MOS1, MOS2,
and PN instruments, respectively. As the short exposure does not allow for a
detailed analysis of the intracluster medium (ICM) properties and the cluster
morphology, we restricted our analysis to obtaining a point-source-subtracted,
exposure-corrected image, as well as a global temperature and luminosity for
the cluster.

We performed the data reduction using the XMM-Newton Science Analysis System
(SAS) version 19.0.0. We excluded periods of high soft-proton background by
imposing a cutoff of 0.4 (0.8) on the soft-proton rate for the MOS (PN)
detectors333We define MOS (PN) soft-proton events as those with energy
$E{>}10$ ($10{<}E{<}12$) keV. The higher-than-usual baseline soft-proton rate
for this observation may result in significant residual soft-proton
contamination even after the exclusion of flare events. This is at least
partially mitigated by the background-subtraction strategy., which resulted in
filtered exposure times of 6817 s, 6621 s, and 4178 s for MOS1, MOS2, and PN,
respectively. Only single-to-quadruple events from MOS and single-to-double
events from PN were used in our analysis. Point-source detection and masking
were performed by the cheese routine from the ESAS package. The contours in
Figure 1 are from the 0.4-1.25 keV band after point-source masking and
exposure correction, using the procedure described in the XMM ESAS Cookbook
(Snowden & Kuntz, 2014) and adaptively smoothed using the adapt routine from
ESAS.

In order to obtain a background-subtracted spectrum, we used the double-
subtraction method described in Arnaud et al. (2002). We defined the source
region as a circle with a 90′′ radius centered on the cluster, whereas the
background was extracted from a slightly larger circular region away from the
cluster. Blank-sky files (Carter & Read, 2007) were used to mitigate the
effects of the spatial variation of background components across the detector.
For a detailed description of the method, we refer to Arnaud et al. (2002).
Using XSPEC (Arnaud, 1996), we fit the spectrum to an apec model, multiplied
by a phabs model to account for galactic absorption. We obtained a total
unabsorbed luminosity in the $0.5-10.0$ keV range of $L_{X}=3.8\pm 0.2\times
10^{44}$ erg/s and a temperature $T_{X}=5.9^{+1.1}_{-0.8}$ keV, where the
uncertainties represent the $90\%$ confidence intervals.

##  3 Redshift survey and clustering kinematics

###  3.1 Redshift survey

Observational setup. We observed Abell 56 with the DEIMOS multi-object
spectrograph (Faber et al., 2003) at the W. M. Keck Observatory on July 1,
2022 (UT). The DEIMOS field of view is approximately 16${}^{\prime}\times
4^{\prime}$, making it well suited to merging clusters when the long axis is
placed along the subcluster separation vector. We prepared two slitmasks with
approximately sixty 1′′ wide slits in each. Galaxies were selected for
targeting based on (i) a preference for brighter targets; and (ii) a
preference for galaxies likely to be in the cluster based on Pan-STARRS
photometric redshifts (Beck et al., 2021). Because the photometric redshifts
are imprecise, this approach naturally helps probe for potential
foreground/background structures that could affect the modeling of Abell 56.
Specifically, each Pan-STARRS photometric redshift $z_{\rm PS}$ has a
corresponding uncertainty $\sigma_{\rm PS}$ such that the likelihood of the
galaxy being in a cluster at redshift $z_{\rm cl}$ is

$\mathcal{L}\propto\frac{1}{\sigma_{\rm PS}}\exp\frac{(z_{\rm PS}-z_{\rm
cl})^{2}}{2\sigma_{\rm PS}^{2}}$ (1)

The median value of $\sigma_{\rm PS}$ was 0.16, so a broad range of redshifts
was included. We then upweighted brighter galaxies by applying a
multiplicative weight $(24-r)$, where $r$ is the apparent $r$ magnitude, to
quantify the priority of each galaxy as input to the slitmask design software
dsimulator; larger numbers indicate higher priority. We manually raised the
priority of a few galaxies that potentially formed a foreground group at the
north end of the field.

We used the 1200 line mm-1 grating, which results in a pixel scale of 0.33 Å
pixel-1 and a resolution of ${\sim}1$ Å (50 km/s in the observed frame). The
grating was tilted to observe the wavelength range $\approx$ 4200–6900 Å(the
precise range depends on the slit position), which at $z\approx 0.3$ includes
spectral features from the [OII] 3727 Å doublet to the magnesium line at 5177
Å. The total exposure time was 45 (77) minutes on the first (second) mask,
divided in three (four) exposures. The seeing was roughly 1″, with minor
variations over time.

Data reduction and redshift extraction. We calibrated and reduced the data to
a series of 1-D spectra using PypeIt (Prochaska et al., 2020, 2020). We
double-checked the arc lamp wavelength calibration against sky emission lines,
and found good agreement.

To extract redshifts from the 1-D spectra we wrote custom Python software to
emulate major elements of the approach used by the DEEP2 (Newman et al., 2013)
survey using the same instrument. The throughput as a function of wavelength
varies from slit to slit, hindering direct comparison to template spectra.
Because throughput is generally a slowly varying function of wavelength, the
spectra are compared to templates only after removing the slowly varying
trends from each. First, telluric absorption features are reversed using Mauna
Kea models from the PypeIt development suite. Next, we create a smooth model
or unsharp mask by convolving the 1-D spectrum with a kernel 150 Å wide, which
is uniform but for a 10 Å diameter hole in the center. Finally, the intensity
of each pixel in the 1-D spectrum is expressed as a fraction of the intensity
in the smooth model.

The same operations are performed on redshifted versions of the galaxy
templates from the Sloan Digital Sky Survey444Available at
https://classic.sdss.org/dr7/algorithms/spectemplates/; we used templates 23
through 27., and a $\chi^{2}$ value is computed for each template-redshift
combination. A user then inspects the match between the data and the model
with the global minimum $\chi^{2}$, or other models at local minima, before
determining whether a redshift is secure. A secure redshift may not appear at
the global minimum due to poorly subtracted sky lines or other artifacts
(e.g., a spurious “line” with spuriously small uncertainties may appear at the
gap between CCDs). Furthermore, some slits suffer from vignetting at the red
end, which appears as a drop in intensity too steep to be removed by the
unsharp masking process. In these cases the user specifies a maximum
wavelength to consider for template matching. A negligible fraction of slits
contained stars, so we did not include stellar templates in the automated
search; users can manually classify a spectrum as a star without extracting a
redshift.

The uncertainty in the redshift is initially computed from the curvature of
the $\chi^{2}$ surface about the minimum, and is typically $\lesssim 10^{-4}$
(23 km/s in the frame of the cluster). We compared redshifts obtained by
different users on different computing hardware, operating systems, and Python
installations. We found that user-dependent uncertainty is also $\lesssim
10^{-4}$, mostly due to specification of the maximum wavelength. We therefore
add $10^{-4}$ in quadrature to the uncertainty derived from the curvature of
the $\chi^{2}$ surface to derive a final uncertainty estimate. We found 54 and
48 secure redshifts in the two masks respectively, for a total of 102. These
are listed in Table 3.1.

Table 1: Galaxy redshifts

RA (deg) | DEC (deg) | z | uncertainty
---|---|---|---
8.410133 | -7.758144 | 0.370305 | 0.000107
8.412567 | -7.747167 | 0.370005 | 0.000102
8.411033 | -7.742000 | 0.370055 | 0.000103
8.412454 | -7.750711 | 0.411192 | 0.000182
8.414025 | -7.767139 | 0.361506 | 0.000100
8.418721 | -7.738369 | 0.124575 | 0.000103
8.419221 | -7.752617 | 0.309164 | 0.000132
8.424929 | -7.768175 | 0.550809 | 0.000390
8.439262 | -7.728628 | 0.342070 | 0.000149
8.449196 | -7.770197 | 0.300759 | 0.000103
8.451654 | -7.772564 | 0.301399 | 0.000100
8.483425 | -7.776081 | 0.305369 | 0.000106
8.430800 | -7.878178 | 0.285649 | 0.000108
8.433154 | -7.865233 | 0.302410 | 0.000111
8.436871 | -7.800781 | 0.303911 | 0.000101
8.437129 | -7.820883 | 0.307063 | 0.000103
8.439008 | -7.874475 | 0.310816 | 0.000104
8.442646 | -7.833144 | 0.295355 | 0.000102
8.442646 | -7.833144 | 0.304912 | 0.000273
8.442892 | -7.838442 | 0.296056 | 0.000101
8.443604 | -7.795317 | 0.368704 | 0.000100
8.446533 | -7.918775 | 0.510949 | 0.000554
8.446804 | -7.847947 | 0.302537 | 0.000120
8.446900 | -7.865711 | 0.292854 | 0.000103
8.448988 | -7.939147 | 0.301960 | 0.000111
8.451133 | -7.910875 | 0.305312 | 0.000107
8.448746 | -7.942206 | 0.764905 | 0.000207
8.452063 | -7.952244 | 0.364001 | 0.000107
8.455208 | -7.976750 | 0.306263 | 0.000109
8.454325 | -7.842231 | 0.291603 | 0.000101
8.457279 | -7.888228 | 0.309365 | 0.000103
8.457279 | -7.888228 | 0.309039 | 0.000103
8.457421 | -7.879375 | 0.300609 | 0.000104
8.460504 | -7.866011 | 0.304762 | 0.000103
8.460571 | -7.821856 | 0.312767 | 0.000104
8.462892 | -7.842867 | 0.299458 | 0.000100
8.463775 | -7.837614 | 0.304938 | 0.000101
8.463446 | -7.852147 | 0.312048 | 0.000101
8.465396 | -7.800806 | 0.209615 | 0.000103
8.466696 | -7.805417 | 0.166386 | 0.000104
8.466279 | -7.910831 | 0.300395 | 0.000119
8.466696 | -7.805417 | 0.303377 | 0.000103
8.468108 | -7.882522 | 0.301556 | 0.000102
8.469058 | -7.866528 | 0.308831 | 0.000108
8.469108 | -7.846475 | 0.300255 | 0.000108
8.470600 | -7.894619 | 0.300445 | 0.000115
8.469867 | -7.922953 | 0.299295 | 0.000100
8.472408 | -7.810456 | 0.368080 | 0.000101
8.479588 | -7.966142 | 0.301866 | 0.000104
8.487733 | -7.926608 | 0.301159 | 0.000101
8.491658 | -7.896106 | 0.297784 | 0.000100
8.494708 | -7.913458 | 0.307990 | 0.000874
8.497242 | -7.883722 | 0.298694 | 0.000107
8.512950 | -7.938439 | 0.533748 | 0.000110
8.440079 | -7.972850 | 0.304762 | 0.000106
8.438658 | -7.847369 | 0.302260 | 0.000147
8.440954 | -7.950100 | 0.761086 | 0.000358
8.441692 | -7.832436 | 0.295155 | 0.000102
8.441692 | -7.832436 | 0.304628 | 0.000120
8.444413 | -7.918158 | 0.301403 | 0.000102
8.444925 | -7.742297 | 0.298120 | 0.000119
8.445521 | -7.845039 | 0.643521 | 0.000100
8.446550 | -7.870319 | 0.303561 | 0.000101
8.443504 | -7.981761 | 0.390419 | 0.000106
8.449642 | -7.842747 | 0.371022 | 0.000101
8.448617 | -7.785400 | 0.306663 | 0.000103
8.453133 | -7.920878 | 0.293771 | 0.000101
8.455708 | -7.876181 | 0.412010 | 0.000107
8.458538 | -7.838819 | 0.316319 | 0.000101
8.458929 | -7.864583 | 0.303911 | 0.000103
8.460292 | -7.783261 | 0.263697 | 0.000101
8.459621 | -7.840900 | 0.295088 | 0.000100
8.467704 | -7.995767 | 0.449141 | 0.000103
8.463604 | -7.836289 | 0.304228 | 0.000104
8.464392 | -7.744386 | 0.127610 | 0.000100
8.464413 | -7.885931 | 0.298040 | 0.000102
8.466088 | -7.898686 | 0.299508 | 0.000104
8.465567 | -7.866972 | 0.302907 | 0.000101
8.467179 | -7.830192 | 0.302136 | 0.000114
8.470292 | -7.732906 | 0.457597 | 0.000110
8.467813 | -7.806775 | 0.166486 | 0.000103
8.471492 | -7.765292 | 0.458958 | 0.000101
8.471050 | -7.773569 | 0.276543 | 0.000101
8.474892 | -7.873158 | 0.307967 | 0.000100
8.476692 | -7.907672 | 0.368621 | 0.000100
8.477433 | -7.955953 | 0.369154 | 0.000103
8.480146 | -7.883764 | 0.296739 | 0.000107
8.478713 | -7.778783 | 0.305796 | 0.000104
8.481079 | -7.796372 | 0.370439 | 0.000102
8.481037 | -7.852150 | 0.303194 | 0.000102
8.479967 | -7.927103 | 0.292623 | 0.000104
8.484304 | -7.771003 | 0.469204 | 0.000102
8.485038 | -7.902189 | 0.412787 | 0.000105
8.486742 | -7.826083 | 0.369608 | 0.000111
8.487963 | -7.821725 | 0.367643 | 0.000108
8.487963 | -7.821725 | 0.367550 | 0.000115
8.488279 | -7.849608 | 0.304134 | 0.000101
8.490088 | -7.816639 | 0.302260 | 0.000101
8.494679 | -7.892922 | 0.144521 | 0.000101
8.497133 | -7.756675 | 0.290352 | 0.000108
8.500467 | -7.915733 | 0.458747 | 0.000103
8.503554 | -7.811536 | 0.371723 | 0.000105

Comparison to archival redshifts. We searched NED for archival spectroscopic
redshifts within a radius of 5′, and found 16 galaxies, largely from BOSS
(Dawson et al., 2013). Of these, five were galaxies that we had targeted. The
mean redshift difference between independent measurements of the same target
is $9$ km/s, with an rms scatter of $9$ km/s. We then removed the duplicates
and merged the catalogs from NED and from our two masks to produce a final
catalog of 113 galaxies. Figure 2 shows a histogram of these redshifts.

Figure 2: Redshift histogram, with inset showing the redshift interval around
Abell 56.

###  3.2 Subclustering and kinematics

Non-Abell 56 structures. Figure 2 reveals a potential background cluster at
$z=0.37$, and possibly another at $z=0.46$. To assess how strongly clustered
these galaxies are on the sky, Figure 3 shows a sky map color-coded by
redshift. Galaxies in the putative cluster at $z=0.37$ (0.46) are shown as
open (closed) green circles, while galaxies near $z=0.30$ (i.e., associated
with Abell 56) are shown with a continuous color map that contains no green.
Neither set of background galaxies shows signs of clustering in space.
Furthermore, we estimate the velocity dispersion of each set using the
biweight estimator (Beers et al., 1990) and find only $279\pm 45$ ($105\pm
48$) for the structure at $z=0.37$ (0.46). Uncertainties on biweight
estimators are obtained by the jackknife method throughout this paper. These
velocity dispersions are far less than the velocity dispersion of Abell 56
(below), suggesting that they are an order of magnitude less massive and
unlikely to substantially contaminate the weak-lensing and X-ray maps
presented below.

Figure 3: Redshift map. Abell 56 galaxies are coded with a continuous color
map, while galaxies in the putative background cluster at $z=0.37$ (0.46) are
shown as open (closed) green circles. XMM-Newton contours are shown in red.

Abell 56. Of 67 galaxies in the $0.285\leq z\leq 0.314$ window, the biweight
estimate for the systemic redshift is $0.30256\pm 0.00058$. At this redshift,
the physical scale is 4.521 kpc/arcsec given our adopted cosmological model
(Wright, 2006). The redshift distribution is compatibile with a single
Gaussian, according to a Kolmogorov-Smirnov test. This is consistent with the
low line-of-sight velocity component $\Delta v_{\rm los}$ suggested by the
archival redshifts of the north and south BCGs. The biweight estimate of
velocity dispersion is $1264\pm 145$ km/s—rather large, but likely to be
inflated by merger activity as noted below.

We use the mc3gmm code (Golovich et al., 2019) to assign galaxy membership to
subclusters. This code models the distribution of galaxies in (RA, Dec, z)
space as a mixture of $N$ elliptical Gaussian profiles (i.e., subclusters),
with physically motivated priors on subcluster variance in each dimension as
well as covariance (i.e., ellipticity and rotation). $N$ is determined by the
user; we set $N{=}2$ based on the optical imaging and further supported by the
lensing data presented in §(1).555Golovich et al. (2019) varied $N$ and for
each merging system found the value of $N$ that best satisfied the Bayesian
Information Criterion (BIC), thus deriving $N$ from the spectroscopic data
alone. Here, the subcluster separation is smaller than typically seen in
Golovich et al. (2019), and the spectroscopic data points are fewer, making it
more difficult to meet BIC criteria for $N{>}1$ based on the spectroscopy
alone. The user sets nonoverlapping bounds for the central (RA, Dec, z) of
each subcluster to avoid degeneracies; mc3gmm maximizes the likelihood by
adjusting the parameters within those bounds. We run mc3gmm on the galaxies in
the redshift window $0.285\leq z\leq 0.314$ and the result is shown in Figure
4. The velocities of the subclusters are nearly identical, suggesting that the
relative motion of the subclusters is in a direction close to the plane of the
sky. The biweight estimate for the systemic redshift of the 33 (34) north
(south) members is $0.30298\pm 0.00099$ ($0.30231\pm 0.00071$). This yields
$\Delta v_{\rm los}=153\pm 281$ km/s.

The biweight velocity dispersion is $1283\pm 236$ km/s for the north
subcluster, and $1251\pm 191$ for the south. Simulations of merging clusters
(e.g., Pinkney et al., 1996; Takizawa et al., 2010) show that a pericenter
passage in the plane of the sky boosts the observed velocity dispersion by a
factor of ${\approx}1.5$ for hundreds of Myr afterward. Hence, one should not
interpret these large velocity dispersions as indicative of extremely massive
clusters.

Figure 4: Corner plot showing distribution of subcluster members in RA, DEC,
and velocity space relative to the overall mean. Members of the north (south)
subcluster are shown in blue (red).

##  4 Weak lensing analysis

We perform a weak-lensing analysis on the HST F814W imaging. Galaxies are
detected in the F814W image with SExtractor (Bertin & Arnouts, 1996). For each
galaxy, PSF models are generated following the method of Jee et al. (2007) by
utilizing their publicly available PSF catalog. Each galaxy is fit with a PSF-
convolved Gaussian distribution and the complex ellipticities are recorded
(our ACS weak-lensing pipeline is outlined in Finner et al., 2017, 2021,
2023). Objects with ellipticity greater than 0.8, ellipticity uncertainty
greater than 0.3, and intrinsic size (pre-psf) less than 0.5 pixels are
removed to prevent spurious sources such as diffraction spikes around bright
stars and poorly fit objects from entering the source catalog.

The next step is to eliminate as many foreground and cluster galaxies as
possible, while still retaining a sizeable sample of background sources. With
only single-band imaging available, we select galaxies with F814W AB
magnitudes fainter than 24. We apply this magnitude cut to the GOODS-S
photometric redshift catalog (Dahlen et al., 2013) and find that the
contamination by foreground galaxies is expected to be $\sim 2\%$. Cluster
galaxies may contribute additionally to the contamination. As their
contamination should be radially dependent, we test the radial dependence of
the source density. We find it to be flat, which suggests cluster galaxies are
not significantly contaminating our source catalog. The final source catalog
contains $\sim 43$ galaxies arcmin-2. The source catalog is then provided to
the FIATMAP code (Wittman et al., 2006) to create a surface mass density map.
FIATMAP convolves the observed shear field with a kernel of the form r^-2
(1-exp(-r^2 2ri2)) exp(-r^2 2ro2) where $r_{i}$ and $r_{o}$ are inner and
outer cutoffs, respectively. The inner cutoff is necessary to prevent
amplification of shape noise in sources at small $r$, and was set to 50
arcsec. The outer cutoff suppresses noise that may come from unrelated
structures along the line of sight at large projected separations, and is of
limited value in a small field; we set it to 100 arcsec, which is comparable
to the radius of the field. The results were pixelized onto a map with 1.5
arcsec pixels. In addition to this fiducial map, a family of viable
reconstructions can be made by bootstrap resampling the shear catalog (see
below). Figure 5 shows the fiducial map as a set of contours overlaid on a
Pan-STARRS multiband image666Retrieved from http://ps1images.stsci.edu/cgi-
bin/ps1cutouts. (Waters et al., 2020). Two weak lensing peaks are evident,
associated with (albeit slightly offset from) each galaxy subcluster, with the
X-ray peak in between. This confirms the basic merger scenario developed
above.

Figure 5: Surface mass density contours from weak lensing (white) overlaid on
a Pan-STARRS multiband image and red XMM-Newton surface brightness contours.
The small closed contour between the subclusters is a trough. Green GMRT 650
MHz contours (§6) start at 70 $\mu$Jy/beam with increments of 680 $\mu$Jy/beam
and a 4′′ beam.

We estimate the mass of each subcluster by fitting a two-halo NFW model with a
fixed mass-concentration relation from Diemer & Joyce (2019). To achieve the
best-fit model, the shear of the two-halo model is derived at the position of
each background galaxy and the chi-square is minimized. The effective distance
ratio of the sources is set by the effective distance ratio of GOODS-S sources
fainter than 24th magnitude. The best-fit two-halo model has a mass of
$M_{200}=4.5\pm 0.8\times 10^{14}$ M⊙ and $M_{200}=2.8\pm 0.7\times 10^{14}$
M⊙ for the south and north subclusters, respectively. We allow the centroid of
each halo to be fit and they converge to the projected mass distribution
peaks. On the other hand, if we fix the halo centroids to the BCGs, we find
the south and north subcluster masses decrease by $10\%$ and $60\%$,
respectively. To test the dependence of the mass estimate on our choice of
magnitude cut, we vary the magnitude constraint on the background catalog from
22nd to 25th magnitude and find that the mass estimates decrease for brighter
magnitude cuts but within the mass uncertainty. To estimate the total mass of
the cluster, we simulate two NFW halos at the projected separation of the two
mass peaks. Integrating the model from the center of mass to $R_{200}$, we
estimate the total mass of the cluster to be $M_{200}=9.7\pm 2.0\times
10^{14}$ M⊙ ($M_{500}=7.1\pm 1.6\times 10^{14}$ M⊙).

To further quantify the detection significance, we bootstrap resampled the
source catalog to generate 1000 realizations of the mass map. As expected, the
mean map yielded by these resamplings matches the fiducial map yielded by the
original catalog. At any given sky position, we can measure the rms variation
of surface mass density across the map realizations to obtain a noise map. The
ratio of the fiducial map to this noise map is then a significance map. The
peak of the southern (northern) subcluster is detected at a significance of
$6.3$ ($5.5$).

The projected separation between the mass peaks, $d_{\rm proj}=438$ kpc, is
important for the dynamical modeling in §5. To estimate the uncertainty on the
peak locations, the peak from each of the 1000 realizations was recorded. The
1000 peaks were then passed to a $k$-means algorithm with the number of
distributions fixed to two. The $k$-means algorithm iteratively calculates the
centroid of the peaks and assigns peaks to each centroid until the centroid
converges. This procedure yields two distributions of mass-peak locations,
which are then processed with a kernel density estimator to find the $1\sigma$
and $2\sigma$ uncertainties. We find that the southern mass peak is consistent
with its BCG at the $1\sigma$ level. In contrast, the northern mass peak is
offset $19.2\pm 4.9$ arcsec ($87\pm 22$ kpc) to the south of the northern BCG.
We address this offset further in §7. Our immediate goal here is to define a
68% confidence interval on the projected separation between mass peaks, which
we find to be $96.9\pm 45.6$ arcsec ($438\pm 206$ kpc).

To further check the halo position uncertainties, we consider again the two-
halo fit. As a model-driven procedure, this should be more robust against edge
effects than the mapping procedure, which convolves the observed shear field.
Nevertheless, as noted above, the halo center parameters converge to the
projected mass distribution peaks. The positional uncertainties from the two-
halo fit are smaller than those from the resampled mapping method. Hence our
adoption of the values from the latter method is the more cautious approach.

##  5 Simulated analogs and dynamical parameters

We find analog systems in the Big Multidark Planck (BigMDPL) Simulation
(Klypin et al., 2016) using the method of Wittman et al. (2018a) and Wittman
(2019). The observables used to constrain the likelihood of any given analog
and viewing angle are:

  * •

the projected separation between mass peaks $d_{\rm proj}$, for which we use
$438\pm 206$ kpc from §(1).

  * •

the line-of-sight relative velocity $\Delta v_{\rm los}$, for which we use
$153\pm 281$ km/s from §3.

  * •

the subcluster masses, for which we use $M_{200}=4.5\pm 0.8\times 10^{14}$ M⊙
and $M_{200}=2.8\pm 0.7\times 10^{14}$ M⊙ for the south and north subclusters,
respectively, from §(1). Note that dynamical timescales and velocities depend
only weakly on the masses.

Table 2 lists the resulting highest probability density confidence intervals
for time since pericenter (TSP), pericenter speed $v_{\rm max}$, viewing angle
$\theta$ (defined as the angle between the subcluster separation vector and
the line of sight, i.e. 90∘ when the separation vector is in the plane of the
sky), and the angle $\varphi$ between the current separation and velocity
vectors. $\varphi$ is potentially an indicator of how head-on the trajectory
is, as well as of merger phase (surpassing 90∘ at apocenter). The likelihood
ratio of analogs in the outbound vs. returning phase is 19:1. Table 2 also
lists the confidence intervals for the dynamical parameters when the analysis
is restricted to the outbound scenario. These particular parameters are not
sensitive to the current merger phase.

Table 2: Dynamical parameters from analogs Scenario | TSP (Myr) |  $v_{\rm max}$ (km/s) |  $\theta$ (deg) |  $\varphi$ (deg)
---|---|---|---|---
68% CI
All | 60-271 | 1960-2274 | 68-90 | 6-33
Outbound | 90-291 | 1952-2282 | 68-90 | 0-26
95% CI
All | 0-451 | 1729-2510 | 42-90 | 0-86
Outbound | 0-366 | 1681-2489 | 44-90 | 1-67

##  6 Radio observations and results

Pericenter speeds in cluster mergers are typically greater than the sound
speed in the gaseous ICM, so each subcluster launches a shock in the ICM of
the other subcluster (Ha et al., 2018). In hydrodynamic simulations of the
Bullet (Springel & Farrar, 2007) the shock begins at pericenter speed and
loses very little speed over time, while the corresponding subcluster falls
behind due to the gravity of the other subcluster. Our Abell 56 analogs do not
include gas, but we use the pericenter speed, gravitational subcluster
slowing, and analog time of observation to predict the separation between a
subcluster and a hypothetical constant-velocity shock. We find ${\sim}200$ kpc
separation in the outbound phase. The analogs indicate that an additional
${\approx}1.2$ Gyr passes before the subclusters return to the same projected
separation en route to a second pericenter. In this time, a hypothetical
constant-velocity shock would have proceeded over 2 Mpc further out.
Therefore, observing the shock location could further disambiguate between
outbound and returning scenarios. This toy model glosses over the complexities
of ICM properties affecting the shock speed, but the timescale of the
returning scenario is so long that the subcluster-shock separation remains
${>}1$ Mpc even with factor-of-two variations in shock speed, or complete
stalling of the shock after ${\sim}500$ Myr.

Shocks are often detected as discontinuities in the X-ray surface brightness,
but in this case the archival X-ray data are too shallow to support such a
detection. Shocks may also inject sufficient non-thermal energy into charged
particle motion that electrons emit synchrotron radiation, detectable as an
extended radio source known as a radio relic (van Weeren et al., 2019).
Archival 150 MHz data from the TIFR GMRT Sky Survey (TGSS) Alternative Data
Release (Intema et al., 2017) show extended emission 270 kpc south of the
southern BCG. Due to the large synthesized beam size (25′′) and an
accompanying point source, it is difficult to further characterize this
emission using the TGSS data alone. Cuciti et al. (2021) observed the cluster
at 1.5 GHz and $\approx 12^{\prime\prime}$ beam using the Jansky Very Large
Array (JVLA). The source in question appears at the southern edge of their
Figure A.1. However, they classified this cluster as having no extended
emission, presumably because they pointed at the original Abell coordinates,
about $7^{\prime}$ north of the source in question, and because they were
primarily searching for radio halos rather than relic candidates. We also
checked the VLASS (Lacy et al., 2020) and GLEAM (Wayth et al., 2015) surveys,
and found no evidence of a halo or relic.

We were granted 15 hours on the upgraded GMRT (uGMRT, Gupta et al., 2017) for
Band 4 (550-900 MHz) observations of Abell 56 (proposal code 42_069) with much
smaller synthesized beam size (4′′). Observations were taken on 20 June 2022
and 24 June 2022. We used the SPAM pipeline (Intema, 2014) to calibrate the
visibilities, and used wsclean (Offringa et al., 2014; Offringa & Smirnov,
2017) to create an image. The source 270 kpc south of the southern BCG extends
for $\approx 420$ kpc (93′′) in the east-west direction and is barely resolved
in the north-south direction. Its contours are overlaid in green on the Pan-
STARRS image in Figure 5. This makes it clear that the bright point source at
the western end of the radio emission is coincident with a galaxy; our
redshift survey confirms that this galaxy is in the cluster.

The most likely explanation for most of this emission is an AGN tail. Given
the orientation of this feature which matches that expected of a merger shock,
it is worth considering that AGN tails play a role in the formation of some
relics by providing seed electrons that are re-accelerated by the passage of a
shock (e.g., van Weeren et al., 2017). In such cases there is spectral aging
across the narrow axis of the tail in addition to the expected aging from head
to tail. Exploring this possibility would require high angular resolution
spectral maps. Finally, we note that there is no evidence of a relic much
further south as expected in the returning scenario, nor of a relic on the
north side of the north subclusters.

##  7 Dark matter cross section estimate

Spergel & Steinhardt (2000) first suggested that dark matter (DM) particles
may scatter off each other in a process distinct from the interactions with
standard model particles that are probed by direct detection experiments. The
cross section for such scattering is usually quoted in terms of
$\frac{\sigma_{\rm DM}}{m}$, the cross section per unit mass, because the mass
of the DM particle is unknown. Markevitch et al. (2004) laid out multiple
physical arguments for inferring this parameter, at least at a back-of-the-
envelope level, from merging cluster observations. Simulations (e.g., Randall
et al., 2008; Robertson et al., 2017) are required to properly interpret such
observations. However, as a first estimate to motivate deeper observations and
perhaps simulations of Abell 56, we present an initial back-of-the-envelope
estimate.

One physical argument is that momentum exchange will slow the DM halos
relative to the galaxies, resulting in a DM-galaxy offset. Markevitch et al.
(2004) developed an argument based on finding no significant offset: requiring
that the scattering depth be ${<1}$ leads to an upper limit on
$\frac{\sigma_{\rm DM}}{m}$. In this case, there is a significant offset in
the north, so we turn to the method of Harvey et al. (2014) and Harvey et al.
(2015), which uses the ratio of gas-galaxy and DM-galaxy offsets. This method
relies on an analogy between DM and the much more interactive gas, so it has
some limitations, but it also reduces some sources of observational
uncertainty. Foremost, it eliminates the assumption that the surface mass
density relevant to DM scattering—the volume density integrated along the
merger axis—equals the surface mass density we can measure, which is nearly
perpendicular to the merger axis. In fact clusters are triaxial (Harvey et
al., 2021) and align to some extent with their neighboring clusters (Joachimi
et al., 2015), hence one may expect greater column density along the merger
axis. To the extent this pattern is echoed by the gas, the gas analogy may
reduce this systematic error. Second, the gas analogy eliminates any
uncertainty due to viewing angle, as that angle applies equally to the gas-
galaxy and DM-galaxy separations.

The chief limitation of the gas analogy is that it breaks down over time. SIDM
simulations show that, given enough time, the galaxies within each subcluster
fall back to their associated DM—and continue oscillating (Kim et al., 2017).
Around the time of apocenter between subclusters, the DM-galaxy offset in each
subcluster has a sign opposite that predicted by the gas analogy. Hence, the
gas analogy should not be applied if the system is observed long after
pericenter. The analogs indicate that Abell 56 is observed much closer to
pericenter than apocenter, so the gas analogy is appropriate here for a first
estimate.

In the southern subcluster, the DM-BCG separation777All separations in this
paragraph are quoted after projecting them onto the merger axis, but we note
that the components perpendicular to the merger axis are generally negligible.
is $7\pm 16$ kpc and the gas-BCG separation is $111\pm 38$ kpc, yielding
$\frac{\sigma_{\rm DM}}{m}=0.35\pm 1.03$ cm2/g, consistent with zero. In the
northern subcluster, the DM-BCG separation is $87\pm 22$ kpc, while the gas-
BCG separation is unclear because it is difficult to identify a gas peak
specifically associated with the northern subcluster. To be conservative we
use the offset to the main gas peak, $424\pm 38$ kpc. This yields
$\frac{\sigma_{\rm DM}}{m}=1.43\pm 0.61$ cm2/g. Multiplying the two
likelihoods yields $\frac{\sigma_{\rm DM}}{m}=1.10\pm 0.64$ cm2/g.

We performed a few checks on the statistical significance of the offset in the
north. None of the 1000 bootstrap realizations of the convergence map in §(1)
placed the overall mass peak as far north as the northern BCG, and only three
of them placed a local mass peak (defined as a peak in the northern half of
the field) that far north.

We emphasize the tentative nature of the dark matter constraint. More work
will be needed to understand why the northern subcluster has a significant DM-
BCG offset while the south does not. Ground-based weak lensing, or more space-
based pointings, may be helpful to reduce any systematic uncertainties related
to the relative small footprint of the ACS data. Deeper imaging may reveal
strongly lensed sources that could lead to more precise mass models. X-ray or
radio confirmation of a shock position could further build confidence in the
merger scenario. Even without detection of a shock, deeper data on the overall
X-ray morphology combined with hydrodynamical simulations would greatly
advance understanding of this merger.

##  8 Summary and discussion

We have presented a new binary, dissociative merging galaxy cluster discovered
by cross-referencing archival X-ray data with locations of bimodal redMaPPer
clusters. The selection technique has the potential to be applied more widely,
as optical surveys continue to cover more area more deeply than ever before.
In particular, the southern sky may provide new targets via the 5000 deg2 Dark
Energy Survey (DES; Abbott et al., 2018) and eventually the deeper 20,000 deg2
Legacy Survey of Space and Time (LSST; LSST Science Collaborations et al.,
2009). Finding the rare merger through pointed X-ray followup of selected
candidates will require very careful selection. The forthcoming eROSITA X-ray
survey could enable more of a cross-correlation approach where candidates are
selected based on joint optical and X-ray properties.

This particular cluster promises to be useful for constraints on
$\frac{\sigma_{\rm DM}}{m}$, given that its merger axis is close to the plane
of the sky and its trajectory was sufficiently head-on to provide a
substantial separation between the gas peak and the main BCG. The lensing map
presented here is based on a single orbit of ACS time, and should be
supplemented with deeper and wider data to better understand why there is a
significant offset in the north but not in the south. Hydrodynamical
simulations could shed light on whether this could happen in a Cold Dark
Matter (CDM) scenario, perhaps with projection effects or other complications
not identified here. Such simulations should also be compared to deeper X-ray
maps to confirm that we understand the merger scenario.

To place this system in context with other merging clusters with the potential
to probe $\frac{\sigma_{\rm DM}}{m}$, we refer to Table 1 of Wittman et al.
(2018b), which ranked the importance of various subclusters used in their
ensemble analysis and that of Harvey et al. (2015). In Harvey et al. (2015),
the measurement uncertainties on the “star-gas” separation $\delta_{\rm SG}$
and the “star-interacting DM” separation $\delta_{\rm SI}$ were assumed to be
the same for all subclusters in the ensemble. Wittman et al. (2018b) noted
that this resulted in a particularly simple analytic expression for the
(unnormalized) inverse-variance weight of a given subcluster in an ensemble:
${\delta_{\rm SG}^{2}\over 1+\delta_{\rm SI}^{2}/\delta_{\rm SG}^{2}}$. By
glossing over the measurement uncertainties in any given observation, this
quantifies the importance of a subcluster in a hypothetical ensemble where all
subclusters are equally well observed. After normalizing this weight in the
same way as did Wittman et al. (2018b) for their Table 1, we find that the
southern subcluster of Abell 56 would appear in eighth place on the list of
usable subclusters (additional subclusters with formally greater weight were
marked as unusable in that table due to various complications). The northern
subcluster of Abell 56 is difficult to place on this table because only an
upper limit, not a measurement, is available for $\delta_{\rm SG}$. More X-ray
data will be needed to determine the constraining potential of this
substructure.

RJvW acknowledges support from the ERC Starting Grant ClusterWeb 804208. We
thank Nissim Kanekar for help with GMRT exposure time calculations, and Huib
Intema for help with the SPAM pipeline. We thank the staff of the GMRT that
made these observations possible. GMRT is run by the National Centre for Radio
Astrophysics of the Tata Institute of Fundamental Research. Some of the data
presented in this paper were obtained from the Mikulski Archive for Space
Telescopes (MAST) at the Space Telescope Science Institute. The specific
observations analyzed can be accessed via https://doi.org/10.17909/d922-3v11
(catalog 10.17909/d922-3v11).

## References

  * Abbott et al. (2018) Abbott, T. M. C., Abdalla, F. B., Allam, S., et al. 2018, ApJS, 239, 18
  * Abell et al. (1989) Abell, G. O., Corwin, Harold G., J., & Olowin, R. P. 1989, ApJS, 70, 1
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Arnaud et al. (2002) Arnaud, M., Majerowicz, S., Lumb, D., et al. 2002, Astronomy & Astrophysics, 390, 27
  * Beck et al. (2021) Beck, R., Szapudi, I., Flewelling, H., et al. 2021, MNRAS, 500, 1633
  * Beers et al. (1990) Beers, T. C., Flynn, K., & Gebhardt, K. 1990, AJ, 100, 32
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
  * Bradač et al. (2006) Bradač, M., Clowe, D., Gonzalez, A. H., et al. 2006, ApJ, 652, 937
  * Carter & Read (2007) Carter, J. A., & Read, A. M. 2007, Astronomy & Astrophysics, 464, 1155
  * Clowe et al. (2006) Clowe, D., Bradač, M., Gonzalez, A. H., et al. 2006, ApJ, 648, L109
  * Cuciti et al. (2021) Cuciti, V., Cassano, R., Brunetti, G., et al. 2021, A&A, 647, A50
  * Dahlen et al. (2013) Dahlen, T., Mobasher, B., Faber, S. M., et al. 2013, ApJ, 775, 93
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10
  * Diemer & Joyce (2019) Diemer, B., & Joyce, M. 2019, ApJ, 871, 168
  * Faber et al. (2003) Faber, S. M., Phillips, A. C., Kibrick, R. I., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1657–1669
  * Finner et al. (2017) Finner, K., Jee, M. J., Golovich, N., et al. 2017, ApJ, 851, 46
  * Finner et al. (2021) Finner, K., HyeongHan, K., Jee, M. J., et al. 2021, ApJ, 918, 72
  * Finner et al. (2023) Finner, K., Randall, S. W., Jee, M. J., et al. 2023, ApJ, 942, 23
  * Golovich et al. (2019) Golovich, N., Dawson, W. A., Wittman, D. M., et al. 2019, ApJ, 882, 69
  * Gupta et al. (2017) Gupta, Y., Ajithkumar, B., Kale, H. S., et al. 2017, Current Science, 113, 707
  * Ha et al. (2018) Ha, J.-H., Ryu, D., & Kang, H. 2018, ApJ, 857, 26
  * Harvey et al. (2015) Harvey, D., Massey, R., Kitching, T., Taylor, A., & Tittley, E. 2015, Science, 347, 1462
  * Harvey et al. (2021) Harvey, D., Robertson, A., Tam, S.-I., et al. 2021, MNRAS, 500, 2627
  * Harvey et al. (2014) Harvey, D., Tittley, E., Massey, R., et al. 2014, MNRAS, 441, 404
  * Intema (2014) Intema, H. T. 2014, SPAM: Source Peeling and Atmospheric Modeling, Astrophysics Source Code Library, record ascl:1408.006, , , ascl:1408.006
  * Intema et al. (2017) Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A, 598, A78
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Jee et al. (2007) Jee, M. J., Blakeslee, J. P., Sirianni, M., et al. 2007, PASP, 119, 1403
  * Joachimi et al. (2015) Joachimi, B., Cacciato, M., Kitching, T. D., et al. 2015, Space Sci. Rev., 193, 1
  * Kim et al. (2017) Kim, S. Y., Peter, A. H. G., & Wittman, D. 2017, MNRAS, 469, 1414
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340
  * Lacy et al. (2020) Lacy, M., Baum, S. A., Chandler, C. J., et al. 2020, PASP, 132, 035001
  * LSST Science Collaborations et al. (2009) LSST Science Collaborations, Abell, P. A., Allison, J., et al. 2009, ArXiv e-prints, arXiv:0912.0201
  * Markevitch et al. (2004) Markevitch, M., Gonzalez, A. H., Clowe, D., et al. 2004, ApJ, 606, 819
  * NASA/IPAC Extragalactic Database (2019) (NED) NASA/IPAC Extragalactic Database (NED). 2019, NASA/IPAC Extragalactic Database (NED), IPAC, doi:10.26132/NED1
  * Newman et al. (2013) Newman, J. A., Cooper, M. C., Davis, M., et al. 2013, ApJS, 208, 5
  * Offringa et al. (2014) Offringa, A. R., McKinley, B., Hurley-Walker, et al. 2014, MNRAS, 444, 606
  * Offringa & Smirnov (2017) Offringa, A. R., & Smirnov, O. 2017, MNRAS, 471, 301
  * Pinkney et al. (1996) Pinkney, J., Roettiger, K., Burns, J. O., & Bird, C. M. 1996, ApJS, 104, 1
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A27
  * Prochaska et al. (2020) Prochaska, J. X., Hennawi, J. F., Westfall, K. B., et al. 2020, Journal of Open Source Software, 5, 2308
  * Prochaska et al. (2020) Prochaska, J. X., Hennawi, J., Cooke, R., et al. 2020, pypeit/PypeIt: Release 1.0.0, v.v1.0.0, Zenodo, doi:10.5281/zenodo.3743493
  * Racca et al. (2016) Racca, G. D., Laureijs, R., Stagnaro, L., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9904, Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, ed. H. A. MacEwen, G. G. Fazio, M. Lystrup, N. Batalha, N. Siegler, & E. C. Tong, 99040O
  * Randall et al. (2008) Randall, S. W., Markevitch, M., Clowe, D., Gonzalez, A. H., & Bradač, M. 2008, ApJ, 679, 1173
  * Repp & Ebeling (2018) Repp, A., & Ebeling, H. 2018, MNRAS, 479, 844
  * Robertson et al. (2017) Robertson, A., Massey, R., & Eke, V. 2017, MNRAS, 465, 569
  * Rozo & Rykoff (2014) Rozo, E., & Rykoff, E. S. 2014, ApJ, 783, 80
  * Rykoff et al. (2014) Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104
  * Sereno & Ettori (2017) Sereno, M., & Ettori, S. 2017, MNRAS, 468, 3322
  * Simet et al. (2017) Simet, M., McClintock, T., Mandelbaum, R., et al. 2017, MNRAS, 466, 3103
  * Snowden & Kuntz (2014) Snowden, S. L., & Kuntz, K. D. 2014, XMM ESAS Cookbook, ,
  * Spergel & Steinhardt (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760
  * Springel & Farrar (2007) Springel, V., & Farrar, G. R. 2007, MNRAS, 380, 911
  * Takizawa et al. (2010) Takizawa, M., Nagino, R., & Matsushita, K. 2010, PASJ, 62, 951
  * van Weeren et al. (2019) van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, Space Sci. Rev., 215, 16
  * van Weeren et al. (2017) van Weeren, R. J., Andrade-Santos, F., Dawson, W. A., et al. 2017, Nature Astronomy, 1, 0005
  * Waters et al. (2020) Waters, C. Z., Magnier, E. A., Price, P. A., et al. 2020, ApJS, 251, 4
  * Wayth et al. (2015) Wayth, R. B., Lenc, E., Bell, M. E., et al. 2015, PASA, 32, e025
  * Wittman (2019) Wittman, D. 2019, ApJ, 881, 121
  * Wittman et al. (2018a) Wittman, D., Cornell, B. H., & Nguyen, J. 2018a, ApJ, 862, 160
  * Wittman et al. (2006) Wittman, D., Dell’Antonio, I. P., Hughes, J. P., et al. 2006, ApJ, 643, 128
  * Wittman et al. (2018b) Wittman, D., Golovich, N., & Dawson, W. A. 2018b, ApJ, 869, 104
  * Wright (2006) Wright, E. L. 2006, PASP, 118, 1711
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579


# Organic photovoltaic bulk heterojunctions with spatially varying
composition.

Paul M. Haney1 1Center for Nanoscale Science and Technology, National
Institute of Standards and Technology, Gaithersburg, Maryland 20899-6202, USA

###### Abstract

Models of organic bulk heterojunction photovoltaics which include the effect
of spatially varying composition of donor/acceptor materials are developed and
analyzed. Analytic expressions for the current-voltage relation in simplified
cases show that the effect of varying blend composition on charge transport is
minimal. Numerical results for various blend compositions, including the
experimentally relevant composition of a donor-rich region near the cathode (a
“skin layer” of donor material), show that the primary effect of this
variation on device performance derives from its effect on photocharge
generation. The general relation between the geometry of the blend and its
effect on performance is given explicitly. The analysis shows that the effect
of a skin layer on device performance is small.

##  I introduction

Photovoltaic devices consisting of two types of organic materials (referred to
as donor (D) and acceptor (A)) have attracted considerable scientific interest
in recent years. Their operation consists of the generation of an exciton in
the donor molecule, which is then disassociated into free carriers at the D-A
interface (the electron is transferred to the acceptor molecule’s lowest
unoccupied molecular orbital (LUMO), leaving a hole in the donor molecule’s
highest occupied molecular orbital (HOMO)). Carriers which avoid recombination
are then collected by contacts. The geometry first studied consisted of single
D and A layers, with a single planar interface tang . The resulting
efficiencies were low (1 %), owing in part to the short exciton diffusion
length (10 nm) - only excitons within this short distance from the interface
lead to free carriers. It was subsequently discovered that blending D and A
together throughout the device thickness led to increased efficiencies yu ,
now above 5 % eff1 ; eff2 ; eff3 . This increase in efficiency is attributed
to an increase in D-A interfacial area; carrier transport is sufficiently
robust to the disorder present in the blend to accommodate reasonable quantum
efficiencies. If the organic blend is completely homogeneous, the contacts on
the device must be different in order to break spatial symmetry and permit a
nonzero short-circuit current in a preferred direction. The key difference
between the contacts is their work function: a lower (higher) work function
ensures that the contact preferentially collects and injects electrons
(holes). Hence it is understood that the cathode collects electrons, and the
anode collects holes.

A major thrust of experimental efforts has been to attain control over blend
morphology in order to optimize both exciton disassociation and charge
transport. Recent examples include using nanoimprint lithography to control
the structure of the donor-acceptor molecules’ interfacial profile morph1 , or
using a graded donor-acceptor blend to optimize both carrier collection and
transport holmes . A key challenge of engineering blend morphology is the
measurement and characterization of the structure of the organic blend.
Techniques for accomplishing this include atomic force microscopy hamadani ,
ellipsometry germack , and X-ray photoelectron spectroscopy xu . These
techniques have revealed that typical methods for fabricating devices lead to
a layer of enhanced donor molecule density at the cathode, which has been
attributed to surface energy differences between the active layer and other
components xu . This would seem to present an impediment to good device
performance: the cathode collects electrons, but in its vicinity is mostly
holes! Nevertheless, internal quantum efficiencies of 90 % have been observed
in these materials schilinksy , indicating that charge collection is still a
relatively efficient process germack .

In this work, I theoretically study the effect of a nonuniform blend on
organic photovoltaic (OPV) device performance. I employ a drift-diffusion
equation to describe electron and hole transport, a field and temperature
dependent generation and recombination rate that captures the exciton physics,
and the Poisson equation for electrostatics. To this model I add the effect of
a spatially varying effective density of states (EDOS) (note that “density of
states” refers to the number of states per volume per unit energy, whereas
“effective density of states” refers to the number of states per volume). Part
I describes details of the model. In part II, I present analytic solutions for
the transport under certain approximations; these point to the fact that the
effect of a spatially varying EDOS on charge transport is small. In part III,
I present numerical results which indicate that the primary effect of a
spatially varying EDOS is on the charge generation and ensuing $J_{\rm sc}$.
It is shown that this can be understood in terms of the overall geometry of
the composition profile. I conclude that since the skin layer near the cathode
is geometrically small on the scale of the device thickness, its effect on
performance is similarly small.

##  II Model

The model used to describe the system is similar to that found in Ref. koster,
. Its basic equations are presented here in dimensionless form. Table I shows
the variable scalings used. The dimensionless drift-diffusion/Poisson
equations including a spatially varying EDOS are given as fonash :

$\displaystyle J_{n}$ $\displaystyle=$ $\displaystyle
n\left(-\frac{\partial}{\partial x}V-\frac{\partial}{\partial x}\left[{\rm
ln}N\right]\right)+\frac{\partial}{\partial x}n~{},$ $\displaystyle J_{p}$
$\displaystyle=$ $\displaystyle f_{\mu}\left[p\left(-\frac{\partial}{\partial
x}V+\frac{\partial}{\partial x}\left[{\rm
ln}P\right]\right)-\frac{\partial}{\partial x}p\right]~{},$
$\displaystyle-\frac{\partial}{\partial x}J_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{f_{\mu}}\frac{\partial}{\partial x}J_{p}=G-R~{},$ (1)
$\displaystyle\frac{\partial^{2}}{\partial x^{2}}V$ $\displaystyle=$
$\displaystyle p-n~{},$ (2)

where $f_{\mu}=\mu_{h}/\mu_{e}$ is the ratio of hole to electron mobility, $G$
is the carrier density generation rate, and $R$ is the recombination. $N(x)$
and $P(x)$ are the spatially-dependent electron and hole effective density of
states, respectively. $n$ and $N$ are related via:
$n=Ne^{-\left(E_{c}-E_{F,n}\right)/kT}$, where $E_{F,n}$ is the electron
quasi-Fermi level, $E_{c}$ is the conduction band edge, and all quantities are
position-dependent (the densities are assumed to be such that the system is in
a nondegenerate regime). $p$ and $P$ are related similarly. $N(x)$ and $P(x)$
are fixed material parameters, while $n$ and $p$ are system variables that
depend on applied voltage and illumination. For a single band semiconductor,
the effective density of states $N$ is given by
$\frac{1}{\sqrt{2}}\left(\frac{m^{*}_{n}k_{\rm
B}T}{\pi\hbar^{2}}\right)^{3/2}$, where $m^{*}_{n}$ is the effective electron
mass. In the present context of organic materials, $N$ is more properly
understood as the number of HOMO states per unit volume, and is proportional
to the donor molecule density.

Table 1: Normalization to dimensionless variables. In the below $N_{0}$ is the characteristic density (typically chosen to be on the order of $10^{-25}~{}{\rm m}^{-3}$) $D_{n}$ is the electron diffusivity, $\epsilon$ is the dielectric constant of the organic blend, $q$ is the magnitude of the electron charge, $T$ is the temperature, and $k_{\rm B}$ is Boltzmann’s constant. Quantity | Normalization
---|---
density | $N_{0}$
position | $\sqrt{\epsilon k_{B}T/(q^{2}N_{0})}\equiv x_{0}$
charge current | $qD_{n}N_{0}/x_{0}$
electric potential | $k_{\rm B}T/q$
rate density | $x_{0}^{2}/N_{0}D_{n}$

The boundary conditions are given as:

$\displaystyle n\left(0\right)$ $\displaystyle=$ $\displaystyle
N\left(0\right)e^{-E_{g}+\phi_{L}},$ $\displaystyle p\left(0\right)$
$\displaystyle=$ $\displaystyle P\left(0\right)e^{-\phi_{L}},$
$\displaystyle~{}n\left(L\right)$ $\displaystyle=$ $\displaystyle
N\left(L\right)e^{-\phi_{R}},$ $\displaystyle~{}p\left(L\right)$
$\displaystyle=$ $\displaystyle P\left(L\right)e^{-E_{g}+\phi_{R}},$ (3)

where $L$ is the device thickness (this represents placing the anode at $x=0$,
and the cathode at $x=L$). $\phi_{L(R)}$ is the absolute value of the
difference between HOMO (LUMO) and left (right) contact Fermi level. The
boundary condition for the Poisson equation is:

$\displaystyle V(L)-V(0)=\left(E_{g}-\phi_{L}-\phi_{R}\right)-V_{\rm A},$ (4)

where $V_{\rm A}$ is the applied voltage (with the sign convention above,
$V_{\rm A}>0$ corresponds to forward bias).

I consider only bimolecular recombination, with (dimensionless) form:

$\displaystyle R=\left(np-n_{i}^{2}\right)~{},$ (5)

where $n_{i}^{2}=n_{0}p_{0}$, and $n_{0}~{}(p_{0})$ is the equilibrium
electron (hole) density. The carrier generation rate density is taken to be
spatially uniform. As described in Ref. koster , adding the exciton density as
a system variable modifies the source term in Eq. (1):

$\displaystyle\left(G-R\right)\rightarrow\tilde{P}\times
G_{0}-\left(1-\tilde{P}\right)\times R,$ (6)

where $G_{0}$ is the exciton density generation rate, and $\tilde{P}$ is a
field and temperature dependent factor which represents the probability for an
exciton to disassociate into free electron and hole koster ; braun . The field
and temperature dependence is described by Braun’s extension of Onsager’s
theory of ion disassociation in electrolytes onsager .

Charge recombination and generation also generally depend on the donor and
acceptor effective density of states. The total source term of Eq. (1)
(denoted here by $U(x)$) is therefore of the generic form:

$\displaystyle U(x)$ $\displaystyle=$ $\displaystyle\tilde{P}\times
G_{0}\times g\left[N(x),P(x)\right]-$ (7)
$\displaystyle~{}~{}\left(1-\tilde{P}\right)\times R\times
r\left[N(x),P(x)\right].$

The appropriate forms for $g$ and $r$ depend on several factors, such as the
dependence of the optical absorption and D-A interface area on relative D-A
composition.

Figure 1: (a) Energy diagram for device model; cartoon of particle flow
depicts dark current in forward bias. (b) Spatial dependence of EDOS: linear
variation (shaded region), and step-like change (dotted line) in both $N(x)$
and $P(x)$.

##  III Analytic cases

The set of equations described in Eq. (1) can be solved analytically for
limiting cases, which can provide some insight into the effect of a spatially
varying EDOS on the transport. Two cases are considered here: the first is an
exponentially varying EDOS (which can be extrapolated to a linearly varying
EDOS), and the second is an abrupt, step-like change in the EDOS. I present
both solutions first and discuss the physics they describe second.

In both cases the electric field $E$ is taken to be spatially uniform (so that
$V(x)=-Ex$), and recombination is ignored. I suppose further that $G$ is
constant, and independent of $N,P$ (that is, $g(N,P)=1$ in Eq. (7)). The
exponentially varying EDOS is parameterized as:

$\displaystyle N\left(x\right)$ $\displaystyle=$ $\displaystyle
P\left(x\right)=A_{0}~{}e^{ax/L},$ (8)

where $A_{0}=a/\left(e^{a}-1\right)$ ensures that the total number of states
is independent of $a$. Substituting the expressions for electron (hole)
current density $J_{n(p)}$ into the equation of continuity (Eq. (1)) results
in a second order differential equation for the electron (hole) density $n$
($p$). For the EDOS of Eq. (8), the resulting general solution is:

$\displaystyle n\left(x\right)$ $\displaystyle=$ $\displaystyle
c_{1}e^{\left(a-f\right)x}+c_{2}+\frac{Gx}{a-f}~{},$ $\displaystyle
p\left(x\right)$ $\displaystyle=$ $\displaystyle
c_{1}e^{\left(a+f\right)x}+c_{2}+\frac{Gx}{a+f},$ (9)

where $c_{1},~{}c_{2}$ are determined by the boundary conditions of Eq. (3).
From this solution the current density can be obtained directly.

I express the resulting current-voltage relation as a sum of dark current and
light current:

$\displaystyle J\left(V_{\rm A}\right)=J_{D}+GJ_{L}~{},$ (10)

Both light and dark currents are well described by expanding to lowest order
in the spatial variation of EDOS parameter $a$; I take $\phi_{L}=\phi_{R}=0$,
and express the applied voltage dependence in terms of $f=\left(E_{g}-qV_{\rm
A}\right)/k_{\rm B}T$. $f$ is bigger than 1 in the region of interest
footnote1 , leading to the further approximation that $\sinh f\approx\cosh
f\approx 1/2~{}e^{f}$. It’s useful to express current-voltage relation in
terms of that for a uniform EDOS and electric field:

$\displaystyle J_{D}^{0}$ $\displaystyle=$
$\displaystyle\frac{2f\left(e^{V_{\rm
A}}-1\right)}{L\left(e^{f}-1\right)e^{V_{\rm A}}}$ $\displaystyle J_{L}^{0}$
$\displaystyle=$ $\displaystyle
L\left(\frac{2}{f}-\coth\left(\frac{f}{2}\right)\right).$ (11)

The dark and light current for the exponentially varying profile of Eq. (8) is
then found to be:

$\displaystyle J_{D}^{\rm exp}$ $\displaystyle\approx$ $\displaystyle
J_{D}^{0}\left(1+a^{2}\left(\frac{1}{12}-\frac{1}{2f}\right)+O\left(a^{4}\right)+...\right)~{},$
$\displaystyle J_{L}^{\rm exp}$ $\displaystyle\approx$ $\displaystyle
J_{L}^{0}\left(1+a^{2}\left(-\frac{2}{f^{3}}+e^{-f}\right)+O\left(a^{4}\right)+...\right).$
(12)

I next consider a step function form of $N(x),P(x)$. I use the following form:

$\displaystyle
N\left(x\right)=P\left(x\right)=\left\\{\begin{array}[]{rl}1-a/2&\text{if
}x<L/2,\\\ 1+a/2&\text{if }x\geq L/2.\end{array}\right.$ (15)

The general solutions for each region ($x<L/2,~{}x>L/2$) are of the form given
by Eq. (9) with $a$=0. In addition to the boundary condition Eq. (3), there’s
an additional boundary condition for this EDOS of continuity of charge and
current density at $x=L/2$. Making the same approximations as above leads to
the following dark and light current:

$\displaystyle J_{D}^{\rm step}$ $\displaystyle\approx$ $\displaystyle
J_{D}^{0}\left(1-a^{2}e^{-f/2}+O\left(a^{4}\right)+...\right)~{},$
$\displaystyle J_{L}^{\rm step}$ $\displaystyle\approx$ $\displaystyle
J_{L}^{0}\left(1-\frac{a^{2}}{2}\frac{fe^{-f/2}}{f-2}+O\left(a^{4}\right)+...\right).$
(16)

A number of interesting and relevant features emerge from these solutions:
first, only even powers of $a$ appear in the expansions. This is a consequence
of the symmetry built in to the system: when $\phi_{L}=\phi_{R}$ and
$f_{\mu}=1$, electron particle transport from left to right is equal to hole
particle transport from right to left. In both the exponential and step-like
cases above, holes encounter an expansion in the EDOS along their transport
direction, which increases the hole current. Conversely, electrons encounter a
constriction, which decreases the electron current. To linear order (and all
odd orders) in the expasion/contraction parameter $a$, these effects cancel
each other so that the total charge current only appears with even powers of
$a$. If the electron/hole symmetry is broken, or the symmetry of the EDOS is
reduced (by shifting the step away from the center of the device), then odd
powers of $a$ are present (with prefactors whose magnitude reflects the degree
of symmetry breaking).

The other relevant feature of Eqs. (12) and (16) is the small magnitude of the
$a^{2}$ prefactor. Noting again that $f$ is generally larger than 1 for the
applied voltages of interest to solar cells, it’s clear by inspection that the
prefactors are much smaller than 1. This indicates that the effect on
transport of a spatial variation of the EDOS is quite weak.

Figure 2: Extinction of current when EDOS goes to zero. This is for the step-
like change in EDOS, for parameters $G=10^{-9},~{}V_{\rm A}=0.7$. Both
approximate and exact values are shown (where the approximate expression is
given by Eq. (16)). It is seen that the current decreases substantially only
when the EDOS is nearly zero (or when $a$ is nearly 2).

The intuitive picture that emerges from this analysis is that electrons and
holes can very easily “squeeze” through regions of reduced density. A natural
question concerns the way in which transport is ultimately “pinched off” by
letting the density vanish at a point in space. This is shown in Fig. (2),
which shows the current in the step-like structure as $a\rightarrow 2$. The
way in which the current vanishes is very steep; it is only at very small
values of EDOS at the cathode that the current drops appreciably (in this
limit, the approximation $a\ll 1$ used in deriving Eq. (16) is not satisfied,
hence the discrepancy between exact solution and Eq. (16)). However, for very
small values of HOMO and LUMO density in real systems, the model presented
here is likely not appropriate. This point is discussed more fully in the
conclusion.

##  IV numerical studies

I next consider the effect of spatially varying EDOS when the Poisson equation
for the electric potential and bimolecular recombination are included. Recall
that the dependence of the generation and recombination on EDOS of electron
$N$ and hole $P$ is described generically as:

$\displaystyle\tilde{P}\times G_{0}\times
g\left(N,P\right)-\left(1-\tilde{P}\right)\times R\times
r\left(N,P\right)~{}.$ (17)

I make the following ansatz for $g$ (the main conclusion can be formulated in
a way that’s independent of this specific choice for $g$):

$\displaystyle g\left(N,P\right)=P^{2}N~{}.$ (18)

This is motivated by the observation that the D-A interfacial area requires
both $P$ and $N$, hence $g$ has a factor of each; an extra factor of $P$ is
added since the exciton is initially generated in the donor. $r$ is taken
simply to be 1, since $R$ already has $N$ and $P$ dependence built in through
$n$ and $p$. Adding a factor of $P$ to the recombination (so that the
$N,P$-dependence of both generation and recombination is the same) has only a
weak effect on the results.

A range of composition profiles has been explored for the numeric evaluation
of device performance, and I present two representative examples here:

$\displaystyle N_{1}(x)=1-P_{1}(x)$ $\displaystyle=$
$\displaystyle\frac{a}{\left(e^{a}-1\right)}e^{ax/L}~{},$ (19) $\displaystyle
N_{2}(x)=1-P_{2}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+\left(1-2a\right)\tanh\left(\frac{x-x_{0}}{\lambda}\right)\right).$

Fig. (3) shows the $J-V$ curves for the $\left(N_{2},P_{2}\right)$ case (Eq.
(IV)) for the uniform profile ($a=1/2$), and a sharp S-shaped profile
($a=0.95$). Note that the effect of the EDOS profile on the short circuit
current $J_{\rm sc}$ is substantial, while the effect on open circuit voltage
$V_{\rm oc}$ is small.

Figure 3: Current density-Voltage relation for two spatial profiles of D-A
EDOS profiles. Blue dotted line is for uniform EDOS profile, red line is for
S-shaped EDOS profile, given by Eq. (IV)

The previous analysis can explain the relative insensitivity of $V_{\rm oc}$
to a nontrivial EDOS profile: the effect of a varying EDOS profile on
transport is weak, so that the injected current required to offset the
photogenerated current (and the corresponding required voltage - $V_{\rm oc}$)
is only weakly sensitive to changes in EDOS footnote2 .

The change in $J_{\rm sc}$ can be understood as a direct consequence of the
model construction. $J_{\rm sc}$ is the current collected in the absence of an
applied voltage, that is, in the absence of charge injected from the contacts.
As such it is simply equal to the total charge generation rate in the device:
$J_{\rm sc}=\int dx\left({\rm Generation}(x)-{\rm Recombination}(x)\right)$.
As described above, this is directly parameterized as:

$\displaystyle J_{\rm sc}$ $\displaystyle=$ $\displaystyle\int
dx\left(\tilde{P}\times G_{0}\times g\left[N(x),P(x)\right]-\right.$ (21)
$\displaystyle\left.~{}~{}~{}~{}~{}~{}~{}~{}(1-\tilde{P})\times R\times
r\left[N(x),P(x)\right]\right).$

In analyzing the effect of $N(x),P(x)$ on $J_{\rm sc}$, it is instructive to
separate the $N,P$ dependence of the generation from the above integral. This
leaves a quantity $\delta U$ which depends only on the geometry of the D-A
EDOS profile:

$\displaystyle\delta U$ $\displaystyle=$ $\displaystyle\int
dx~{}g\left[N(x),P(x)\right]$ (22) $\displaystyle=$ $\displaystyle\int
dx~{}N(x)P^{2}(x).$ (23)

Strictly speaking the integral in Eq. (21) does not factorize in a manner that
leads lead directly to a $\delta U$ term. However, as I show in the following,
$\delta U$ is a good predictor of the effect of the geometry of the EDOS on
the device performance.

For each EDOS profile, I also vary other system parameters. The three
different parameterizations are shown in Fig. (4). In system 1, HOMO/LUMO
levels are aligned with cathode/anode Fermi levels ($\phi_{L}=\phi_{R}=0$),
and electron and hole mobility are equal. For system 2, $\phi_{L}=\phi_{R}=0$,
but electron and hole mobilities are not equal ($\mu_{e}=10\mu_{h}$). In
system 3, the HOMO/LUMO are offset from cathode/anode by $0.2~{}{\rm eV}$
($\phi_{L}=\phi_{R}=0.2~{}{\rm eV}$), and electron/hole mobilities are equal.

Figure 4: Cartoon of the three system parameterizations: system 1:
$\phi_{L}=\phi_{R}=0,~{}\mu_{h}=\mu_{e}$, system 2:
$\phi_{L}=\phi_{R}=0,~{}\mu_{h}=10\mu_{e}$, system 3:
$\phi_{L}=\phi_{R}=0.2~{}{\rm eV},~{}\mu_{h}=\mu_{e}$

Fig. (5a) and (5b) shows $\delta U$ as the profile parameter $a$ is varied,
for various EDOS configurations given by Eqs. (19) and (IV), respectively.
This is shown for the three system parameterizations. The overall device
efficiency $\eta$ tracks $\delta U$ very closely for all of these cases (the
efficiency is proportional to the maximum absolute value of $(JV)$ in the 4th
quadrant of the $J-V$ plane). For this reason I conclude that the primary
effect of a spatially varying EDOS on device performance is to change the
total carrier generation rate and ensuing $J_{\rm sc}$. $\delta U$ in Fig. (5)
is calculated using Eq. (23), however the conclusion is valid for any choice
of $g$ I’ve tried. Hence the effect of a nonuniform blend on performance can
be approximately specified in the generic form given by Eq. (22).

Figure 5: The efficiency $\eta$ and geometrical factor $\delta U$ (normalized
by their maximum value) versus geometrical parameter $a$ for (a) exponentially
varying profile ($N_{1}(x)$ of Eq. (19)) (b) S-shaped profile $N_{2}(x)$ of
Eq. (IV), with $x_{0}=L/2$, $\lambda=L/8$), (c) “skin” layer geometry.
Representations of the spatial variation of EDOS as a function of $a$ are
shown above the figure. The gray and white regions represent $N(x)$ and
$P(x)$, respectively. The efficiency closely follows the geometrical factor
$\delta U$ for most cases. For each geometry I use the three system
parameterizations described in Fig. (4) (the subscript of $\eta$ specifies the
system parameterization).

Next I turn to the experimentally motivated geometry of a skin layer of D near
the cathode. It’s parameterized as:

$\displaystyle
N(x)=1-P(x)=\frac{2+a}{4}+\frac{2-a}{4}\tanh\left(\frac{x-x_{0}}{\lambda}\right),$

with $\lambda=0.0075~{}L$, $x_{0}=0.05~{}L$. Fig. (5c) shows how the
efficiency evolves as the skin layer goes from mostly D-like (small $a$), to
an even D-A mix, to mostly A-like (large $a$) (the experimentally realistic
case is smaller $a$). The change in efficiency is a rather small effect for
all three cases (a maximum of 10 % change). Also shown is the geometrical
factor $\delta U$ (solid line). The efficiency of system 1 conforms most
closely to the geometrical factor profile dependence. Inspection of the $J-V$
curves for the three systems reveals subtle differences in the fill-factor
between the three; there is no simple or obvious source for the difference in
behavior between the three system parameterizations. The difference in
behavior between the three systems is more conspicuous for the skin layer
geometry because the effect of a nonuniform blend is smaller for the skin
layer, so that the overall performance is more sensitive to other system
parameters. (When the blend profile leads to larger effects, for example that
shown in Fig. (5b), there is a similar dependence of device performance on
profile for all system parameterizations.) Nevertheless, the important
conclusion common to all three system parameterizations of the skin layer
geometry is that the effect of the skin layer is small. Its smallness can be
understood in terms of the analysis of the previous sections. The analytic
work points to the fact that the effect of blend non-uniformity on charge
transport is generically small (except in extreme cases). The numerical work
of the previous test cases indicates that the effect of blend non-uniformity
can be understood in terms of its effect on charge generation and resulting
$J_{\rm sc}$ \- and that this effect is essentially geometrical (see Eq.
(22)). Since a skin layer is by definition geometrically small, its effect is
similarly small.

##  V Conclusion

In this work I presented a simple model for the effect of nonuniform blend
profiles on OPV device performance. The main effect of a nonuniform D-A blend
is on the the charge generation and resulting short-circuit current: in
regions where the blend is primarily of one type at the expense of the other,
there is less charge generation due to a reduced D-A interfacial area. The
details of how charge generation depends on local blend mix are complicated,
and involve almost all aspects of OPV device operation (e.g optics moule ,
exciton diffusion holmes , etc.). The influence of a nonuniform blend on
electron and hole transport is a weaker effect.

It’s important to appreciate the simplicity of the model presented here
relative to the complexity of real OPV devices. Two simplifications of the
model are: its treatment of the metal-organic interface, and its restriction
to 1 spatial dimension. I make no attempt to capture the effect of a skin
layer geometry on the metal-organic contact. The physics at this interface is
included most simply as a finite recombination velocity scott  (which can also
depend on temperature and field lacic ). A hallmark of less effective charge
collection/injection at this interface is S-shaped $J-V$ curves deibel . This
feature is correlated to metal contact deposition techniques deibel , and is
not ubiquitously observed in devices. I therefore conclude that the details of
the metal-organic contact is not directly tied to the phase segregation in the
organic blend.

A more severe approximation of this model is its restriction to 1-d. When the
EDOS is small, the charge and current density is also small. However,
experiments reveal localized hot-spots of conducting paths hamadani . A 1-d
model necessarily averages these localized hot-spots or large current density
over the entire cross-sectional area, leading to a diffuse current. As the
overall area of hot-spots decreases, the charge and current density they must
accommodate increases, and current may become space-charge limited. A 1-d
model is unable to capture the physics described in this scenario. However for
less extreme cases, the treatment described here offers the simplest account
for a spatially varying blend structure.

I acknowledge very useful discussions with Behrang Hamadani and Lee Richter.

## References

  * (1) C. W. Tang. Two-layer organic photovoltaic cell. Appl. Phys. Lett. 48, 183 (1986).
  * (2) G. Yu, J. Gao, J. C. Hummelen, F. Wudl, and A. J. Heeger, Science 270, 1789 (1995).
  * (3) Jiangeng Xue, Barry P. Rand, Soichi Uchida, and Stephen R. Forrest. J. Appl. Phys., 9, 124903 (2005).
  * (4) J. Peet, J. Y. Kim, N. E. Coates, W. L. Ma, D. Moses, A. J. Heeger, and G. C. Bazan. Nat. Mater. 6 497 (2007).
  * (5) S. H. Park, A. Roy, S. Beaupre, S. Cho, N. Coates, J. S. Moon, D. Moses, M. Leclerc, K. Lee, and A. J. Heeger. Nat. Photon. 3 297 (2009).
  * (6) D. Cheyns, K. Vasseur, C. Rolin, J. Genoe, J. Poortmans, and P. Heremans. Nanotechn. 19, 424016 (2008).
  * (7) 1 R. Pandey and R.J. Holmes, Adv. Mater. 22, 5301-5305 (2010).
  * (8) B.H. Hamadani, S. Jung, P. M. Haney, L. J. Richter, and N. B. Zhitenev, N.B. Nano letters 10, 1611 1617 (2010).
  * (9) D.S. Germack, C.K. Chan, R.J. Kline, D.A. Fischer, D.J. Gundlach, M.F. Toney, L.J. Richter, and D.M. DeLongchamp, Macromolecules 43, 3828 (2010).
  * (10) Z. Xu, L. Chen, G. Yang, C. Huang, J. Hou, Y. Wu, G. Li, C. Hsu, and Y. Yang, Advanced Functional Materials, 19, 1227 (2009).
  * (11) P. Schilinsky, C. Waldauf, and C.J. Brabec, App. Phys. Lett. 81, 3885 (2002).
  * (12) S. J. Fonash, Solar Cell Device Physics (Academic Press, Inc., London, 1981).
  * (13) B. V. Andersson, A. Herland, S. Masich, and Olle Ingan as. Nano Lett. 9 853, (2009).
  * (14) L. J. A. Koster, E. C. P. Smits, V. D. Mihailetchi, and P. W. M. Blom, Phys. Rev. B 72, 085205 (2005).
  * (15) The most relevant region of the current-voltage relation is in the 4th quadrant. In this region, $V_{\rm A}<E_{g}$; the different between $E_{g}$ and $qV_{\rm A}$ is scaled by $1/k_{\rm B}T=40\left({\rm eV}\right)^{-1}$, so that $f$ is generally much larger than 1.
  * (16) The change in $J_{\rm sc}$ induced by a spatially varying EDOS has some effect on $V_{\rm oc}$ as well, but this is also a small effect, as $V_{\rm oc}$ generically varies only logarithmically with $J_{\rm sc}$.
  * (17) H. K. Gummel, IEEE Transactions on Electron Devices, 11, 455 (1964).
  * (18) R. Sokel and R. C. Hughes, J. Appl. Phys, 53, 7414 (1982).
  * (19) Adam J. Moulé, Jörg B. Bonekamp, and Klaus Meerholz, J. Appl. Phys. 100, 094503 (2006).
  * (20) L. Braun, J. Chem. Phys. 80, 4157 (1984).
  * (21) L. Onsager, J. Chem. Phys. 2, 599 (1934).
  * (22) J.C. Scott and G.G. Malliaras, Chem. Phys. Lett. 299, 115-119 (1999).
  * (23) S. Lacic and O. Ingana?s, J. Appl. Phys. 97, 124901 (2005).
  * (24) A. Wagenpfahl, D. Rauh, M. Binder, C. Deibel, and V. Dyakonov, Phys. Rev. B 82, 115306 (2010).