# Accurate dust temperature determination in a $z=7.13$ galaxy

Tom J. L. C. Bakx1,2, Laura Sommovigo3, Stefano Carniani3, Andrea
Ferrara3,Hollis B. Akins4, Seiji Fujimoto5,6, Masato Hagimoto1, Kirsten K.
Knudsen7,Andrea Pallottini3, Yoichi Tamura1 and Darach Watson5,6
1 Division of Particle and Astrophysical Science, Graduate School of Science,
Nagoya University, Aichi 464-8602, Japan.
2 National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo
181-8588, Japan.
3 Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126, Pisa, Italy
4 Department of Physics, Grinnell College, 1116 Eighth Ave., Grinnell, IA
50112, USA
5 Cosmic Dawn Center (DAWN), Copenhagen, Denmark
6 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, Copenhagen,
Denmark
7 Department of Earth and Space Sciences, Chalmers University of Technology,
Onsala Space Observatory, SE-43992 Onsala, Sweden
E-mail<EMAIL_ADDRESS>(Nagoya University)

(Accepted 2021 September 3. Received 2021 September 3; in original form 2021
August 10.)

###### Abstract

We report ALMA Band 9 continuum observations of the normal, dusty star-forming
galaxy A1689-zD1 at $z=7.13$, resulting in a $\sim$4.6 $\sigma$ detection at
$702$ GHz. For the first time these observations probe the far infrared (FIR)
spectrum shortward of the emission peak of a galaxy in the Epoch of
Reionization (EoR). Together with ancillary data from earlier works, we derive
the dust temperature, $T_{\rm d}$, and mass, $M_{\rm d}$, of A1689-zD1 using
both traditional modified blackbody spectral energy density fitting, and a new
method that relies only on the [C ii] $158\ \mathrm{\mu m}$ line and
underlying continuum data. The two methods give $T_{\rm
d}=(42^{+13}_{-7},40^{+13}_{-7}$) K, and $M_{\rm
d}=(1.7^{+1.3}_{-0.7},2.0^{+1.8}_{-1.0})\,\times{}\,10^{7}\,M_{\odot}$. Band 9
observations improve the accuracy of the dust temperature (mass) estimate by
$\sim 50$% (6 times). The derived temperatures confirm the reported increasing
$T_{\rm d}$-redshift trend between $z=0$ and $8$; the dust mass is consistent
with a supernova origin. Although A1689-zD1 is a normal UV-selected galaxy,
our results, implying that $\sim$85% of its star formation rate is obscured,
underline the non-negligible effects of dust in EoR galaxies.

###### keywords:

galaxies: high-redshift — galaxies: evolution — galaxies: formation —
galaxies: individual (A1689-zD1) — ISM: dust, extinction — submillimetre:
galaxies

††pubyear: 2021††pagerange: Accurate dust temperature determination in a
$z=7.13$ galaxy–Acknowledgements

##  1 Introduction

Atacama Large Millimeter/submillimeter Array (ALMA) observations have revealed
the presence of dust in galaxies approaching the epoch of reionization (EoR;
e.g., Capak et al. 2015; Willott et al. 2015; Barisic et al. 2017; Laporte et
al. 2017). This was somewhat surprising, since UV studies mapping out the
cosmic star-formation rate density (SFRD) to $z\sim 10$ suggested a dearth of
dust at the high-redshift end based on the blue UV slopes of low-stellar mass
high-$z$ galaxies ($\beta_{\rm UV}$; e.g., Finkelstein et al. 2015; Bouwens et
al. 2016). Initially, the strong FIR emission at $z>7$ revealed by ALMA
observations was attributed to the presence of unexpectedly large dust masses
($M_{\rm d}$) in the observed high-$z$ galaxies, which was hard to reconcile
with known dust production mechanisms that operate on that timescale
(predominantly SN and grain growth; see Leśniewska & Michałowski 2019 and
references therein for the latest constraints).

This resulted in the so-called dust budget crisis, which also impacted star-
formation history (SFH) estimates of high-redshift galaxies (e.g., Mawatari et
al., 2020; Roberts-Borsani et al., 2020). The stringent constraints on SNe
dust production, coupled with the large deduced dust masses at $z>7$, required
very early stellar populations originating at $z\sim 14$ (Tamura et al.,
2019). However, the conclusions on the dust masses were heavily dependent on
the assumed (cold) dust temperatures ($T_{\rm d}\sim 30-40\ \mathrm{K}$) for
these high-$z$ sources, since in most cases only a single data point was
available in the FIR continuum. Recent observations (e.g., Bakx et al. 2020)
and theoretical studies (e.g. Behrens et al. 2018; Sommovigo et al. 2020) have
suggested the presence of warm dust in several high-$z$ galaxies ($T_{\rm
d}>60\ \mathrm{K}$), alleviating the large dust mass requirements set by their
observed $L_{\rm FIR}$ ($M_{\rm d}\propto T_{\rm d}^{\rm-(4+\beta_{d})}$ at
fixed $L_{\rm FIR}$, where typically $1.0<\beta_{\rm d}<3.0$). Unfortunately,
the large uncertainties on derived $T_{\rm d}$ at high-$z$ still hinder
accurate SFH studies.

Partially due to the lack of knowledge on the dust temperature at high-$z$,
the total fraction of obscured star-formation beyond $z>4$ is also largely
unknown (Novak et al., 2017; Casey et al., 2018; Bouwens et al., 2020;
Gruppioni et al., 2020; Schouws et al., 2021; Talia et al., 2021; Zavala et
al., 2021). This has strong implications for the cosmic SFRD; for example,
some of these recent works suggest that  there is no steep drop-off in SFRD at
$z>3$ (e.g., Gruppioni et al. 2020), which could indicate that we might be
underestimating the contribution of highly obscured systems to the SFRD at
$z>3$ due to the bias towards UV bright objects. On top of that, most studies
calculate the obscured star-formation rates and far-infrared luminosities of
single sources either by assuming a dust temperature, and/or by scaling
directly from the infrared excess (${\rm IRX}=L_{\rm FIR}/L_{\rm
UV}$)-$\beta_{\rm UV}$ relation. Both approaches suffer from the inherent
uncertainty in dust temperature (since obscured SFR and IRX both scale with
$T^{4+\beta_{\rm d}}$).

Moreover, the validity of IRX-$\beta_{\rm UV}$ relation at high-$z$ demands
that the UV and dust-emitting regions to be co-spatial, relying on the
absorbed UV emission to be re-emitted at far-infrared wavelengths. However,
observations suggest the possibility of spatial separation between these
regions in several sources  at $z=4-6$ (e.g., Faisst et al. 2017) and at
$z\sim 7-8$ sources (e.g., Carniani et al. 2017a; Laporte et al. 2019 and
Tamura et al. in prep.). In fact, this spatial separation scenario between UV
and IR is also supported by theoretical studies and simulations (Behrens et
al., 2018; Cochrane et al., 2019; Liang et al., 2019; Sommovigo et al., 2020).
A deviating IRX-$\beta_{\rm UV}$ relation would impact the results of galaxies
at high-$z$ (Le Fèvre et al., 2020; Fudamoto et al., 2020), and will impact
re-emission studies (e.g., MAGPHYS; da Cunha et al. 2008 and CIGALE; Boquien
et al. 2019) which will be prevalent in the ALMA + JWST era.

In this letter, we use the band 9 observations to estimate the dust properties
of a $z=7.1$ galaxy from the spectrum directly in order to measure the
obscured star-formation directly. We describe the source and data in Section
2, the fitting techniques in Section 3 and the implications in Section
4.111Throughout this paper, we assume a flat $\Lambda$-CDM cosmology with the
best-fit parameters derived from the Planck results (Planck Collaboration et
al., 2016), which are $\Omega_{\mathrm{m}}=0.307$,
$\Omega_{\mathrm{\Lambda}}=0.693$ and $h=0.678$.

Table 1: Continuum and fitting properties of A1689-zD1

|  $\lambda$ [mm] | F${}_{\nu}^{\rm int}$ [$\mu$Jy]† | Reference
---|---|---|---
Band 9 | 0.427 | 154 $\pm$ 37 | This work
Band 8 | 0.728 | 180 $\pm$ 39 | Inoue et al. (2020)
Band 7 | 0.873 | 143 $\pm$ 15 | Knudsen et al. (2017)
Band 6 | 1.33 | 60 $\pm$ 11 | Watson et al. (2015)
| SED fit | [C ii]-based | SED fit (no B9)
$T_{\rm d}$ [K] | 42${}^{+13}_{-7}$ | $40_{-7}^{+13}$ | $38_{-13}^{+35}$
$\beta_{\rm d}$ | 1.61${}^{+0.60}_{-0.75}$ |  2.03 ‡ | 1.78${}^{+0.55}_{-0.97}$
$M_{\rm d}$ [$10^{7}\ \mathrm{M_{\odot}}$]† | 1.7${}^{+1.3}_{-0.7}$ | $2.0^{+1.8}_{-1.0}$ | $2.4^{+11}_{-1.9}$
$L_{\rm FIR}$ [$10^{11}\ \mathrm{L_{\odot}}$]† | 1.9${}^{+0.5}_{-0.4}$ | 2.2${}^{+4.1}_{-1.0}$ | 1.5${}^{+3.0}_{-0.6}$
${\rm log_{10}\ IRX}$ | 1.0${}^{+0.1}_{-0.1}$ | 1.0${}^{+0.5}_{-0.3}$ | 0.8${}^{+0.5}_{-0.2}$
$M_{\rm d}$/SN [$\mathrm{M_{\odot}}$]† | 0.4${}^{+0.3}_{-0.1}$ | $0.6^{+0.6}_{-0.3}$ | 0.8${}^{+3.0}_{-0.7}$

Notes: $\,{}^{\dagger}$ Corrected for the magnification assuming $\mu=9.3$
from Knudsen et al. (2017). ‡ $\beta_{\rm d}$ is fixed to 2.03.

##  2 Target and Observations

A1689-zD1 was identified in Bradley et al. (2008) as a bright (m${}_{\rm
AB}\sim 25$) $z>7$ galaxy. Due to the foreground galaxy cluster (A1689;
Struble & Rood 1999), it is magnified222While $\mu$ is high, there is only
little shear, and we do not account for any differential lensing effects in
this paper by $\mu\simeq 9.3$ (Knudsen et al., 2017). Its intrinsic UV
magnitude indicates it is a sub-$L$* galaxy representing the bulk of galaxies
at $z=7$ (Ono et al., 2018). Band 6 observations at 1.3 mm by Watson et al.
(2015) reported the first detection of dust beyond redshift 7, and indicated
an intrinsic star-formation rate of $\sim 12$ $\rm M_{\odot}$ $\rm yr^{-1}$.
Notably, the estimated dust mass of this normal galaxy (assuming 35 K) was
found in tension to star-formation history and dust production estimates in
Leśniewska & Michałowski (2019).

In this letter, we combine the existing data on A1689-zD1 reported in Watson
et al. (2015), Knudsen et al. (2017) and Inoue et al. (2020) with archival
band 9 data from (Program ID: 2019.1.01778.S, P.I. D. Watson), see Table 1. We
use the [C ii] luminosity as reported in Knudsen et al. in prep., which is
(6.1 $\pm$ 0.7)$\ \times\ 10^{8}$ L⊙, and use their value for spectroscopic
redshift, $z_{\rm spec}=7.13$.

For the band 9 (Baryshev et al., 2015) data, the source was observed for 95
min. in baselines ranging from 14 to 312 m. The lower and upper sidebands
covered the contiguous frequency ranges of 690.4 to 697.6 and of 706.5 to
713.6 GHz. We assume a typical flux accuracy of 10%. The continuum image is
produced with CASA pipeline version 5.6.1-8 (McMullin et al., 2007), using
natural weighting, a taper of 0.5 arcseconds, and excluding any channels
within 1000 km/s of the [O iii] 52$\mu$m emission at 711.4 GHz. Figure 1 shows
the resulting image with a 0.61 by 0.67 arcsecond beam with a beam position
angle of 75 degrees, with an r.m.s. level of 210 $\mu$Jy beam-1. Using CASA’s
IMFIT routine, we spatially integrate the emission using a 2D Gaussian
profile. This results in a flux of 1.43 $\pm$ 0.31 mJy ($\sim 4.6\sigma$;
excluding calibration flux), with an apparent (or lensed) beam-deconvolved
size of 0.81 $\pm$ 0.26 by 0.38 $\pm$ 0.22 arcsec at a position angle of 44
$\pm$ 38 degrees.  The emission appears co-spatial to the UV-emission seen in
Knudsen et al. (2017), although we leave further discussion of this to Knudsen
et al. in prep.

Figure 1: The tapered band 9 data (background and black contours; drawn at -3,
-2, 2 and 3 $\sigma$) is shown against the band 8 continuum emission (white
contours; drawn at 5, 7 and 10 $\sigma$). The continuum emissions appears co-
spatial, and we find a 4.6$\sigma$ dust detection in band 9.

##  3 Methods

###  3.1 Spectral fitting

Figure 2 shows the modified black-body  (eq. 8 in Sommovigo et al. 2021)
fitted to the continuum points reported in Table 1. We use equations 12 and 18
from da Cunha et al. (2013) to account for the heating of dust by and
decreasing contrast against the CMB, respectively. We approximate the dust
mass absorption coefficient ($\kappa_{\rm\nu}$) as
$\kappa_{\rm\star}\left(\nu/\nu_{\star}\right)^{\beta_{\rm d}}$, with
($\kappa_{\rm\star},\nu_{\star})$ as (10.41 cm2/g, 1900 GHz) from Draine
(2003). We use the emcee MCMC-fitting routine, and allow $M_{\rm d}$, $T_{\rm
d}$ and $\beta_{\rm d}$ to vary freely using flat priors, resulting in a
magnification-corrected dust mass of 1.7${}^{+1.3}_{-0.7}$ $\times{}10^{7}$
M⊙, a dust temperature of 42${}^{+13}_{-7}$ K and a $\beta_{\rm d}$ of
1.61${}^{+0.60}_{-0.75}$. We note that the spectrum appears well-represented
by a single modified black-body. For comparison, we also fit the spectral
energy distribution (SED) without band 9 data, with an upper limit on
$\beta_{\rm d}$ of 2.5 to ensure convergence, and find significantly larger
errors across the board.  If we take a fiducial $\beta_{\rm d}=1.8$ (e.g.,
Casey 2012; Faisst et al. 2020), we find a more accurate dust temperature of
39${}^{+4}_{-5}$ K, however there is no improvement on the error of dust mass
(2.0${}^{+1.4}_{-0.7}\times 10^{7}$ M⊙) nor luminosity
(1.7${}^{+0.5}_{-0.4}\times 10^{11}$ L⊙). The accuracy of these later
parameters thus depends solely on observational uncertainties, indicating that
we fully trace the dust emission in this source.

Figure 2: We fit a modified black-body (red line and fill) to the observed
data points of A1689-zD1, including the band 9 data point (star). The [C
ii]-based spectrum (blue line and fill) is fit solely to the 158 $\mu$m
continuum data point (blue), and it predicts a consistent galaxy spectrum,
providing confidence in the [C ii]-based method for this specific source even
at 50 $\mu$m rest-frame. For comparison, the dashed red lines show the spread
in SEDs fitted without band 9 data, which results in a twice larger error in
dust temperature,  and a sixfold increase in the error in dust mass.

###  3.2 Dust temperature from [C ii] emission

We use the novel method proposed in Sommovigo et al. (2021) to derive the dust
temperature in galaxies, based on the combination of $1900\ \mathrm{GHz}$
continuum and the overlying [C ii] line emission. We provide a brief summary
of this method below; for further details and verification of this method on
19 local galaxies, three galaxies at $z\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}4$, and a $z\sim 6.7$ simulated
galaxy, we refer to Sommovigo et al. (2021).

We relate the observed [C ii] luminosity to the total dust mass via the gas
mass and a gas-to-dust ratio (assumed to scale linearly with the metallicity,
which is justified down to $Z\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}0.1\ Z_{\odot}$, see James et al.
2002; Draine & Li 2007; Galliano et al. 2008; Leroy et al. 2011). The gas mass
and [C ii] luminosity are related through a conversion factor $M_{\rm
gas}=\alpha_{\rm CII}\ L_{\rm CII}$. This conversion factor $\alpha_{\rm CII}$
is analytically derived from the combination of the de Looze relation (De
Looze et al., 2014) and the Kennicutt–Schmidt relation (Kennicutt, 1998,
hereafter, KS). Two parameters are added in the expression for $\alpha_{\rm
CII}$ in order to account for both (i) the expected offset from the KS-
relation (i.e., the burstiness of the SF of a galaxy parametrized by
$\kappa_{\rm s}$) and (ii) the observed larger extension of [C ii] with
respect to stellar emission at high-$z$ (up to 1.5-3 times larger; Carniani et
al., 2017b, 2018, 2020; Matthee et al., 2017, 2019; Fujimoto et al., 2019,
2020, 2021; Ginolfi et al., 2020; Herrera-Camus et al., 2021).

We fit a modified blackbody to derive the dust temperature using the
neighbouring continuum emission at $\sim$1900 GHz rest-frame wavelength,
assuming a fixed $\beta_{\rm d}=2.03$, which is based on the Draine (2003)
predictions for the Milky Way and the Small Magellanic Cloud. Within this
fitting routine, both the burstiness parameter ($\kappa_{\rm s}$) and the
metallicity are largely uncertain. In order to constrain the dust temperature,
two broad physical constraints are placed: (i) The dust mass cannot exceed the
maximal dust mass producible by supernovae (SNe), assuming all the SNe metal
yield ($\sim 2\ \mathrm{M_{\odot}}$ per SN) ends up locked into dust grains;
(ii) The dust-obscured star formation333This IR luminosity-to-SFR conversion
factor, $1.73\times 10^{-10}\ \mathrm{M_{\odot}yr^{-1}/L_{\odot}}$, is valid
for a Salpeter $1-100\ \mathrm{M_{\odot}}$ IMF, which we assume consistently
throughout the paper (Kennicutt, 1998; Madau & Dickinson, 2014), cannot
significantly (by 1 order of magnitude) exceed the SFR deduced from [C ii]
using the relation from De Looze et al. (2014) for starbursts. Applying our
method to A1689-zD1, we find a dust temperature $T_{\rm d}=40_{-7}^{+13}\
\mathrm{K}$ and mass of $M_{\rm d}=2.0^{+1.8}_{-1.0}\ \times 10^{7}\
\mathrm{M_{\odot}}$. These values are obtained assuming a wide range of values
for the metallicity $Z=0.2-1\ \mathrm{Z_{\odot}}$, and the burstiness
parameter $\kappa_{\rm s}=1-50$ (Vallini et al., 2021).  For our further
discussion of dust production mechanisms, we note that the removal of the dust
production constraint does not influence the derived quantities.

##  4 Implications

Figure 3: Observed peak wavelength $\lambda_{\rm peak,obs}$ is shown against
dust temperature $T_{\rm d}$ for a given dust emissivity index, $\beta_{\rm
d}$. The grey shaded region shows $\lambda_{\rm peak,obs}(T_{\rm d})$ at
redshift $z=7.1$ for $\beta_{\rm d}=2.0-1.0$. We show the sources with
reported dust temperatures beyond $z>5$ (Faisst et al. 2020; Harikane et al.
2020; Sugahara et al. 2021; Bakx et al. 2020; Laporte et al. 2019). The shaded
regions show the wavelength ranges probed by ALMA bands $7-10$. Without band
9,  we cannot probe the FIR peak on both sides and accurately estimate $T_{\rm
d}$ through SED fitting, while for lower-redshift observations band 10 might
even be required to accurately trace the SED.

The dust temperature and mass estimates from the [C ii]-based method agree
with the results from the direct SED fitting, which adds confidence to the
method from Sommovigo et al. (2021). As shown in Table 1, band 9 observations
reduce the uncertainty in the dust temperature by $\sim 50\%$, which
translates to much-improved estimate on the far-infrared luminosity and dust
mass estimate. In Figure 3, we show the observed peak emission wavelength
($\lambda_{\rm peak,obs}$) of galaxies at $z>5$ against our current best-
estimates for $T_{\rm d}$ and $\beta_{\rm d}$. To guide the eye, we include
the trend of $\lambda_{\rm peak,obs}$ with $T_{\rm d}$ for $\beta_{\rm d}=1-2$
at $z=7.1$. We also overlay the wavelengths of the ALMA bands 7 through 10. We
calculate this $\lambda_{\rm peak,obs}$ using

$\begin{split}\frac{\lambda_{\rm peak,obs}}{\mathrm{mm}}&=\frac{14.42\ (1+z)\
(T_{\rm d}/K)^{-1}}{W(-a\ e^{-a})+a}\ ,\end{split}$

where $a=3+\beta_{\rm d}$ and $W$ is the Lambert $W$ function. This is the
wavelength at which the continuum spectrum $F_{\rm\nu}$ peaks in frequency
units (e.g., Figure 2). This is an important distinction to keep in mind when
visualising $\lambda_{\rm obs,peak}$ from the analogy to Wien’s law, which
provides the peak of the spectrum when reported in wavelength units
$F_{\lambda}$.

Particularly for galaxies at lower redshifts and at higher temperatures,
short-wavelength observations are crucial to estimate the dust temperature,
whereas band 8 might be able to probe the emission peak for cold $z>8$
galaxies. In the foreseeable future, the high bands of ALMA (9 and 10) are the
only instrument capable of probing this regime, until such missions as the
Origins Space Telescope444https://asd.gsfc.nasa.gov/firs/docs/ (Meixner et
al., 2019).

Figure 4: Dust temperature $T_{\rm d}$ in “normal” (main-sequence) galaxies as
a function of redshift. The newest $T_{\rm d}$ estimates for A1689-zD1 are
shown in red (star for SED-fit and hexagon for the [C ii]-based result). Dust
temperatures obtained from stacked SEDs (blue and green circles) increase
linearly with redshift up to $z=6$. We highlight all the continuum detected
sources at $z>7$ with small post-stamps, and include their estimated $T_{\rm
d}$ based on both SED fits (triangles) and the [C ii]-based method (hexagons;
Faisst et al. 2020; Harikane et al. 2020; Sugahara et al. 2021; Bakx et al.
2020; Laporte et al. 2019). The addition of band 9 data significantly reduces
the uncertainty on the dust temperature of this source with respect to the
other high-$z$ sources, which are not observed in that band.

Accurate estimates of the the dust-obscured fraction of the star-formation
rate require strong constraints on the dust temperature, as $\mathrm{SFR}_{\rm
obsc}\propto L_{\rm FIR}\propto M_{\rm d}T_{\rm d}^{4+\beta_{\rm d}}$. Our
band 9 observations confirm that this relatively-cold ($T_{\rm d}$ $\sim
40-60$ K) system has a very large obscured fraction of the SFR around
$\sim$85% ($\mathrm{SFR}_{\rm obs}=33\pm 9\ \mathrm{M_{\odot}yr^{-1}}$,
whereas555derived using the magnification-corrected $L_{\rm UV}/10^{10}\
\mathrm{L_{\odot}}=2.28\pm 0.1$ (Hashimoto et al., 2019), and the UV
luminosity-to-unobscured SFR conversion factor in Madau & Dickinson (2014).
$\mathrm{SFR}_{\rm UV}=5.7\pm 0.3\ \mathrm{M_{\odot}yr^{-1}}$), even though it
was selected to be UV-bright.  The dust-obscured ratio is higher than the 61%
found for the typically more-massive ALPINE survey (Le Fèvre et al., 2020;
Faisst et al., 2020; Fudamoto et al., 2020; Béthermin et al., 2020; Khusanova
et al., 2021) at $z=5.5$, although this 61% is expected to decrease with
higher redshift. Albeit extreme, our dust-obscured ratio is in line with
recent results (both theoretical and observational) suggesting that we might
have been underestimating the dust-obscured contribution to the total SFR in
$z>4$ galaxies (see e.g. Novak et al., 2017; Gruppioni et al., 2020). On the
other hand, some studies of similarly-massive, UV-bright sources at very
high-$z$ ($z\sim 7$; e.g., Bouwens et al. 2021; Schouws et al. 2021) have so
far failed to detect dust continuum at 158 $\mu$m in half of their sources,
even though their average stellar mass is similar to those of A1689-zD1. These
undetected sources might have low dust contents, but that does not guarantee
low obscured star-formation fractions, since it is possible that this dust is
warm and is mainly emitting at wavelengths shorter than the currently observed
ones (mainly $158\ \mathrm{\mu m}$ rest-frame). In fact, while the continuum
around $158\ \mathrm{\mu m}$ of MACS0416_Y1 (Tamura et al., 2019; Bakx et al.,
2020) has yet to be seen, its spectrum is indicative of a similar obscured
fraction to A1689-zD1 (i.e., 94 to 85% for $\beta_{\rm d}=1.5$ to 2), even
though Y1’s UV-observed stellar mass is one order of magnitude lower than
A1689-zD1.

The selection towards UV-bright sources might also bias towards lower
fractions of obscured-to-total star-formation rate. With the discovery of so-
called optically-dark galaxies (e.g., Simpson et al. 2014; Franco et al. 2018;
Hatsukade et al. 2018; Wang et al. 2019; Yamaguchi et al. 2019; Williams et
al. 2019; Algera et al. 2020; Romano et al. 2020; Talia et al. 2021; Toba et
al. 2020; Umehata et al. 2020; Zhou et al. 2020; Shibuya et al. 2021; Smail et
al. 2021), we know of the existence of galaxies without detections in optical
wavelengths at high redshift. These sources, by definition, have exceedingly
high obscured fractions and might well account for a substantial fraction of
the SFRD at high redshift (Alcalde Pampliega et al., 2019; Gruppioni et al.,
2020; Zavala et al., 2021). The typical obscured star-formation rate fraction
across all $z>7$ galaxies might thus be higher than predicted by UV-selected
samples alone, with for example Gruppioni et al. (2020) predicting an increase
in SFRD of 17% at $z=5$ by this population.

Recently, attempts to quantify dust obscuration at high-$z$ have used a linear
scaling between dust temperature (and $L_{\rm IR}$ given a fixed $\beta_{\rm
d}$) and redshift (see e.g. Schreiber et al., 2018; Bouwens et al., 2020;
Vijayan et al., 2021). Other recent works have suggested that this linearly
increasing $T_{\rm d}-z$ trend flattens at $z>4$ (Liang et al., 2019; Faisst
et al., 2020). In Figure 4 we show the reported linear evolution of the dust
temperature with redshift, adding our latest results for A1689-zD1, and where
available, include the results from the method in Sommovigo et al. (2021). The
observed dust temperature for A1689-zD1 is compatible with both a flattening
(Liang et al., 2019; Faisst et al., 2020) and a linear (Schreiber et al.,
2018) $T_{\rm d}-z$ evolution. Meanwhile, the exceedingly-large scatter in
$T_{\rm d}$ at the highest redshifts (particularly at $z>7$) prevents us from
reaching a definitive conclusion on this observed evolution. Much of this
scatter is due to observational limitations, and only through further short-
wavelength observations of galaxies beyond $z>7$ can we distinguish the
possible scenarios. Part of the scatter could also be due to a larger source-
to-source variation in $T_{\rm d}$,  which is for example seen by the large
diversity of galaxies among the typically more-massive ALPINE sources (Le
Fèvre et al., 2020). Such source-to-source variation can only be identified by
larger unbiased samples looking at the dust-obscured star-formation at high
redshift. Here we note that an increased intrinsic scatter in dust temperature
would significantly boost the resulting dust-obscured star-formation rate,
given their strong dependence of star-formation rate on dust temperature,
similar to an Eddington-type bias.

Due to the large obscured faction of the SFR in A1689-zD1, one might naively
expect that this galaxy also contains an exceedingly large dust mass. Instead,
the dust mass derived from SED fitting implies a dust yield of $y_{\rm
d}=0.4^{+0.3}_{-0.1}\ \mathrm{M_{\odot}}$ per SN. This estimate is almost an
order of magnitude more accurate than the one derived without band 9 data, and
most importantly, it is consistent with latest SN dust production constraints
by Leśniewska & Michałowski (2019) based on the expected number of SNe given
its stellar mass estimate. They find at most a $y_{\rm d}=1.1\
\mathrm{M_{\odot}}$ per SN, derived in the extreme case of no dust
destruction/ejection. We note that SN yield is still highly debated, with
other works suggesting that dust destruction processes might only spare $0.1\
\mathrm{M_{\odot}}$ per SN (e.g., Matsuura et al., 2015, 2019; Slavin et al.,
2020).  In this extreme case, inter-stellar medium grain growth (Mancini et
al., 2015; Michałowski, 2015) or more exotic dust production mechanisms might
well be required at $z>7$, such as dust produced in supershells (e.g.,
Martínez-González et al., 2021) or in the wake of Wolf-Rayet stars (e.g., Lau
et al., 2021).

## References

  * Alcalde Pampliega et al. (2019) Alcalde Pampliega B., et al., 2019, ApJ, 876, 135
  * Algera et al. (2020) Algera H. S. B., et al., 2020, ApJ, 903, 139
  * Bakx et al. (2020) Bakx T. J. L. C., et al., 2020, MNRAS, 493, 4294
  * Barisic et al. (2017) Barisic I., et al., 2017, ApJ, 845, 41
  * Baryshev et al. (2015) Baryshev A. M., et al., 2015, A&A, 577, A129
  * Behrens et al. (2018) Behrens C., Pallottini A., Ferrara A., Gallerani S., Vallini L., 2018, MNRAS, 477, 552
  * Béthermin et al. (2020) Béthermin M., et al., 2020, A&A, 643, A2
  * Boquien et al. (2019) Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue A. K., Salas H., 2019, A&A, 622, A103
  * Bouwens et al. (2016) Bouwens R. J., et al., 2016, ApJ, 833, 72
  * Bouwens et al. (2020) Bouwens R., et al., 2020, ApJ, 902, 112
  * Bouwens et al. (2021) Bouwens R. J., et al., 2021, arXiv e-prints, p. arXiv:2106.13719
  * Bradley et al. (2008) Bradley L. D., et al., 2008, ApJ, 678, 647
  * Capak et al. (2015) Capak P. L., et al., 2015, Nature, 522, 455
  * Carniani et al. (2017a) Carniani S., et al., 2017a, A&A, 605, A42
  * Carniani et al. (2017b) Carniani S., et al., 2017b, A&A, 605, A42
  * Carniani et al. (2018) Carniani S., et al., 2018, MNRAS, 478, 1170
  * Carniani et al. (2020) Carniani S., et al., 2020, MNRAS, 499, 5136
  * Casey (2012) Casey C. M., 2012, MNRAS, 425, 3094
  * Casey et al. (2018) Casey C. M., et al., 2018, ApJ, 862, 77
  * Cochrane et al. (2019) Cochrane R. K., et al., 2019, MNRAS, 488, 1779
  * De Looze et al. (2014) De Looze I., et al., 2014, A&A, 568, A62
  * Draine (2003) Draine B., 2003, Annual Review of Astronomy and Astrophysics, 41, 241
  * Draine & Li (2007) Draine B. T., Li A., 2007, ApJ, 657, 810
  * Faisst et al. (2017) Faisst A. L., et al., 2017, ApJ, 847, 21
  * Faisst et al. (2020) Faisst A. L., Fudamoto Y., Oesch P. A., Scoville N., Riechers D. A., Pavesi R., Capak P., 2020, MNRAS, 498, 4192
  * Finkelstein et al. (2015) Finkelstein S. L., et al., 2015, ApJ, 810, 71
  * Franco et al. (2018) Franco M., et al., 2018, A&A, 620, A152
  * Fudamoto et al. (2020) Fudamoto Y., et al., 2020, A&A, 643, A4
  * Fujimoto et al. (2019) Fujimoto S., et al., 2019, ApJ, 887, 107
  * Fujimoto et al. (2020) Fujimoto S., et al., 2020, ApJ, 900, 1
  * Fujimoto et al. (2021) Fujimoto S., et al., 2021, ApJ, 911, 99
  * Galliano et al. (2008) Galliano F., Dwek E., Chanial P., 2008, ApJ, 672, 214
  * Ginolfi et al. (2020) Ginolfi M., et al., 2020, A&A, 633, A90
  * Gruppioni et al. (2020) Gruppioni C., et al., 2020, A&A, 643, A8
  * Harikane et al. (2020) Harikane Y., et al., 2020, ApJ, 896, 93
  * Hashimoto et al. (2019) Hashimoto T., et al., 2019, Pub. Astron. Soc. Japan, p. 70
  * Hatsukade et al. (2018) Hatsukade B., et al., 2018, Pub. Astron. Soc. Japan, 70, 105
  * Herrera-Camus et al. (2021) Herrera-Camus R., et al., 2021, A&A, 649, A31
  * Inoue et al. (2020) Inoue A. K., Hashimoto T., Chihara H., Koike C., 2020, MNRAS, 495, 1577
  * James et al. (2002) James A., Dunne L., Eales S., Edmunds M. G., 2002, MNRAS, 335, 753
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, ApJ, 498, 541
  * Khusanova et al. (2021) Khusanova Y., et al., 2021, A&A, 649, A152
  * Knudsen et al. (2017) Knudsen K. K., Watson D., Frayer D., Christensen L., Gallazzi A., Michałowski M. J., Richard J., Zavala J., 2017, MNRAS, 466, 138
  * Laporte et al. (2017) Laporte N., et al., 2017, ApJL, 837, L21
  * Laporte et al. (2019) Laporte N., et al., 2019, MNRAS, 487, L81
  * Lau et al. (2021) Lau R. M., et al., 2021, ApJ, 909, 113
  * Le Fèvre et al. (2020) Le Fèvre O., et al., 2020, A&A, 643, A1
  * Leroy et al. (2011) Leroy A. K., et al., 2011, The Astrophysical Journal, 737, 12
  * Leśniewska & Michałowski (2019) Leśniewska A., Michałowski M. J., 2019, A&A, 624, L13
  * Liang et al. (2019) Liang L., et al., 2019, MNRAS, 489, 1397
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Mancini et al. (2015) Mancini M., Schneider R., Graziani L., Valiante R., Dayal P., Maio U., Ciardi B., Hunt L. K., 2015, MNRAS, 451, L70
  * Martínez-González et al. (2021) Martínez-González S., Silich S., Tenorio-Tagle G., 2021, MNRAS, 507, 1175
  * Matsuura et al. (2015) Matsuura M., et al., 2015, ApJ, 800, 50
  * Matsuura et al. (2019) Matsuura M., et al., 2019, MNRAS, 482, 1715
  * Matthee et al. (2017) Matthee J., et al., 2017, The Astrophysical Journal, 851, 145
  * Matthee et al. (2019) Matthee J., et al., 2019, The Astrophysical Journal, 881, 124
  * Mawatari et al. (2020) Mawatari K., et al., 2020, ApJ, 889, 137
  * McMullin et al. (2007) McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Astronomical Data Analysis Software and Systems XVI. p. 127
  * Meixner et al. (2019) Meixner M., et al., 2019, arXiv e-prints, p. arXiv:1912.06213
  * Michałowski (2015) Michałowski M. J., 2015, A&A, 577, A80
  * Novak et al. (2017) Novak M., et al., 2017, A&A, 602, A5
  * Ono et al. (2018) Ono Y., et al., 2018, Pub. Astron. Soc. Japan, 70, S10
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Roberts-Borsani et al. (2020) Roberts-Borsani G. W., Ellis R. S., Laporte N., 2020, MNRAS, 497, 3440
  * Romano et al. (2020) Romano M., et al., 2020, MNRAS, 496, 875
  * Schouws et al. (2021) Schouws S., et al., 2021, arXiv e-prints, p. arXiv:2105.12133
  * Schreiber et al. (2018) Schreiber C., Elbaz D., Pannella M., Ciesla L., Wang T., Franco M., 2018, A&A, 609, A30
  * Shibuya et al. (2021) Shibuya T., Miura N., Iwadate K., Fujimoto S., Harikane Y., Toba Y., Umayahara T., Ito Y., 2021, arXiv e-prints, p. arXiv:2106.03728
  * Simpson et al. (2014) Simpson J. M., et al., 2014, ApJ, 788, 125
  * Slavin et al. (2020) Slavin J. D., Dwek E., Mac Low M.-M., Hill A. S., 2020, ApJ, 902, 135
  * Smail et al. (2021) Smail I., et al., 2021, MNRAS, 502, 3426
  * Sommovigo et al. (2020) Sommovigo L., Ferrara A., Pallottini A., Carniani S., Gallerani S., Decataldo D., 2020, MNRAS, 497, 956
  * Sommovigo et al. (2021) Sommovigo L., Ferrara A., Carniani S., Zanella A., Pallottini A., Gallerani S., Vallini L., 2021, MNRAS, 503, 4878
  * Struble & Rood (1999) Struble M. F., Rood H. J., 1999, ApJS, 125, 35
  * Sugahara et al. (2021) Sugahara Y., et al., 2021, arXiv e-prints, p. arXiv:2104.02201
  * Talia et al. (2021) Talia M., Cimatti A., Giulietti M., Zamorani G., Bethermin M., Faisst A., Le Fèvre O., Smolçić V., 2021, ApJ, 909, 23
  * Tamura et al. (2019) Tamura Y., et al., 2019, ApJ, 874, 27
  * Toba et al. (2020) Toba Y., et al., 2020, ApJ, 899, 35
  * Umehata et al. (2020) Umehata H., et al., 2020, A&A, 640, L8
  * Vallini et al. (2021) Vallini L., Ferrara A., Pallottini A., Carniani S., Gallerani S., 2021, MNRAS, 505, 5543
  * Vijayan et al. (2021) Vijayan A. P., et al., 2021, arXiv e-prints, p. arXiv:2108.00830
  * Wang et al. (2019) Wang T., et al., 2019, Nature, 572, 211
  * Watson et al. (2015) Watson D., Christensen L., Knudsen K. K., Richard J., Gallazzi A., Michałowski M. J., 2015, Nature, 519, 327
  * Williams et al. (2019) Williams C. C., et al., 2019, ApJ, 884, 154
  * Willott et al. (2015) Willott C. J., Carilli C. L., Wagg J., Wang R., 2015, ApJ, 807, 180
  * Yamaguchi et al. (2019) Yamaguchi Y., et al., 2019, ApJ, 878, 73
  * Zavala et al. (2021) Zavala J. A., et al., 2021, ApJ, 909, 165
  * Zhou et al. (2020) Zhou L., et al., 2020, A&A, 642, A155
  * da Cunha et al. (2008) da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595
  * da Cunha et al. (2013) da Cunha E., et al., 2013, ApJ, 766, 13

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## Acknowledgements

We would like to thank the anonymous referee for their insightful comments and
suggested additions. This paper makes use of the following ALMA data:
ADS/JAO.ALMA 2011.1.00319.S, 2012.1.00216.S, 2013.1.01064.S, 2016.1.00954.S,
and 2019.1.01778.S. TB and YT acknowledge funding from NAOJ ALMA Scientific
Research Grant Numbers 2018-09B and JSPS KAKENHI No. 17H06130. LS, SC, AF, AP
acknowledge support from the ERC Advanced Grant INTERSTELLAR H2020/740120 (PI:
Ferrara). The Cosmic Dawn Center is funded by the Danish National Research
Foundation under grant No. 140. This project has received funding from the
European Union’s Horizon 2020 research and innovation program under the Marie
Sklodowska-Curie grant agreement No. 847523 ‘INTERACTIONS’. DW and SF are
supported by Independent Research Fund Denmark grant DFF–7014-00017 Any
dissemination of results must indicate that it reflects only the author’s view
and that the Commission is not responsible for any use that may be made of the
information it contains.


# Circular Rosenzweig-Porter random matrix ensemble

Wouter Buijsman<EMAIL_ADDRESS>Department of Physics, Ben-Gurion
University of the Negev, Beer-Sheva 84105, Israel  Yevgeny Bar Lev  Department
of Physics, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel

###### Abstract

The Rosenzweig-Porter random matrix ensemble serves as a qualitative
phenomenological model for the level statistics and fractality of eigenstates
across the many-body localization transition in static systems. We propose a
unitary (circular) analogue of this ensemble, which similarly captures the
phenomenology of many-body localization in periodically driven (Floquet)
systems. We define this ensemble as the outcome of a Dyson Brownian motion
process. We show numerical evidence that this ensemble shares some key
statistical properties with the Rosenzweig-Porter ensemble for both the
eigenvalues and the eigenstates.

##  I Introduction

Many-body localization is known as a robust mechanism of ergodicity breaking
in disordered interacting quantum many-body systems that can be thought of as
Anderson localization in Fock space [1, 2]. It is characterized by non-trivial
properties, such as the violation of the eigenstate thermalization hypothesis
[3] and the persistent memory of an initial state, among many others [4].
Significant progress towards an understanding of many-body localization from
both theoretical and experimental approaches has been witnessed during the
last one and a half decade [5, 6, 7].

Many-body localization can be mapped to single-particle localization on a Fock
space lattice [8], an approach which has been successfully utilized both
analytically [9, 10, 11, 12, 13, 14] and numerically [15, 16, 17, 18, 19].
Within the framework of this approach, obtaining analytical results is
notoriously difficult due to strong correlations between disorder potentials
[20, 21]. It was therefore realized that it is of interest to have simpler
phenomenological models for many-body localization that are analytically
tractable while still capturing the essential properties [22, 23, 24]. A
profound example of such a model is the Rosenzweig-Porter random matrix
ensemble [25]. This model was suggested as a qualitative phenomenological
model for the level statistics and fractality of eigenstates across the many-
body localization transition [23], which separates a many-body localized from
a thermal phase. This proposal triggered considerable interest for the
Rosenzweig-Porter ensemble over the last years [26, 27, 28, 29, 30, 31, 32,
33, 34, 35].

The notion of many-body localization extends to interacting periodically
driven (Floquet) systems [36, 37, 38, 39, 40], for which localization in Fock
space has been investigated analytically very recently [41]. While static
systems are described by Hamiltonians, which are Hermitian, periodically
driven systems are described by Floquet operators, which are unitary. This
raises the arguably natural question whether a unitary (or _circular_ since
the eigenvalues lie on the unit circle in the complex plane) analogue of the
Rosenzweig-Porter ensemble can be constructed.

In this work, we propose a circular analogue of the Rosenzweig-Porter
ensemble, which we define as the result of a Dyson Brownian motion process
[42]. We provide numerical evidence that this circular analogue has similar
features as the Rosenzweig-Porter ensemble by investigating some key
statistical properties of the eigenvalues and the eigenstates. Our work might
serve as a first step in constructing generalizations covering e.g.
multifractality, analog to recent physically motivated proposals for
generalizations of the Rosenzweig-Porter ensemble [29, 43, 44, 45, 46] and
other Floquet models with multifractal eigenstates [47, 48, 49, 50].

The structure of the paper is as follows. Sec. II reviews the Rosenzweig-
Porter ensemble and some of its properties. Sec. III discusses the
construction of the circular Rosenzweig-Porter ensemble. Sec. IV numerically
investigates a number of statistical properties of the circular analogue. Sec.
V closes the paper with conclusions and an outlook.

##  II Rosenzweig-Porter ensemble

The Rosenweig-Porter (RP) ensemble [23], which was originally proposed in the
context of complex atomic nuclei [25], can be seen as generalization of the
Gaussian orthogonal ensemble (GOE) [51, 52] with a preferential basis. The
ensemble consists of real symmetric matrices $H$ of the form

$H=H_{0}+\frac{\epsilon}{\sqrt{N^{\gamma}}}V,$ (1)

where $N$ is the matrix dimension, $H_{0}$ is diagonal with elements sampled
independently from the Gaussian distribution with mean $\mu=0$ and variance
$\sigma^{2}=1$, and $V$ is sampled from the GOE. The GOE consists of real-
valued symmetric matrices with the diagonal and upper triangular elements
sampled independently from the Gaussian distribution with mean $\mu=0$ and
variance $\sigma^{2}=2$ (diagonal elements) or $\sigma^{2}=1$ (upper
triangular elements). Interpreting $H$ as an Hamiltonian, the parameter
$\epsilon\sim\mathcal{O}(1)$ can be viewed of as a perturbation strength. The
positive-valued parameter $\gamma$ controls the relative strength of the
terms, and thus the properties of the ensemble. Below, we outline a number of
these properties, which we numerically explore for the circular analogue in
Sec. IV.

The statistics of spacings between consecutive energy levels provides a
standard diagnostic for quantum chaos [51, 3]. In the large-$N$ limit, level
spacings of the RP ensemble obey Wigner-Dyson level statistics for $\gamma<2$
[53], which are typically observed for chaotic quantum systems. For
$\gamma\geq 2$ and in the same limit, these statistics are Poissonian [53], as
typically observed for integrable (non-chaotic) systems. It is convenient to
probe level spacing statistics by the average ratio of consecutive level
spacings [54, 55]. Let $E_{n}$ denote the eigenvalues of $H$ sorted in
ascending order, such that the spacings $s_{n}$ of consecutive levels are
given by $s_{n}=E_{n+1}-E_{n}$. The ratios $r_{n}$ are defined as

$r_{n}=\min\bigg{(}\frac{s_{n+1}}{s_{n}},\frac{s_{n}}{s_{n+1}}\bigg{)}.$ (2)

On average, the ratios of consecutive level spacings acquire the values
$\overline{r}\approx 0.386$ for Poissonian and $\overline{r}\approx 0.530$ for
Wigner-Dyson (GOE) level statistics [55].

Next, we consider eigenstates of $H$ taken from the middle of the spectrum. We
focus on the basis in which $H_{0}$ is diagonal. We denote eigenstates of $H$
and $H_{0}$ by respectively $|\psi_{n}\rangle$ and $|n\rangle$. Localization
of an eigenstate $|\psi_{n}\rangle$ can be quantified by the inverse
participation ratios

$\text{IPR}_{q}=\sum_{m}\left|\langle m|\psi_{n}\rangle\right|^{2q}$ (3)

with $q>1/2$. Asymptotically in $N$, the inverse participation ratios scale as
$N^{-(q-1)D_{q}}$, where $D_{q}$ is known as the fractal dimension with
parameter $q$ [30]. Eigenstates of the RP ensemble are characterized by
$D_{q}=1$ for $\gamma\leq 1$, indicating that these are spread out over a
finite fraction of the Hilbert space as $\text{IPR}_{2}\sim N^{-1}$. For
$1<\gamma<2$, the eigenstates are fractal with $D_{q}=2-\gamma$. As $D_{q}$
does not depend on $q$, the eigenstates are fractal but not multifractal. For
$\gamma\geq 2$, the eigenstates are characterized by $D_{q}=0$, which reflects
that these eigenstates have significant overlap with only $\mathcal{O}(1)$
basis states [23].

Eigenstates of the RP ensemble obey Breit-Wigner statistics [30], which are
expected to apply rather generically to quantum many-body systems [56, 57, 58,
59]. In the large-$N$ limit, it follows that

$\overline{\left|\langle
m|\psi_{n}\rangle\right|^{2}}\sim\frac{1}{\big{(}E_{n}-E_{m}^{(0)}\big{)}^{2}+\Gamma^{2}(E_{n})},$
(4)

where the bar denotes an average over $V$ in Eq. (1), $E_{n}$ is the
eigenvalue corresponding to the eigenstate $|\psi_{n}\rangle$,
$E_{n}^{(0)}=(H_{0})_{nn}$, and the so-called spreading width $\Gamma(E_{n})$
is obtained through Fermi’s golden rule as

$\Gamma(E_{n})=\frac{\pi\epsilon^{2}}{N^{\gamma}}\rho(E_{n}),$ (5)

where $\rho(E_{n})=\sum_{i}\delta(E_{n}-E_{i})$ is the density of states in
the vicinity of the eigenvalue $E_{n}$. Eq. (4) describes what is known as the
_shape_ of the eigenstates [59]. For $1<\gamma<2$, the spreading width gives
the width of the so-called mini-band [60, 23, 31]. In this energy window, the
eigenstate amplitudes $|\langle\psi_{n}|m\rangle|^{2}$ are on average of order
$N^{-D_{2}}$ (where $D_{2}$ is the fractal dimension $D_{q}$ for $q=2$), which
is much larger than the value $N^{-1}$ obtained when averaging over all basis
states.

##  III Circular analogue

In this Section, we propose a unitary analogue of the RP ensemble, which we
refer to as the circular RP (CRP) ensemble. The starting point of our
construction is the observation made in Ref. [26] that the RP ensemble can be
seen as the outcome of a finite-time stochastic process for the matrix
elements, referred to as a Dyson Brownian motion process [42]. Let $M(t)$
denote a time-dependent, real symmetric matrix of dimension $N$, which evolves
stochastically according to

$dM(t)=\sqrt{dt}X,$ (6)

where $dM(t)=M(t+dt)-M(t)$ and $X$ is a GOE matrix sampled independently at
each infinitesimal time step $dt$. For the initial condition $M(0)=H_{0}$ with
$H_{0}$ as given in Eq. (1), one finds $M(t)=H_{0}+\sqrt{t}V$, where $V$ is a
sample from the GOE. From the solution $M(t)$ is it clear that the RP ensemble
results at $t=\epsilon^{2}/N^{\gamma}$. The GOE is invariant under
transformations of the basis [51, 52]. Using this, for an infinitesimaly small
time difference $dt$ one can obtain stochastic equations for the eigenvalues
$E_{n}(t)$ of $M(t)$ perturbatively as

$dE_{n}(t)=\sqrt{dt}X_{nn}+dt\sum_{m\neq
n}\frac{(X_{nm})^{2}}{E_{n}(t)-E_{m}(t)},$ (7)

and the stochastic equations for the corresponding eigenstates
$|\psi_{n}\rangle$ as

$\begin{split}d|\psi_{n}(t)\rangle&=\sqrt{dt}\sum_{m\neq
n}\frac{X_{mn}}{E_{n}(t)-E_{m}(t)}|\psi_{m}(t)\rangle\\\
&-\frac{dt}{2}\sum_{m\neq
n}\frac{(X_{mn})^{2}}{[E_{n}(t)-E_{m}(t)]^{2}}|\psi_{n}(t)\rangle,\end{split}$
(8)

where the increments are given by $dE_{n}(t)=E_{n}(t+dt)-E_{n}(t)$ and
$d|\psi_{n}(t)\rangle=|\psi_{n}(t+dt)\rangle-|\psi_{n}(t)\rangle$. In the
literature, the stochastic evolution of the eigenstates is known as the
_eigenvector moment flow_ [61, 62].

To construct a circular ensemble, we follow an approach introduced by Dyson
[42]. We introduce $S(t)$, which is a time-dependent symmetric unitary matrix
of dimension $N$ with eigenvalues $e^{i\theta_{n}(t)}$. Such a matrix can be
interpreted as the Floquet operator of a periodically driven system obeying
time-reversal symmetry. The matrix evolves stochastically according to

$dS(t)=i\sqrt{dt}U^{T}(t)XU(t),$ (9)

where $dS(t)=S(t+dt)-S(t)$, and $U(t)$ is a unitary matrix defined via the
decomposition $S(t)=U^{T}(t)U(t)$. The decomposition is defined up to
transformations $U\to OU$ for real orthogonal matrices $O$111The set
$U_{S}(N)$ of symmetric unitary matrices of dimension $N$ is isomorphic to
$U(N)/O(N)$, where $U(N)$ denotes the set of imaginary unitary matrices of
dimension $N$, and $O(N)$ denote the set of real orthogonal matrices of
dimension $N$. For a proof, see e.g. Proposition 2.2.4 of Ref. [52]., which
can be absorbed in $X$ as the GOE is invariant under basis transformations. It
is useful to take $U^{T}=U$, in which case $U(t)$ becomes $\sqrt{S(t)}$. Eq.
(9) conserves the unitarity of $S(t)$ only in the limit $dt\to 0$. Other,
numerically equivalent methods which conserve unitarity of $S(t)$ for any $dt$
can be also used (see Sec. IV). We propose the CRP ensemble as the outcome of
the stochastic process described by Eq. (9) at $t=\epsilon^{2}/N^{\gamma}$.
For the RP ensemble, $M(0)$ is diagonal with uncorrelated elements. This
motivates us to analogously take

$S(0)=\operatorname{diag}\left(e^{i\theta_{1}(0)},\dots,e^{i\theta_{N}(0)}\right)$
(10)

with the phases $\theta_{n}(0)$ sampled independently from the uniform
distribution ranging over $[-\pi,\pi)$. This non-unique choice respects the
uniform density of states as observed generically for Floquet operators. The
phases evolve according to a stochastic equation that can be obtained as

$d\theta_{n}(t)=\sqrt{dt}X_{nn}+dt\sum_{m\neq
n}\frac{\left(X_{nm}\right)^{2}}{2\tan\frac{1}{2}\left[\theta_{n}(t)-\theta_{m}(t)\right]},$
(11)

where $d\theta_{n}(t)=\theta_{n}(t+dt)-\theta_{n}(t)$.

Eqs. (7) and (11) are known to yield the same statistics in the middle of the
spectrum and in the large-$N$ limit, provided that the densities of states of
$M(t)$ and $S(t)$ are the same [64]. Heuristically, this can be understood by
noticing that the summations are dominated by contributions from the few most
nearby levels, which are of the order of the density of states,
$\rho(E_{n})\sim N$. While there are $\mathcal{O}(N)$ additional terms in
these summations, at the middle of the spectrum the contributions from levels
$E_{m}>E_{n}$ ($\theta_{m}>\theta_{n}$) and $E_{m}<E_{n}$
($\theta_{m}<\theta_{n}$) are respectively negative and positive, and
therefore provide a contribution with an overall sub-leading magnitude of
$\mathcal{O}(\sqrt{N})$. The density of states of $S(t)$ can be made equal to
those of $M(t)$ by scaling the latter by $1/\sqrt{2\pi}$. We therefore expect
the spectral statistics of the RP for a given $\epsilon$ at the middle of the
spectrum to match those of the CRP with $\epsilon\to\epsilon/\sqrt{2\pi}$. We
remark that $M(t)$ has a density of states that remains constant during the
Brownian motion process for $\gamma>1$ [30].

The eigenstates $|\psi_{n}\rangle$ of $S(t)$ associated with the eigenvalues
$e^{i\theta_{n}(t)}$ evolve stochastically according to

$\begin{split}d|\psi_{n}(t)\rangle&=\sqrt{dt}\sum_{m\neq
n}\frac{X_{mn}}{2\tan\tfrac{1}{2}[\theta_{n}(t)-\theta_{m}(t)]}|\psi_{m}(t)\rangle\\\
&-\frac{dt}{2}\sum_{m\neq
n}\frac{(X_{mn})^{2}}{4\tan^{2}\tfrac{1}{2}[\theta_{n}(t)-\theta_{m}(t)]}|\psi_{n}(t)\rangle,\end{split}$
(12)

where again $d|\psi_{n}(t)\rangle=|\psi_{n}(t+dt)\rangle-|\psi_{n}(t)\rangle$.
In line with the arguments for the eigenvalues, one might expect that the RP
and CRP have similar eigenstate statistics. We remark that no explicit
construction for the matrices resulting from circular Dyson Brownian motion
processes is known [64]. Moreover, we note that the CRP can not be constructed
as $S=e^{iH}$ for $H$ sampled from the RP ensemble. Indeed, for example then
the density of states of $S$ would not be uniform.

##  IV Numerical evaluation

The results below are obtained through a numerical evaluation of Eq. (9) up to
$t=\epsilon^{2}/N^{\gamma}$ using the algorithm described Ref. [52] (Sec.
11.2.1). This algorithm utilizes a different and numerically equivalent
expression for of the stochastic equations (9), given by

$S(t+dt)=S(t)e^{i\sqrt{dt}X}.$ (13)

This description has the computational advantages that no decomposition
$S=U^{T}U$ has to be made, and that $S(t)$ remains unitary even for finite
$dt$. In the numerical evaluations, we average over at least $10^{6}$
eigenvalues or eigenstates. We take $\epsilon=1$ and $10^{4}$ time steps to
reach $t=\epsilon^{2}/N^{\gamma}$. Empirically, this time step turns out to be
small enough for convergence of the results. We have validated that our
results are visibly indistinguishable from results obtained when taking the
number of time steps half as large, and found that no qualitative differences
can be observed when taking the number of time steps an order of magnitude
smaller. The numerical evaluation of Eq. (9) requires diagonalization or
matrix-matrix multiplications in each step, making it is computationally
expensive. We therefore restricted to matrix dimensions $N\leq 1000$, which
still provide sufficient numerical evidence for the claims of this work.

We denote the eigenphases for samples from the CRP ensemble by $\theta_{n}$,
which we label in ascending order ($\theta_{n+1}\geq\theta_{n}$). First, we
consider the average ratio of consecutive level spacings $\overline{r}$. The
ratios $r_{n}$ of consecutive level spacings are defined in Eq. (2), with
$E_{n}$ replaced by the phases $\theta_{n}$. The upper panel of Fig. 1 shows
$\overline{r}$ as a function of $\gamma$ for several matrix dimensions. As
discussed in Sec. II, the RP ensemble is characterized by $\overline{r}\approx
0.530$ for $\gamma<2$ (Wigner-Dyson) and $\overline{r}\approx 0.386$
(Poissonian) for $\gamma\geq 2$ in the large-$N$ limit, which is consistent
with our results. The level statistics exhibit a phase transition with
critical exponent $\nu=1$ from Wigner-Dyson to Poissonian at
$\gamma=\gamma_{c}=2$. For the RP ensemble, a scaling collapse of
$\overline{r}$ as a function of $(\gamma-\gamma_{c})\ln(N)^{1/\nu}$ has been
observed in Ref. [33]. The lower panel shows that this observation can be made
as well for the CRP ensemble.

Figure 1: The average ratio of consecutive level spacing spacings
$\overline{r}$ as a function of $\gamma$ (upper panel) and
$(\gamma-\gamma_{c})\ln(N)^{1/\nu}$ with $\gamma_{c}=2$ and $\nu=1$ (lower
panel) for matrix dimensions $N=250$, $500$, and $1000$. Poissonian and
Wigner-Dyson level statistics are characterized by respectively
$\overline{r}\approx 0.386$ and $\overline{r}\approx 0.530$ (dashed lines).

Next, we consider eigenstates of the CRP ensemble. We focus on the basis in
which $S(0)$ is diagonal. We denote eigenstates of $S$ and $S_{0}$ by
respectively $|\psi_{n}\rangle$ and $|n\rangle,$ with the eigenvalue
$e^{i\theta_{n}}$ corresponding to $|\psi_{n}\rangle$. We first focus on the
eigenstate inverse participation ratio $\text{IPR}_{2}$ as given in Eq. (3).
Asymptotically in $N$, the RP ensemble is characterized by $\text{IPR}_{2}\sim
N^{-1}$ for $\gamma\leq 1$, and $\text{IPR}_{2}=\mathcal{O}(1)$ for
$\gamma\geq 2$ as discussed in Sec. II. In the intermediate region
$1<\gamma<2$, the eigenstates are fractal, and are characterized by
$\text{IPR}_{2}\sim N^{-(2-\gamma)}$. Fig. 2 shows the average
$\overline{\text{IPR}}_{2}$ (upper panel) and
$N\times\overline{\text{IPR}}_{2}$ (lower panel) as a function of $\gamma$,
for the same matrix dimensions as studied above. The results we obtain for the
CRP ensemble are again consistent with the large-$N$ behavior for the RP
ensemble.

Figure 2: The average of the inverse participation ratio $\text{IPR}_{2}$
(upper panel) and the value scaled by $N$ (lower panel) as a function of
$\gamma$ for matrix dimensions $N=250$, $500$, and $1000$.

Finally, we consider the shape of the eigenstates. The eigenstates of the RP
ensemble obey Breit-Wigner statistics [30]. For the CRP ensemble we aim to
show [cf. Eq. (4)] that

$\overline{\langle
m|\psi_{n}\rangle|^{2}}\sim\frac{1}{\big{(}\theta_{n}-\theta_{m}^{(0)}\big{)}^{2}+\Gamma^{2}},$
(14)

where $\theta_{m}^{(0)}=\theta_{m}(0)$ as introduced in Eq. (10). Fig. 3 shows
the average $\overline{\langle m|\psi_{n}\rangle|^{2}}$ as a function of
$\theta_{n}-\theta_{m}^{(0)}$ for $\gamma=0.5$ (upper panel) and $\gamma=1.5$
(lower panel) for the largest matrix dimension $N=1000$. The smooth curves
show the evaluation of Eq. (14) with $\Gamma$ obtained from Eq. (5), and the
density of states is given by $\rho=1000/(2\pi)\approx 159.2$.

Figure 3: The average of $|\langle m|\psi_{n}\rangle|^{2}$ as a function of
$\theta_{n}-\theta_{m}^{(0)}$ (markers), combined with the evaluation of Eq.
(14) with $\rho=1000/(2\pi)$ (red solid lines) for $N=1000$ with $\gamma=0.5$
(upper panel) and $\gamma=1.5$ (lower panel). Each of the curves is normalized
to unit area over the plotted range. Blue lines connecting the markers serve
as a guide to the eye.

For $\gamma=0.5$ the expected spreading width has a value $\Gamma\approx
15.8$, such that the resulting Lorentzian shape has a width exceeding the
width $2\pi$ of the spectrum, leading to a slight mismatch with the observed
results. For $\gamma=1.5$ the spreading width evaluates to $\Gamma\approx
1.6\times 10^{-3}$, which is comparable to the mean level spacing
$\rho^{-1}\approx 6.3\times 10^{-4}$. Eq. (4) relies on a continuum
approximation for the density of states [30]. Presumably due to this, the
agreement between our data and the expected results is not perfect, although
qualitative agreement, and a pronounced difference between the structures of
the eigenstates for $\gamma=0.5$ and $\gamma=1.5$ is clearly seen. To our
knowledge, this characteristic of eigenstates has not been previously
investigated for unitary matrix ensembles.

##  V Conclusions and outlook

We have proposed a unitary (circular) analogue of the Rosenzweig-Porter (RP)
ensemble, defined as the outcome of a Dyson Brownian motion process. We have
numerically verified that some key statistics of both the eigenvalues and the
eigenstates of the circular analogue match the behaviour of the RP ensemble.
The circular analogue, therefore, can serve as a phenomenological model for
the level statistics and fractality of eigenstates as observed across the
many-body localization transition for periodically driven systems.

Motivated by, among other observations, the suggestion that eigenstates near
the many-body localization transition are multifractal [65, 66, 67], several
generalizations of the RP ensemble have been proposed [29, 43, 44, 45, 46]. It
would be interesting to see if similar generalization could be constructed for
the circular RP ensemble, next to other Floquet models with multifractal
eigenstates [47, 48, 49, 50]. This might for example be achievable by
considering stochastic processes with correlated increments, which has been
initialized for the RP ensemble recently [68]. Second, the circular RP
ensemble could potentially be of value in studies on random quantum circuits
[69] as a non-maximally random building block, analog to e.g. the proposal
made in Ref. [70].

###### Acknowledgements.

We thank Ivan M. Khaymovich for the careful reading of the manuscript and
providing us with many useful comments, as well as Lea F. Santos for useful
discussions. YB acknowledges support from the Israel Science Foundation
(grants No. 218/19 and 527/19). WB acknowledges support from the Kreitman
School of Advanced Graduate Studies at Ben-Gurion University.

## References

  * Basko  _et al._ [2006] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Metal-insulator transition in a weakly interacting many-electron system with localized single-particle states, Ann. Phys. 321, 1126 (2006).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting Electrons in Disordered Wires: Anderson Localization and Low-$T$ Transport, Phys. Rev. Lett. 95, 206603 (2005).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Adv. Phys. 65, 239 (2016).
  * Huse  _et al._ [2014] D. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, Phys. Rev. B 90, 174202 (2014).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-Body Localization and Thermalization in Quantum Statistical Mechanics, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, C.R. Phys. 19, 498 (2018).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Altshuler  _et al._ [1997] B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov, Quasiparticle Lifetime in a Finite System: A Nonperturbative Approach, Phys. Rev. Lett. 78, 2803 (1997).
  * Monthus and Garel [2010] C. Monthus and T. Garel, Many-body localization transition in a lattice model of interacting fermions: Statistics of renormalized hoppings in configuration space, Phys. Rev. B 81, 134202 (2010).
  * Serbyn  _et al._ [2013] M. Serbyn, Z. Papić, and D. A. Abanin, Local Conservation Laws and the Structure of the Many-Body Localized States, Phys. Rev. Lett. 111, 127201 (2013).
  * Ros  _et al._ [2015] V. Ros, M. Müller, and A. Scardicchio, Integrals of motion in the many-body localized phase, Nucl. Phys. B 891, 420 (2015).
  * Pietracaprina  _et al._ [2016] F. Pietracaprina, V. Ros, and A. Scardicchio, Forward approximation as a mean-field approximation for the Anderson and many-body localization transitions, Phys. Rev. B 93, 054201 (2016).
  * Roy and Logan [2019] S. Roy and D. E. Logan, Strong ergodicity breaking due to local constraints in a quantum system, SciPost Phys. 7, 042 (2019).
  * Roy and Lazarides [2020] S. Roy and A. Lazarides, Strong ergodicity breaking due to local constraints in a quantum system, Phys. Rev. Research 2, 023159 (2020).
  * Bera  _et al._ [2015] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-Body Localization Characterized from a One-Particle Perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * Villalonga  _et al._ [2018] B. Villalonga, X. Yu, D. J. Luitz, and B. K. Clark, Exploring one-particle orbitals in large many-body localized systems, Phys. Rev. B 97, 104406 (2018).
  * Buijsman  _et al._ [2018] W. Buijsman, V. Gritsev, and V. Cheianov, Many-body localization in the Fock space of natural orbitals, Scipost Phys. 4, 038 (2018).
  * Hopjan and Heidrich-Meisner [2020] M. Hopjan and F. Heidrich-Meisner, Many-body localization from a one-particle perspective in the disordered one-dimensional Bose-Hubbard model, Phys. Rev. A 101, 063617 (2020).
  * Roy and Logan [2021] S. Roy and D. E. Logan, Fock-space anatomy of eigenstates across the many-body localization transition, Phys. Rev. B 104, 174201 (2021).
  * Altland and Micklitz [2017] A. Altland and T. Micklitz, Field Theory Approach to Many-Body Localization, Phys. Rev. Lett. , 127202 (2017).
  * Roy and Logan [2020] S. Roy and D. E. Logan, Fock-space correlations and the origins of many-body localization, Phys. Rev. B 101, 134202 (2020).
  * De Luca  _et al._ [2014] A. De Luca, B. L. Altshuler, V. E. Kravtsov, and A. Scardicchio, Anderson Localization on the Bethe Lattice: Nonergodicity of Extended States, Phys. Rev. Lett. 113, 046806 (2014).
  * Kravtsov  _et al._ [2015] V. E. Kravtsov, I. M. Khaymovich, E. Cuevas, and M. Amini, A random matrix model with localization and ergodic transitions, New J. Phys. 17, 122002 (2015).
  * Tikhonov and Mirlin [2021] K. S. Tikhonov and A. D. Mirlin, From Anderson localization on random regular graphs to many-body localization, Ann. Phys. 435, 168525 (2021).
  * Rosenzweig and Porter [1960] N. Rosenzweig and C. E. Porter, “Repulsion of Energy Levels” in Complex Atomic Spectra, Phys. Rev. 120, 1698 (1960).
  * Facoetti  _et al._ [2016] D. Facoetti, P. Vivo, and G. Biroli, From non-ergodic eigenvectors to local resolvent statistics and back: A random matrix perspective, EuroPhys. Lett. 115, 47003 (2016).
  * Truong and Ossipov [2016] K. Truong and A. Ossipov, Eigenvectors under a generic perturbation: Non-perturbative results from the random matrix approach, New J. Phys. 116, 37002 (2016).
  * Amini [2017] M. Amini, Spread of wave packets in disordered hierarchical lattices, Europhys. Lett. 117, 30003 (2017).
  * Monthus [2017] C. Monthus, Multifractality of eigenstates in the delocalized non-ergodic phase of some random matrix models: Wigner–Weisskopf approach, J. Phys. A: Math. Theor. 50, 295101 (2017).
  * Bogomolny and Sieber [2018] E. Bogomolny and M. Sieber, Eigenfunction distribution for the Rosenzweig-Porter model, Phys. Rev. E 98, 032139 (2018).
  * De Tomasi  _et al._ [2019] G. De Tomasi, M. Amini, S. Bera, I. M. Khaymovich, and V. E. Kravtsov, Survival probability in Generalized Rosenzweig-Porter random matrix ensemble, SciPost Phys. 6, 014 (2019).
  * von Soosten and Warzel [2019] P. von Soosten and S. Warzel, Non-ergodic delocalization in the Rosenzweig-Porter model, Lett. Math. Phys. 109, 905 (2019).
  * Pino  _et al._ [2019] M. Pino, J. Tabanera, and P. Serna, From ergodic to non-ergodic chaos in Rosenzweig-Porter model, J. Phys. A: Math. Theor. 52, 475101 (2019).
  * Berkovits [2020] R. Berkovits, Super-Poissonian behavior of the Rosenzweig-Porter model in the nonergodic extended regime, Phys. Rev. B 102, 165140 (2020).
  * Tarzia [2020] M. Tarzia, Many-body localization transition in Hilbert space, Phys. Rev. B 102, 014208 (2020).
  * Ponte  _et al._ [2015] P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin, Many-Body Localization in Periodically Driven Systems, Phys. Rev. Lett. 114, 140401 (2015).
  * Lazarides  _et al._ [2015] A. Lazarides, A. Das, and R. Moessner, Fate of Many-Body Localization Under Periodic Driving, Phys. Rev. Lett. 115, 030402 (2015).
  * Abanin  _et al._ [2016] D. A. Abanin, W. de Roeck, and F. Huveneers, Theory of many-body localization in periodically driven systems, Ann. Phys. 372, 1 (2016).
  * Zhang  _et al._ [2016] L. Zhang, V. Khemani, and D. A. Huse, A Floquet model for the many-body localization transition, Phys. Rev. B 94, 224202 (2016).
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, Periodically driving a many-body localized quantum system, Nat. Phys. 13, 460 (2017).
  * Roy  _et al._ [2021] S. Roy, R. Moessner, and A. Lazarides, How periodic driving stabilizes and destabilizes Anderson localization on random trees, Phys. Rev. B 103, L100204 (2021).
  * Dyson [1962] F. J. Dyson, A Brownian-Motion Model for the Eigenvalues of a Random Matrix, J. Math. Phys. 3, 1196 (1962).
  * Kravtsov  _et al._ [2020] V. E. Kravtsov, I. M. Khaymovich, B. L. Altshuler, and L. B. Ioffe, Localization transition on the Random Regular Graph as an unstable tricritical point in a log-normal Rosenzweig-Porter random matrix ensemble, arXiv:2002.02979 (2020).
  * Khaymovich  _et al._ [2020] I. M. Khaymovich, V. E. Kravtsov, B. L. Altshuler, and L. B. Ioffe, Fragile extended phases in the log-normal Rosenzweig-Porter model, Phys. Rev. Research 2, 043346 (2020).
  * Biroli and Tarzia [2021] G. Biroli and M. Tarzia, Lévy-Rosenzweig-Porter random matrix ensemble, Phys. Rev. B 103, 104205 (2021).
  * Khaymovich and Kratsov [2021] I. M. Khaymovich and V. E. Kratsov, Dynamical phases in a “multifractal” Rosenzweig-Porter model, SciPost Phys. 11, 45 (2021).
  * Bandyopadhyay  _et al._ [2010] J. N. Bandyopadhyay, J. Wang, and J. Gong, Generating a fractal butterfly Floquet spectrum in a class of driven $SU(2)$ systems: Eigenstate statistics, Phys. Rev. E 81, 066212 (2010).
  * Ray  _et al._ [2018] S. Ray, A. Ghosh, and S. Sinha, Drive-induced delocalization in the Aubry-André model, Phys. Rev. E 97, 010101 (2018).
  * Roy  _et al._ [2018] S. Roy, I. M. Khaymovich, A. Das, and R. Moessner, Multifractality without fine-tuning in a Floquet quasiperiodic chain, SciPost Phys. 4, 025 (2018).
  * Sarkar  _et al._ [2021] M. Sarkar, R. Ghosh, A. Sen, and K. Sengupta, Mobility edge and multifractality in a periodically driven Aubry-André model, Phys. Rev. B 103, 184309 (2021).
  * Mehta [2004] M. L. Mehta,  _Random Matrices_ , 3rd ed., Pure and Applied Mathematics, Vol. 142 (Elsvier, New York, 2004).
  * Forrester [2010] P. J. Forrester,  _Log-Gases and Random Matrices_ , London Mathematical Society Monographs, Vol. 34 (Princeton University Press, Princeton and Oxford, 2010).
  * Altland  _et al._ [1997] A. Altland, M. Janssen, and B. Shapiro, Perturbation theory for the Rosenzweig-Porter matrix model, Phys. Rev. E 56, 1471 (1997).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the Ratio of Consecutive Level Spacings in Random Matrix Ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * Wigner [1955] E. P. Wigner, Characteristic Vectors of Bordered Matrices With Infinite Dimensions, Ann. Math. 62, 548 (1955).
  * Jacquod and Shepelyansky [1995] P. Jacquod and D. L. Shepelyansky, Hidden Breit-Wigner Distribution and Other Properties of Random Matrices with Preferential Basis, Phys. Rev. Lett. 75, 3501 (1995).
  * Fyodorov  _et al._ [1996] Y. V. Fyodorov, O. A. Chubykalo, F. M. Izrailev, and G. Casati, Wigner Random Banded Matrices with Sparse Structure: Local Spectral Density of States, Phys. Rev. Lett. 76, 1603 (1996).
  * Borgonovi  _et al._ [2016] F. Borgonovi, F. M. Izrailev, L. F. Santos, and V. G. Zelevinsky, Quantum chaos and thermalization in isolated systems of interacting particles, Phys. Rep. 626, 1 (2016).
  * Cuevas and Kravtsov [2007] E. Cuevas and V. E. Kravtsov, Two-eigenfunction correlation in a multifractal metal and insulator, Phys. Rev. B 76, 235119 (2007).
  * Allez  _et al._ [2015] R. Allez, J. Bun, and J.-P. Bouchaud, The eigenvectors of Gaussian matrices with an external source, arXiv:1412.7108v4 (2015).
  * Bourgade and Yau [2017] P. Bourgade and H.-T. Yau, The eigenvector moment flow and local quantum unique ergodicity, Commun. Math. Phys. 350, 231 (2017).
  * Note [1] The set $U_{S}(N)$ of symmetric unitary matrices of dimension $N$ is isomorphic to $U(N)/O(N)$, where $U(N)$ denotes the set of imaginary unitary matrices of dimension $N$, and $O(N)$ denote the set of real orthogonal matrices of dimension $N$. For a proof, see e.g. Proposition 2.2.4 of Ref. [52].
  * Forrester and Nagao [1998] P. J. Forrester and T. Nagao, Correlations for the circular Dyson Brownian motion model with Poisson initial conditions, Nucl. Phys. B 532, 733 (1998).
  * Serbyn  _et al._ [2017] M. Serbyn, Z. Papić, and D. A. Abanin, Thouless energy and multifractality across the many-body localization transition, Phys. Rev. B 96, 104201 (2017).
  * Macé  _et al._ [2019] N. Macé, F. Alet, and N. Laflorencie, Multifractal Scalings Across the Many-Body Localization Transition, Phys. Rev. Lett. 123, 180601 (2019).
  * Luitz  _et al._ [2020] D. J. Luitz, I. M. Khaymovich, and Y. Bar Lev, Multifractality and its role in anomalous transport in the disordered XXZ spin-chain, SciPost Phys. Core 2, 006 (2020).
  * Kutlin and Khaymovich [2021] A. G. Kutlin and I. M. Khaymovich, Emergent fractal phase in energy stratified random models, SciPost Phys. 11, 101 (2021).
  * Chan  _et al._ [2018a] A. Chan, A. De Luca, and J. T. Chalker, Solution of a Minimal Model for Many-Body Quantum Chaos, Phys. Rev. X 8, 041019 (2018a).
  * Chan  _et al._ [2018b] A. Chan, A. De Luca, and J. T. Chalker, Spectral Statistics in Spatially Extended Chaotic Quantum Many-Body Systems, Phys. Rev. Lett. 121, 060601 (2018b).