# Quantum information theoretic measures to distinguish fermionized bosons
from non-interacting fermions

Barnali Chakrabarti1,2, Arnaldo Gammal2, N D Chavda3, Mantile Leslie Lekala4
1Department of Physics, Presidency University, 86/1 College Street, Kolkata
700073, India.
2 Instituto de Física, Universidade de São Paulo, CEP 05580-090, São Paulo,
Brazil
3 Department of Applied Physics, Faculty of Technology and Engineering,
The Maharaja Sayajirao University of Baroda, Vadodara-390001, India
4 Department of Physics, University of South Africa, P. O. Box 392, Pretoria
0003, South Africa.

###### Abstract

We study the dynamical fermionization of strongly interacting one-dimensional
bosons in Tonks-Girardeau limit by solving the time dependent many-boson
Schrödinger equation numerically exactly. We establish that the one-body
momentum distribution approaches the ideal Fermi gas distribution at the time
of dynamical fermionization. The analysis is further complemented by the
measures on two-body level. Investigation on two-body momentum distribution,
two-body local and non-local correlation clearly distinguish the fermionized
bosons from non-interacting fermions. The magnitude of distinguishablity
between the two systems is further discussed employing suitable measures of
information theory, i.e., the well known Kullback-Leibler relative entropy and
the Jensen-Shannon divergence entropy. We also observe very rich structure in
the higher-body density for strongly correlated bosons whereas non-interacting
fermions do not possess any higher order correlation beyond two-body.

Fermionization, Glauber correlation function, quantum information

##  I Introduction

The one-dimensional (1D) gas of impenetrable bosons termed as Tonks-Girardeau
(TG) gas corresponds to infinitely strong repulsive bosons in the Lieb-Liniger
(LL) model [1, 2], the bosons interact through a contact potential. The
physical properties of 1D bosons in the TG limit coincide with those of ideal
Fermi gas in the same external potential. The properties determined by the
square of the wave function such as the density profile and density-density
correlation become identical with those of non-interacting fermions. It leads
to the interesting manifestation of fermionization when the strongly
interacting bosons escape their spatial overlap [3, 4, 5, 6, 7]. However one-
body density matrix and the momentum distribution are typically different from
those of Fermi gas due to strong phase correlation originated from the bosonic
statistics at the TG limit.

The dynamics of harmonically confined TG gas exhibits more interesting
scenario ”dynamical fermionization” during expansion on sudden release of the
trap–the momentum distribution of TG gas evolves from bosonic to fermionic
character. Some aspects of dynamical evolution of the TG gas has been
theoretically studied [8, 9, 10, 11]. Experimentally, optical confinement
offers quasi-1D geometry and dynamical fermionization has recently been
observed in the measure of asymptotic momentum distribution [12].

Despite of well established understanding of dynamical fermionization, the
topic still remains challenging as no work is reported on higher-body density
and correlation measures. The existing calculations are based on the one-body
level only. The aim of this work is to present some results beyond one-body
measures in real as well as momentum space and to understand whether
fermionized bosons can really be distinguished from ideal Fermi gas. The
dynamical evolution of quasi-1D harmonically trapped strongly interacting
bosons in the TG limit is investigated in this work utilizing an ab initio
many-body technique. We solve the many bosons Schrödinger equation by using
the multiconfigurational time dependent Hartree method for indistinguishable
particles (MCTDH-X) which calculates numerically exact many-body wave function
[13, 14, 15, 16, 17, 18, 19]. Measures of many-body densities and correlation
functions are paramount interest to understand many-body physics. In such
strongly interacting regime we need very high orbitals to obtain the
convergence in the many-body wave function. In the present work, dynamical
evolution of $N=5$ fermionized bosons is studied utilizing $M=24$ orbitals in
the computation of many-body wave function.

Our work is outlined as : i) To explore the dynamical evolution of the
fermionized bosons through the measures of $x$-space, $k$-space density and
correlation in momentum space beyond one-body. ii) To apply information
theoretic measures of relative entropy like, non-symmetric Kullback-Leibler,
symmetric Kullback-Leibler and the Jensen-Shannon divergence entropy both for
one-body and two-body densities. The main aim is to justify that one-body
measures of density, correlation, relative entropies all are inadequate to
understand the difference between the fermions and strongly interacting bosons
which approach fermionization limit dynamically.

Our observations are as follows: i) As expected, dynamics of one-body
densities in $x$-space are identical between strongly interacting bosons and
those of non interacting fermions. ii) One-body $k$-space density exhibits
that strongly interacting bosons attain fermionic character at a particular
time, known as time for dynamical fermionization. iii) Measures of different
relative entropies using one-body $k$-density asymptotically approach to zero
at the time of fermionization. iv) Dynamics of the two-body $x$-space
densities also exhibit the identical evolution for the two systems and can not
extract any additional physics. v) However two-body $k$-space density of
strongly interacting bosons dynamically tends to approach the fermionic two-
body $k$-space density, but remains distinctly different even at the time of
fermionization. vi) Measures of local and non-local two-body correlation in
the momentum space quantitatively estimate the difference between the two
systems. vii) Measures of relative entropies utilizing two-body density in
$k$-space also clearly establish the mutual information between the
fermionized bosons and fermions. viii) We additionally analyze higher-body
coherence. Fermions do not possess any higher order correlation, whereas
strongly interacting bosons exhibit rich structure in the three-body and four-
body coherence.

All the above observations uniquely confirm that the strongly interacting
bosons at TG limit dynamically acquire fermionic character but can be
distinguished. The conclusion made with measures in one-body level is
insufficient. Two-body correlation measures in $k$-space are the key
quantities to distinguish how the ferminozed bosons differ from ideal fermions
qualitatively, whereas information theoretic measures of relative entropies
quantitatively distinguish the two systems.

The paper is organized in the following way. Section II briefly describes the
theoretical framework. We present dynamical measures in one-, two- and higher-
body levels in Section III and conclusion is made in section IV.

##  II Theoretical Framework

In this section we introduce a brief description of the multiconfigurtaional
time dependent Hartree for bosons (MCTDHB) which is used to explore the
dynamics. We also introduce the basic measures which are utilized to
understand the many-body dynamics.

###  II.1  Many-body Hamiltonian, wave function

The dynamical evolution of $N$-bosons is governed by the time dependent
Schrödinger equation

$H|\psi(t)\rangle=i\frac{\partial}{\partial t}|\psi(t)\rangle$ (1)

The full many-body Hamiltonian is

$\hat{H}(x_{1},x_{2},\dots
x_{N})=\sum_{i=1}^{N}\hat{h}(x_{i})+\sum_{i<j=1}^{N}\hat{W}(x_{i}-x_{j})$ (2)

We used dimensionless unit by dividing the Hamiltonian by
$\frac{\hbar^{2}}{mL^{2}}$, where $m$ is the mass of a boson and $L$ is the
length scale. Where, $h(x)=-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+V(x)$ is the one-body Hamiltonian. $V(x)$ is the trapping potential
which is a harmonic potential for the present work. The interparticle
interaction potential is $W(x_{i}-x_{j})=\lambda\delta(x_{i}-x_{j})$,
$\lambda$ is the interaction strength.

In order to compute the time-evolution of the many-body Schrödinger equation
we use the MCTDHB, which employs variational time-adaptive orbitals in the
expansion of the many-boson wave function. With increase in number of orbitals
it converges to the exact many-body Schrödinger results. In MCTDHB method,
time-dependent orbitals are used in the expansion of the field operator

$\hat{\Psi}(x)=\sum^{M}_{j=1}\hat{b}_{j}(t)\phi_{j}(x,t).$ (3)

Here $\\{\phi_{j}(x,t)\\}$ is a complete orthonormal set of orbitals.
$\hat{b}_{j}$ annihilates a boson in $\phi_{j}$. This is in contrast with the
full many-body Hamiltonian using time-independent orbitals, which require a
large Fock space for convergence even for a few bosons. Introducing time-
dependent orbitals, a faithful convergence can be achieved with a truncated
number of orbitals. The orbitals are determined from the time-dependent
variational principle [20, 21, 22, 23]. When $N$ bosons are distributed over
$M$ time-adaptive orbitals, the time-dependent many-boson function is taken as
a linear combination of time dependent permanents $|\bar{n};t\rangle$

$|\psi(t)\rangle=\sum_{\bar{n}}C_{\bar{n}}(t)|\bar{n};t\rangle,$ (4)

In Eq. (4), the $\\{C_{\bar{n}}(t)\\}$ are the time-dependent expansion
coefficients. The vector $\vec{n}=(n_{1},n_{2},\dots,n_{M})$ represents the
occupation of the orbitals and $n_{1}+n_{2}+\dots+n_{M}=N$ which preserves the
total number of particles. In second quantisation representation, the
permanents are given as

$|\bar{n};t\rangle=\prod_{i=1}^{M}\left(\frac{\left(b_{i}^{\dagger}(t)\right)^{n_{i}}}{\sqrt{n_{i}!}}\right)|vac\rangle.$
(5)

In Eq.(4), the summation runs over all possible configurations and in the
limit of $M\rightarrow\infty$, the set of permanents spans the complete
Hilbert space and the expansion becomes exact. However allowing time dependent
permanents, we use mush shorter expansion which leads to significant
computational advantage.

It is to be noted that both the coefficients $\\{C_{\bar{n}}(t)\\}$ and the
orbitals $\\{\phi_{k}(x,t)\\}$ which comprise the permanents are independent
parameters and are determined by the time-dependent variational principle. For
the equations of motion governing the time evolution of $\\{C_{\bar{n}}(t)\\}$
and $\\{\phi_{k}(x,t)\\}$ we follow the variational principle based on the
Lagrangian formulation [20]. Substituting the many-body ansatz into the
functional action of the time dependent many-body Schrödinger equation leads
to

$\begin{split}S\left[\\{C_{\bar{n}}(t)\\},\\{\phi_{k}(x,t)\\}\right]=\int
dt\\{\langle\Psi(t)|\hat{H}-i\frac{\partial}{\partial t}|\Psi(t)\rangle\\\
-\sum^{M}_{k,j=1}\mu_{kj}(t)\left[\langle\phi_{k}|\phi_{j}\rangle-\delta_{kj}\right]\\}\end{split}$
(6)

Next, stationarity of the action with respect to the independent variations of
$\\{C_{\bar{n}}(t)\\}$ and $\\{\phi_{j}(x,t)\\}$ are required.
$\\{\mu_{kj}(t)\\}$ are time-dependent Lagrange multipliers which guarantee
the orthonormality of the orbitals during time propagation.

This is to be noted that the use of optimized time dependent orbitals leads to
very fast convergence in the simulation compared to the many-body Schrödinger
equation with time-independent orbitals, a given degree of accuracy is reached
with much shorter expansion [24, 25]. We also emphasize that MCTDHB is more
accurate than exact diagonalization which uses the finite basis and is not
optimized. Whereas in MCTDHB, as we use a time adaptive many-body basis set,
it can dynamically follow the building correlation due to inter-particle
interaction [16, 15, 26, 27]. It is already established as a very efficient
many-body method and used for different trap geometry and range of inter-
particle interaction [28, 29, 30, 31]. The MCTDHB method and the algorithm has
been cast into a software package [32, 33].

###  II.2  Quantities of interest

####  II.2.1 One-body and higher-order densities

(i)The reduced one-body density matrix in coordinate space is defined as

$\begin{split}\rho^{(1)}(x_{1}^{\prime}|x_{1};t)=N\int
dx_{2}\,dx_{3}...dx_{N}\,\psi^{*}(x_{1}^{\prime},x_{2},\dots,x_{N};t)\\\
\psi(x_{1},x_{2},\dots,x_{N};t).\end{split}$ (7)

Its diagonal gives the one-body density $\rho(x,t)$ defined as

$\begin{split}\rho(x;t)=\rho^{(1)}(x_{1}^{\prime}=x|x_{1}=x;t)\end{split}$ (8)

(ii) The $p$-th order reduced density matrix in coordinate space is defined by

$\begin{split}\rho^{(p)}(x_{1}^{\prime},\dots,x_{p}^{\prime}|x_{1},\dots,x_{p};t)=\frac{N!}{(N-p)!}\int
dx_{p+1}\,...dx_{N}\\\
\psi^{*}(x_{1}^{\prime},\dots,x_{p}^{\prime},x_{p+1},\dots,x_{N};t)\psi(x_{1},\dots,x_{p},x_{p+1}\dots,x_{N};t).\end{split}$
(9)

The diagonal of $p$-body density can be represented as

$\begin{split}\rho^{(p)}(x_{1},\dots,x_{p};t)=\langle\psi(t)|\hat{\psi}^{\dagger}(x_{1})\dots\hat{\psi}^{\dagger}(x_{p})\hat{\psi}(x_{p})\\\
\dots\hat{\psi}(x_{1})|\psi(t)\rangle\end{split}$ (10)

It provides the $p$-particle density distribution at time t.
To calculate the density in momentum space one needs to follow the
prescription provided in the Ref [34]. Introducing $p$ $D$ dimensional Fourier
transform to $r_{i}$ and $r_{i}^{\prime}$ one arrives to

$\begin{split}\rho^{(p)}(k_{1}.....k_{p}|k_{1}^{\prime}.....k_{p}^{\prime};t)=\\\
\sum_{i}n_{i}^{(p)}(t)\alpha_{i}^{(p)}(k_{1}.....k_{p},t){\alpha_{i}^{(p)}}^{*}(k_{1}^{\prime}.....k_{p}^{\prime},t)\end{split}$
(11)

where $n_{i}^{(p)}(t)$ is the $i^{th}$ eigenvalue of the $p^{th}$ order
reduced density matrix and $\alpha_{i}^{(p)}(k_{1}.....k_{p},t)$ the
corresponding eigenfunction. The eigenfunctions are known as natural $p$
functions and the eigenvalues as natural occupations.

####  II.2.2 Glauber correlation function

The normalized $p$-th order Glauber correlation function would be the most
important quantity to measure the spatial coherence. It is defined as

$\begin{split}g^{(p)}(x_{1}^{\prime},\dots,x_{p}^{\prime},x_{1},\dots,x_{p};t)=\\\
\frac{\rho^{(p)}(x_{1},\dots,x_{p}|x_{1}^{\prime},\dots,x_{p}^{\prime};t)}{\sqrt{\Pi_{i=1}^{p}\rho^{(1)}(x_{i}|x_{i};t)\rho^{(1)}(x_{i}^{\prime}|x_{i}^{\prime};t)}}\end{split}$
(12)

The diagonal of
$g^{(p)}(x_{1}^{\prime},\dots,x_{p}^{\prime},x_{1},\dots,x_{p};t)$ gives a
measure of $p$-th order coherence and is calculated as

$g^{(p)}(x_{1},\dots,x_{p};t)=\frac{\rho^{(p)}(x_{1},\dots,x_{p};t)}{\Pi_{i=1}^{p}|\rho^{(1)}(x_{i})|}$
(13)

If $|g^{(p)}(x_{1},\dots,x_{p};t)|=1$, the system is fully coherent, otherwise
the state is only partially coherent. When $g^{(p)}(x_{1},\dots,x_{p};t)>1$,
the detection probabilities at positions $x_{1},\dots,x_{p}$ are correlated
and $g^{(p)}(x_{1},\dots,x_{p};t)<1$ are anti-correlated. Eq. (13) is further
utilized to calculate the local two-particle correlation $g^{(2)}(0,0)$ and
non-local two-particle correlation $g^{(2)}(0,x)$.

The $p^{th}$ order correlation function in momentum space is given by

$\begin{split}g^{(p)}(k_{1}^{\prime},\dots,k_{p}^{\prime},k_{1},\dots,k_{p};t)=\\\
\frac{\rho^{(p)}(k_{1},\dots,k_{p}|k_{1}^{\prime},\dots,k_{p}^{\prime};t)}{\sqrt{\Pi_{i=1}^{p}\rho^{(1)}(k_{i}|k_{i};t)\rho^{(1)}(k_{i}^{\prime}|k_{i}^{\prime};t)}}\end{split}$
(14)

The diagonal of
$g^{(p)}(k_{1}^{\prime},\dots,k_{p}^{\prime},k_{1},\dots,k_{p};t)$ measures
the $p$-th order coherence in the momentum space. For values
$g^{(p)}(k_{1},\dots,k_{p};t)>1$, the detection probabilities at positions
$k_{1},\dots,k_{p}$ are correlated and $g^{(p)}(k_{1},\dots,k_{p};t)<1$
corresponds to loss of coherence and are anti-correlated. Eq. (14) is further
utilized to calculate the local two-particle correlation $g^{(2)}(0,0)$ and
non-local two-particle correlation $g^{(2)}(0,k)$.

####  II.2.3 Information distance measures

Shannon entropy and Fisher information are considered as the key measures to
exhibit higher order characteristics in position and momentum space density
distribution and has been extensively utilized in atomic systems [35, 36, 37,
38, 39]. However to measure the difference between two probability
distributions over the same variable, a measure called Kullback-Leibler (KL)
relative entropy is the ideal quantity [40, 41, 42]. Relative entropy provides
mutual information between two probability distributions of the same order and
defined in the same space. KL relative entropy in $x$\- space is measured as

$K(t)=\int{\rho_{1}^{(1)}(x,t)\ln\frac{\rho_{1}^{(1)}(x,t)}{\rho_{2}^{(1)}(x,t)}\;dx}$
(15)

Where $\rho_{1}^{(1)}$ and $\rho_{2}^{(1)}$ are the one-body densities for two
different systems described in the same space and they are calculated using
$\rho^{(1)}(x,t)$ =
$\langle\psi(t)|\hat{\psi}^{\dagger}(x)\hat{\psi}(x)|\psi(t)\rangle$.

Thus, the measure of $K$ can be interpreted as a measure of deviation of
$\rho_{1}^{(1)}$ from $\rho_{2}^{(1)}$ i.e, to estimate how the two systems
are identical or close or far apart. Thus for any distributions
$\rho_{1}^{(1)}$ and $\rho_{2}^{(1)}$, $K\geq 0$. However it is a non-
symmetric measure as it depends on which distribution is considered as
‘reference’ and which is considered as ‘comparison’ distribution. The
symmetrized Kullback distance $SK$ is defined as [43]

$\begin{split}SK(t)=\int{\rho_{1}^{(1)}(x,t)\ln\frac{\rho_{1}^{(1)}(x,t)}{\rho_{2}^{(1)}(x,t)}\;dx}\\\
+\int{\rho_{2}^{(1)}(x,t)\ln\frac{\rho_{2}^{(1)}(x,t)}{\rho_{1}^{(1)}(x,t)}\;dx}\end{split}$
(16)

The physical meaning of $SK$ distance is very clear. $SK$ is zero for two
identical species and approach to a large value as the difference between the
one-body densities of two systems $\rho_{1}^{(1)}$ and $\rho_{2}^{(1)}$
increases. The another symmetrized measure of relative entropy is the Jensen
Shannon divergence entropy as [44, 45]

$\begin{split}J(t)=-\int\left(\frac{\rho_{1}^{(1)}(x,t)+\rho_{2}^{(1)}(x,t)}{2}\right)\times\\\
\ln\left(\frac{\rho_{1}^{(1)}(x,t)+\rho_{2}^{(1)}(x,t)}{2}\right)dx\\\
+\frac{1}{2}\int{\rho_{1}^{(1)}(x,t)\ln\rho_{1}^{(1)}(x,t)\;dx}\\\
+\frac{1}{2}\int{\rho_{2}^{(1)}(x,t)\ln\rho_{2}^{(1)}(x,t)\;dx}\end{split}$
(17)

Both the measures of $SK$ and $J$ are symmetrized and compare two probability
distributions regardless of which distribution is considered as the
‘reference’ and which is considered as the ‘comparison’ distribution. The
corresponding definitions of $K$, $SK$ and $J$ measures in $k$-space are

$K(t)=\int{\rho_{1}^{(1)}(k,t)\ln\frac{\rho_{1}^{(1)}(k,t)}{\rho_{2}^{(1)}(k,t)}\;dk}$
(18)

$\begin{split}SK(t)=\int{\rho_{1}^{(1)}(k,t)\ln\frac{\rho_{1}^{(1)}(k,t)}{\rho_{2}^{(1)}(k,t)}dk}\\\
+\int{\rho_{2}^{(1)}(k,t)\ln\frac{\rho_{2}^{(1)}(k,t)}{\rho_{1}^{(1)}(k,t)}\;dk}\end{split}$
(19)

$\begin{split}J(t)=-\int\left(\frac{\rho_{1}^{(1)}(k,t)+\rho_{2}^{(1)}(k,t)}{2}\right)\times\\\
\ln\left(\frac{\rho_{1}^{(1)}(k,t)+\rho_{2}^{(1)}(k,t)}{2}\right)dk\\\
+\frac{1}{2}\int{\rho_{1}^{(1)}(k,t)\ln\rho_{1}^{(1)}(k,t)\;dk}\\\
+\frac{1}{2}\int{\rho_{2}^{(1)}(k,t)\ln\rho_{2}^{(1)}(k,t)\;dk},\end{split}$
(20)

where one-body density in $k$-space $\rho^{(1)}(k,t)$ is calculated using
$\rho^{(1)}(k,t)$ =
$\langle\psi(t)|\hat{\psi}^{\dagger}(k)\hat{\psi}(k)|\psi(t)\rangle$.

However it is also interesting to have measures of the above three
distributions using two-body density in $x$-space $\rho^{(2)}(x_{1},x_{2},t)$
=
$\langle\psi(t)|\hat{\psi}^{\dagger}(x_{1})\hat{\psi}^{\dagger}(x_{2})\hat{\psi}(x_{1})\hat{\psi}(x_{2})|\psi(t)\rangle$.

The corresponding measures of $K$, $SK$ and $J$ in $x$ -space are defined
through

$K(t)=\int{\rho_{1}^{(2)}(x_{1},x_{2},t)\ln\frac{\rho_{1}^{(2)}(x_{1},x_{2},t)}{\rho_{2}^{(2)}(x_{1},x_{2},t)}\;dx_{1}\;dx_{2}},$
(21)

where $\rho_{1}^{(2)}$ and $\rho_{2}^{(2)}$ correspond to two-body densities
of the two different systems but described in the same space. It will
facilitate to understand whether the two systems are close in the measure of
two-body perspective. The corresponding definitions of $SK$ and $J$ in
$x$-space are as follows :

$\begin{split}SK(t)=\int{\rho_{1}^{(2)}(x_{1},x_{2},t)\ln\frac{\rho_{1}^{(2)}(x_{1},x_{2},t)}{\rho_{2}^{(2)}(x_{1},x_{2},t)}dx_{1}dx_{2}}\\\
+\int{\rho_{2}^{(2)}(x_{1},x_{2},t)\ln\frac{\rho_{2}^{(2)}(x_{1},x_{2},t)}{\rho_{1}^{(2)}(x_{1},x_{2},t)}\;dx_{1}\;dx_{2}}\end{split}$
(22)

$\begin{split}J(t)=-\int\left(\frac{\rho_{1}^{(2)}(x_{1},x_{2},t)+\rho_{2}^{(2)}(x_{1},x_{2},t)}{2}\right)\times\\\
\ln\left(\frac{\rho_{1}^{(2)}(x_{1},x_{2},t)+\rho_{2}^{(2)}(x_{1},x_{2},t)}{2}\right)\;dx_{1}\;dx_{2}\\\
+\frac{1}{2}\int\rho_{1}^{(2)}(x_{1},x_{2},t)\ln\rho_{1}^{(2)}(x_{1},x_{2},t)\;dx_{1}\;dx_{2}\\\
+\frac{1}{2}\int\rho_{2}^{(2)}(x_{1},x_{2},t)\ln\rho_{2}^{(2)}(x_{1},x_{2},t)\;dx_{1}\;dx_{2}\end{split}$
(23)

Defining two-body density in $k$-space as $\rho^{(2)}(k_{1},k_{2},t)$ =
$\langle\psi(t)|\hat{\psi}^{\dagger}(k_{1})\hat{\psi}^{\dagger}(k_{2})\hat{\psi}(k_{1})\hat{\psi}(k_{2})|\psi(t)\rangle$
we can have the measures of $S$, $SK$ and $J$ in momentum space from two-body
perspective. They are calculated as

$K(t)=\int{\rho_{1}^{(2)}(k_{1},k_{2},t)\ln\frac{\rho_{1}^{(2)}(k_{1},k_{2},t)}{\rho_{2}^{(2)}(k_{1},k_{2},t)}\;dk_{1}\;dk_{2}}$
(24)

$\begin{split}SK(t)=\int{\rho_{1}^{(2)}(k_{1},k_{2},t)\ln\frac{\rho_{1}^{(2)}(k_{1},k_{2},t)}{\rho_{2}^{(2)}(k_{1},k_{2},t)}dk_{1}dk_{2}}\\\
+\int{\rho_{2}^{(2)}(k_{1},k_{2},t)\ln\frac{\rho_{2}^{(2)}(k_{1},k_{2},t)}{\rho_{1}^{(2)}(k_{1},k_{2},t)}\;dk_{1}\;dk_{2}}\end{split}$
(25)

$\begin{split}J(t)=-\int\left(\frac{\rho_{1}^{(2)}(k_{1},k_{2},t)+\rho_{2}^{(2)}(k_{1},k_{2},t)}{2}\right)\times\\\
\ln\left(\frac{\rho_{1}^{(2)}(k_{1},k_{2},t)+\rho_{2}^{(2)}(k_{1},k_{2},t)}{2}\right)\;dk_{1}\;dk_{2}\\\
+\frac{1}{2}\int\rho_{1}^{(2)}(k_{1},k_{2},t)\ln\rho_{1}^{(2)}(k_{1},k_{2},t)\;dk_{1}\;dk_{2}\\\
+\frac{1}{2}\int\rho_{2}^{(2)}(k_{1},k_{2},t)\ln\rho_{2}^{(2)}(k_{1},k_{2},t)\;dk_{1}\;dk_{2}\end{split}$
(26)

Figure 1: Snapshots of the time evolution of reduced one-body density in the
$x$-space ($\rho^{(1)}(x,x^{\prime})$) during the expansion of $N=5$ strongly
interacting bosons with interaction strength $\lambda=25$ are presented (left
panel). Computation is done with $M=24$ orbitals. Comparison is made for the
dynamics of non interacting fermions in the right panel, computation is made
with $M=5$ orbitals. The dynamical evolution of the two systems is identical
in $x$-space. Figure 2: Snapshots of the time evolution of reduced one-body
density in the $k$-space ($\rho^{(1)}(k,k^{\prime})$) during the expansion of
$N=5$ strongly interacting bosons with interaction strength $\lambda=25$ are
presented for different times. Computation is done with $M=24$ orbitals.
Comparison is made with fermionic momentum density distribution (last panel in
the right side), computation is made with $M=5$ orbitals. At time $t=8.0$,
strongly interacting bosons attain the fermionic density distribution.

##  III Results of expansion dynamics

###  III.1 Measures on one-body level

We choose $N=5$ strongly interacting bosons with $\lambda=25$, in the harmonic
oscillator trap $V(x)=\frac{1}{2}x^{2}$. We use MCTDHB for bosons and for the
noninteracting fermions we use MCTDHF (multiconfigurational time-dependent
Hartree for fermions) implemented in the MCTDH-X software package [33]. The
computation is done with $M=24$ orbitals for strongly interacting bosons, the
convergence is assured as the last orbital occupation is insignificant. For
noninteracting fermions, computation is done with $M=5$ orbitals, each orbital
has exactly equal contribution of $20\%$ population. The initial state is
prepared at the ground state in the harmonic oscillator and the 1D expansion
is studied on sudden removal of the trap. We monitor the dynamical evolution
in one and higher-order many-body densities.
Fig.1 depicts a comparison of dynamical evolution of reduced one-body density
in $x$\- space $\rho^{(1)}(x,x^{\prime})$ between the strongly interacting
bosons (left panel) and non interacting fermions (right panel), exhibiting
identical expansion dynamics. Initially at $(t=0)$, the density is clustered
at the center due to the harmonic trap, but becomes flatter and broader with
increase in time. The density gradually acquires modulations and the number of
humps tries to saturate to the number of bosons. At $t=8$, we observe
emergence of five distinct humps. The hump at the centre is the brightest
where the potential was initially zero, the two outer humps are less
pronounced due to the initial non-zero confining potential and the two
outermost humps are least pronounced which are away from the center. It is
clearly seen that both strongly interacting bosons and non-interacting
fermions exhibit identical expansion dynamics in $x$-space. We additionally
note that the density’s maxima in the TG regime are distinct but not isolated,
the five bright spots are interconnected.

Figure 3: Dynamics of relative entropy measures $K$, $SK$ and $J$ in momentum
space as determined by the equations (18), (19) and (20) respectively. The
corresponding one-body momentum density is utilized to evaluate the entropy
measures. All the measures exhibit how the strongly interacting bosons
asymptotically approach fermionization limit at time $t=8.0$ when all the
relative entropy measures reach close to zero value.

Next we study the expansion dynamics through reduced one-body density in
$k$-space, i.e $\rho^{(1)}(k,k^{\prime})$ which evaluates the time for
dynamical fermionization. Fig.2 describes the momentum distribution of the
expanding gas with five bosons at different time ($t$= 0,1,2,3,5,6,8).
Initially at $t=0$, the density is in cluster form at the center of the trap,
with time it gradually develops the peaks, the peaks become prominent and
finally become identical with the fermionic density distribution (last panel
on the right side) at $t=8.0$. We observe dynamical fermionization occurs as
the momentum distribution approaches ideal Fermi distribution at time $t=8.0$.
It is not possible to distinguish the fermionized bosons and the ideal
fermions from the measure of one-body $k$-density. We conclude $t=8$ is the
time for dynamical fermionization for the present system.

We calculate the measures of relative entropies $K,SK$ and $J$ in the
$x$-space using Eqs.(15),(16) and (17) respectively. However as shown in
Fig.1, as the fermionized bosons and noninteracting fermions exhibit identical
one-body expansion dynamics in real space, the corresponding relative entropy
measures are identically zero all throughout the dynamical evolution. The
corresponding measures in momentum space as determined by the Eqs. (18), (19)
and (20) are plotted in Fig.3. All the measures have the same trend and
conclude the same physics. $K$ and $SK$ almost overlap whereas $J$ lies below.
We find, initially the values are significant which infer that the two
systems, strongly interacting bosons and fermions are distinctly different.
With time, all the three measures gradually decrease which signify that the
strongly interacting bosons gradually attain the fermionic characteristics.
Finally at the time of fermionization ($t=8.0$), they asymptotically reach
close to zero which signifies that the density of the two systems become
identical.

All the measures based on the dynamics of one-body density support the well
accepted fact that at the time of fermionization the strongly correlated
bosons asymptotically achieve fermionization. The one-body momentum density of
the two systems become indistinguishable. However it needs to examine whether
the same conclusion can be made from the measures using two-body densities.

###  III.2  Measures on two-body level

We discuss dynamics of the reduced two-body density in $x$-space
$\rho^{(2)}(x,x^{\prime})$ in Fig.4. Comparison between strongly interacting
bosons and non-interacting fermions are presented at the same time scale as of
one-body density [Fig.1]. Left panel corresponds to the strongly interacting
bosons and the right panel corresponds to the non-interacting fermions.
Initially, the atoms are clustered near the center ($x=x^{\prime}=0$) for both
the cases of interacting bosons and non-interacting fermions. With time,
$\rho^{(2)}$ gradually spreads out to the off-diagonal $(x\neq x^{\prime})$.
The diagonal $(x\simeq x^{\prime})$ is gradually depleted and a so-called
”correlation hole” forms on the diagonal. Thus the bosons gradually overcome
their spatial overlap; the probability of finding two bosons at the same
position tends towards zero. With time, the ’correlation hole’ spreads and
finally converges in the fermionization limit at time $t=8.0$. Off-diagonal
spreading also converge in the fermionization time. Similar to the one-body
density, the maxima which are formed across the off-diagonal are distinct but
not isolated. We infer that the correlation hole along the diagonal and the
confined spreading at the time of fermionization are the unique signatures of
the two-body density of a dynamically evolved fermionized state. The dynamics
exhibited by the non-interacting fermions is exactly identical to that of
interacting bosons.

Figure 4: Snapshots of the time evolution of reduced two-body density in
$x$-space ($\rho^{(2)}(x,x^{\prime})$) during the expansion of $N=5$ strongly
interacting bosons with interaction strength $\lambda=25$ are presented (left
panel). Computation is done with $M=24$ orbitals. Comparison is made with
dynamics of expansion for non- interacting fermions with $M=5$ orbitals (right
panel). Identical evolution of the two systems is observed.  Figure 5:
Snapshots of the time evolution of the reduced two-body momentum density
($\rho^{(2)}(k,k^{\prime})$) during the expansion of $N=5$ strongly
interacting bosons with interaction strength $\lambda=25$ and orbitals $M=25$
are presented for different times. The corresponding fermionic distribution
computed with $M=5$ is presented at the last panel on right side. The strongly
interacting bosons which was initially at the center of the trap, gradually
develops very weak correlation hole along the diagonal. However at the time of
dynamical fermionization, the correlation width and intensity strongly differ
from the distinct correlation hole in the fermionic distribution.

In Fig. 5, we plot the reduced two-body density in $k$-space
$\rho^{(2)}(k,k^{\prime})$ for the same time points as chosen for Fig.2.
Initially at time $t=0$, we observe the density in the clustered state at the
center. With time $\rho^{(2)}$ spreads out to the off-diagonal and the
diagonal is depleted. However no distinct correlation hole is created, the
strongly interacting bosons try to make a diagonal gap, but some internal
structure along the diagonal exist. $\rho^{(2)}(k,k^{\prime})$ remains
localized near $k\leq 3.5$, because of the finite energy. The well defined
bright spot indicates the localization of the bosons in momentum space. With
further time evolution the internal structure along the diagonal is greatly
reduced but not completely extinguished like ideal fermionic two-body reduced
momentum distribution.

To understand qualitatively how close the strongly interacting bosons feature
the fermionic properties in the two-body level, we utilize the information
distance measures. The measures utilizing two-body densities in $x$-space from
Eqs. (21), (22) and (23) result to zero as the fermionized bosons and
noninteracting fermions exhibit indistinguishable expansion dynamics as shown
in Fig.4. The corresponding measures utilizing two-body momentum density from
Eqs (24), (25) and (26) are plotted in Fig.6. The initial large value in all
the relative entropy measures clearly distinguish that strongly interacting
boson’s density is completely different from that of fermionic two-body
density. With increase in time, the entropy measures fall which signifies that
the two-body density of bosons try to achieve fermionic two-body density. At
the time of fermionization ($t=8.0$), all the measures settle to a significant
minimum value but far above zero leading to infer that dynamically fermionized
bosons have distinguishable density from that of non-interacting fermions in
the two-body level. It signifies that although at the time of fermionization,
one-body density of strongly interacting bosons becomes identical to that of
ideal fermions, but their corresponding two-body density differs. We conclude
that the bosons in TG limit is close to fermions in ’two-body perspective
analysis’, but not identical as concluded from ’one-body perspective’
analysis.

Figure 6: Dynamics of relative entropy measures between the strongly
interacting bosons and the ideal fermions using two-body momentum density in
the equations (24) to (26). The initial large value clearly signify how the
two systems are quite far from each other. With time all the three measures
gradually decreases with time as the strongly interacting bosons try to
acquire fermionic distribution. However unlike the observation made in Fig. 3,
$K$, $SK$ and $J$ settle to a significant minimum value. Figure 7: The time
evolution of the local two-body correlation $g^{(2)}(0,0)$ for the strongly
interacting bosons in momentum space. $g^{(2)}(0,0)$ = $0$ corresponds to
fermionization limit which infers existence of distinct correlation hole; two-
body correlation is completely extinguished. Strongly interacting bosons
initially have strong local correlation, it dies with time but does not reach
to zero at the time of fermionization. Figure 8: The time evolution of two-
body non-local correlation $g^{(2)}(0,k)$ for the strongly interacting bosons.
$g^{(2)}(0,k)$ $\rightarrow$ $1$ implies the ’no correlation’ limit which is
achieved at the characteristic time of fermionization. We also observe some
peculiar oscillation around the ’no correlation’ zone for the strongly
interacting bosons.

To estimate further the difference between the strongly interacting bosons and
the fermions we calculate the local two body correlation function both in $x$
and $k$-space. From Fig. 4, it is clearly seen that the two-body correlation
is completely extinguished at $x=x^{\prime}=0$ all throughout the dynamics
both for the fermionized bosons and non-interacting fermions, it leads to
$g^{(2)}(x=0,x^{\prime}=0)$ $=0$ as calculated from the Eq. (13). However in
the $k$-space, the local two-body correlation $g^{(2)}(k=0,k^{\prime}=0)$
would contain significant quantitative information about the close proximity
of the two systems in the entire dynamics. From the Eq. (14) we calculate
$g^{(2)}(k=0,k^{\prime}=0)$ for fermions and fermionized bosons and plot it in
Fig. 7. It is expected that for fermions $g^{(2)}(0,0)$ is uniquely zero which
is the typical feature of fermionization. For strongly interacting bosons it
starts from $g^{(2)}(0,0)=1.0$, and then after a small increase, it smoothly
decreases. However at the time of fermionization, it settles to $g^{(2)}=0.4$.
It indicates that for strongly interacting bosons the local two-body
correlation is not completely extinguished, which helps to identify strongly
interacting bosons from non-interacting fermions. The corresponding two-body
non-local correlation $g^{(2)}(0,k)$ is plotted in Fig. 8. Strongly
interacting bosons reach to the ’no correlation’ limit with some initial
characteristic oscillation. However the fermions attains the ’no correlation’
limit sharply and does not exhibit any resolvable dynamical structure.

###  III.3 Higher-order coherence

Fig.9 and Fig. 10 depict the higher-order densities $\rho^{(p)}(p=3,4)$ in $x$
and $k$-space for the strongly interacting bosons; the non-interacting
fermions do not exhibit any higher-order correlation beyond two-body. The
$\rho^{(p)}$ provides the probability of detecting particles at $x$ and
$x^{\prime}$ when the remaining $p-2$ particles are fixed at some reference
positions. For fermionized bosons, the diagonal of the high- order $p$-body
densities $\rho^{(p)}(x,x)$ vanishes as the so-called correlation hole results
from the strong interaction strength that mimics the Pauli principle,
preventing to find two bosons at the same position. The maxima exhibit well
defined peaked structure which indicates the localization of the atoms in
position space. The maxima along the anti-diagonal $x=-x^{\prime}$ infers that
the bosons maximize the distance between each other. Whereas the maxima along
the sub-diagonal infers that the bosons minimize the potential energy. The
additional correlation holes appear at the fixed values of the remaining $p$
$-$ $2$ coordinates of $\rho^{(p)}$, preventing to find other bosons at these
positions.

In Fig.9, we plot three-body density $\rho^{(3)}(x,x^{\prime})$ (keeping the
third particle at the fixed reference point $x_{3}=0$) at the left column and
four-body density $\rho^{(4)}(x,x^{\prime})$ (keeping the third and the fourth
particle at the fixed reference points $x_{3}=0$ and $x_{4}=0.47$) at the
right column. Initially both three-body and four-body density are localized at
the center. With time, the maxima is developed along the anti-diagonal as well
as along the sub-diagonal. Additional correlation holes also appear at the
fixed reference points as probability of finding other bosons at these
positions is zero. The three-body correlation structure exhibit further
delocalization as expected and finally converge at the time of dynamical
fermionization. The corresponding three-body and four-body densities in the
momentum space are plotted in Fig. 10. Like Fig.9, it also exhibits rich
dynamical structure.

Figure 9: Snapshots of the dynamics of three-body density
$\rho^{(3)}(x,x^{\prime})$ (left column) and $\rho^{(4)}(x,x^{\prime})$ (right
column) for $N=5$ strongly interacting bosons. See the text for details.
Figure 10: Snapshots of the dynamics of three-body density
$\rho^{(3)}(k,k^{\prime})$ (left column) and $\rho^{(4)}(k,k^{\prime})$ for
$N=5$ strongly interacting bosons. See the text for details.

##  IV Conclusion

For the strongly interacting TG limit, the wave function of the interacting
bosons in one dimension maps to that of non interacting fermions. The
properties of the two systems which are determined by the amplitude of the
wave function are identical. Whereas the measurement in the momentum space
exhibits the dynamical fermionization. This conclusion is established on the
basis of one-body analysis. In this work we extend the measure to higher order
densities and correlation function using numerically exact solution of the
time dependent many-boson Schrödinger equation. We establish all the existing
physics and additionally include the measures of three different kinds of
relative entropies which utilize the one-body in the momentum space. We
clearly find that all the measures on one-body level unambiguously establish
the onset of fermionization in the dynamical expansion of strongly interacting
bosons. However the measures beyond one-body level is able to distinguish the
strongly interacting bosons from noninteracting fermions in the dynamics. To
understand that how the two systems reach to close proximity at the onset of
fermionization we extend all the measures in two-body level. The two-body
momentum distribution clearly distinguish the two systems in the dynamical
evolution. The measures of relative entropies using two-body momentum density
also quantitatively exhibit the distinction between the two systems. The two-
body local and non local correlations also conclude the same difference.
Additionally we find that the interacting bosons exhibit very rich structure
in the higher body densities whereas for the fermions three- and higher body
correlation are ideally zero.

Acknowledgements – This work was supported by the FAPESP grant Process No.
2023/06550-4. We are thankful to Rhombik Roy for helpful discussion.

## References

  * Lieb and Liniger [1963] E. H. Lieb and W. Liniger, Exact analysis of an interacting bose gas. i. the general solution and the ground state, Phys. Rev. 130, 1605 (1963).
  * Lieb [1963] E. H. Lieb, Exact analysis of an interacting bose gas. ii. the excitation spectrum, Phys. Rev. 130, 1616 (1963).
  * Gangardt and Shlyapnikov [2003a] D. M. Gangardt and G. V. Shlyapnikov, Stability and phase coherence of trapped 1d bose gases, Phys. Rev. Lett. 90, 010401 (2003a).
  * Menotti and Stringari [2002] C. Menotti and S. Stringari, Collective oscillations of a one-dimensional trapped bose-einstein gas, Phys. Rev. A 66, 043610 (2002).
  * Petrov  _et al._ [2000] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, Regimes of quantum degeneracy in trapped 1d gases, Phys. Rev. Lett. 85, 3745 (2000).
  * Girardeau [1965] M. D. Girardeau, Permutation symmetry of many-particle wave functions, Phys. Rev. 139, B500 (1965).
  * Gangardt and Shlyapnikov [2003b] D. M. Gangardt and G. V. Shlyapnikov, Local correlations in a strongly interacting one-dimensional bose gas, New J. Phys 5, 79 (2003b).
  * Minguzzi and Gangardt [2005] A. Minguzzi and D. M. Gangardt, Exact coherent states of a harmonically confined tonks-girardeau gas, Phys. Rev. Lett. 94, 240404 (2005).
  * Rigol and Muramatsu [2005] M. Rigol and A. Muramatsu, Fermionization in an expanding 1d gas of hard-core bosons, Phys. Rev. Lett. 94, 240404 (2005).
  * Muth  _et al._ [2010] D. Muth, B. Schmidt, and M. Fleischhauer, Fermionization dynamics of a strongly interacting one-dimensional bose gas after an interaction quench, New J. Phys. 12, 083065 (2010).
  * Yukalov and Girardeau [2005] V. I. Yukalov and M. D. Girardeau, Fermi-bose mapping for one-dimensional bose gases, Laser Phys. Lett. 2, 375 (2005).
  * Wilson  _et al._ [2020] J. M. Wilson, N. Malvania, Y. Le, Y. Zhang, M. Rigol, and D. S. Weiss, Observation of dynamical fermionization, Science 367, 1461 (2020).
  * Streltsov  _et al._ [2006] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, General variational many-body theory with complete self-consistency for trapped bosonic systems, Phys. Rev. A 73, 063626 (2006).
  * Streltsov  _et al._ [2007] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Role of excited states in the splitting of a trapped interacting bose-einstein condensate by a time-dependent barrier, Phys. Rev. Lett. 99, 030402 (2007).
  * Alon  _et al._ [2007] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Unified view on multiconfigurational time propagation for systems consisting of identical particles, J. Chem. Phys. 127, 154103 (2007).
  * Alon  _et al._ [2008] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Multiconfigurational time-dependent hartree method for bosons: Many-body dynamics of bosonic systems, Phys. Rev. A 77, 033613 (2008).
  * Lode [2016] A. U. J. Lode, Multiconfigurational time-dependent hartree method for bosons with internal degrees of freedom: Theory and composite fragmentation of multicomponent bose-einstein condensates, Phys. Rev. A 93, 063601 (2016).
  * Fasshauer and Lode [2016] E. Fasshauer and A. U. J. Lode, Multiconfigurational time-dependent hartree method for fermions: Implementation, exactness, and few-fermion tunneling to open space, Phys. Rev. A 93, 033635 (2016).
  * Lode  _et al._ [2020] A. U. J. Lode, C. Lévêque, L. B. Madsen, A. I. Streltsov, and O. E. Alon, Colloquium: Multiconfigurational time-dependent hartree approaches for indistinguishable particles, Rev. Mod. Phys 92, 011001 (2020).
  * Kvaal [2013] S. Kvaal, Variational formulations of the coupled-cluster method in quantum chemistry, Molecular Physics 111, 1100 (2013).
  * Kramer and Saraceno [1981] P. Kramer and M. Saraceno,  _Geometry of the Time-Dependent Variational Principle in Quantum Mechanics_ (Springer Berlin, Heidelberg, 1981).
  * McLachlan [1964] A. McLachlan, A variational solution of the time-dependent schrodinger equation, Molecular Physics 8, 39 (1964).
  * Cao  _et al._ [2017] L. Cao, V. Bolsinger, S. I. Mistakidis, G. M. Koutentakis, S. Krönke, J. M. Schurer, and P. Schmelcher, A unified ab initio approach to the correlated quantum dynamics of ultracold fermionic and bosonic mixtures, J. Chem. Phys. 147, 044106 (2017).
  * Lin  _et al._ [2021] R. Lin, C. Georges, J. Klinder, P. Molignini, M. Büttner, A. U. J. Lode, R. Chitra, A. Hemmerich, and H. Keßler, Mott transition in a cavity-boson system: A quantitative comparison between theory and experiment, SciPost Phys. 11, 30 (2021).
  * Nguyen  _et al._ [2019] J. H. V. Nguyen, M. C. Tsatsos, D. Luo, A. U. J. Lode, G. D. Telles, V. S. Bagnato, and R. G. Hulet, Parametric excitation of a bose-einstein condensate: From faraday waves to granulation, Phys. Rev. X 9, 011052 (2019).
  * Lode  _et al._ [2012] A. U. J. Lode, K. Sakmann, O. E. Alon, L. S. Cederbaum, and A. I. Streltsov, Numerically exact quantum dynamics of bosons with time-dependent interactions of harmonic type, Phys. Rev. A 86, 063606 (2012).
  * Lode  _et al._ [2015] A. U. J. Lode, B. Chakrabarti, and V. K. B. Kota, Many-body entropies, correlations, and emergence of statistical relaxation in interaction quench dynamics of ultracold bosons, Phys. Rev. A 92, 033622 (2015).
  * Bera  _et al._ [2019] S. Bera, R. Roy, A. Gammal, B. Chakrabarti, and B. Chatterjee, Probing relaxation dynamics of a few strongly correlated bosons in a 1d triple well optical lattice, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 215303 (2019).
  * Roy  _et al._ [2018] R. Roy, A. Gammal, M. C. Tsatsos, B. Chatterjee, B. Chakrabarti, and A. U. J. Lode, Phases, many-body entropy measrlres, and coherence of interacting bosons in optical lattices, Phys. Rev. A 97, 043625 (2018).
  * Roy  _et al._ [2019] R. Roy, C. Lévêque, A. U. J. Lode, A. Gammal, and B. Chakrabarti, Fidelity and entropy production in quench dynamics of interacting bosons in an optical lattice, Quantum Reports 1, 304 (2019).
  * Roy  _et al._ [2022] R. Roy, B. Chakrabarti, and A. Trombettoni, Quantum dynamics of few dipolar bosons in a double-well potential., European Physical Journal D 76, 215303 (2022).
  * Lin  _et al._ [2020] R. Lin, P. Molignini, L. Papariello, M. C. Tsatsos, C. Lévêque, S. E. Weiner, E. Fasshauer, and R. Chitra, Mctdh-x: The multiconfigurational time-dependent hartree method for indistinguishable particles software, Quantum Sci. Technol. 5, 024004 (2020).
  * Lode  _et al._ [2024] A. U. J. Lode, M. C. Tsatsos, E. Fasshauer, S. E. Weiner, R. Lin, L. Papariello, P. Molignini, C. Lévêque, M. Büttner, J. Xiang, S. Dutta, and Y. Bilinskaya, Mctdh-x: The multiconfigurational time-dependent hartree method for indistinguishable particles software (2024).
  * Sakmann  _et al._ [2008] K. Sakmann, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Reduced density matrices and coherence of trapped interacting bosons, Phys. Rev. A 78, 023615 (2008).
  * Salazar  _et al._ [2020] S. J. C. Salazar, H. G. Laguna, V. Prasad, and R. P. Sagar, Shannon-information entropy sum in the confined hydrogenic atom, Int J Quantum Chem. 120, e26188 (2020).
  * Aquino  _et al._ [2013] N. Aquino, A. Flores-Riveros, and J. Rivas-Silva, Shannon and fisher entropies for a hydrogen atom under soft spherical confinement, Phys. Lett. A 377, 2062 (2013).
  * Jiao  _et al._ [2017] L. G. Jiao, L. R. Zan, Y. Z. Zhang, and Y. K. Ho, Benchmark values of shannon entropy for spherically confined hydrogen atom, Int J Quantum Chem. 117, e25375. (2017).
  * Sen [2005] K. D. Sen, Characteristic features of shannon information entropy of confined atoms, J. Chem. Phys. 123, 074110 (2005).
  * Sriraman  _et al._ [2017] T. Sriraman, B. Chakrabarti, A. Trombettoni, and P. Muruganandam, Characteristic features of the shannon information entropy of dipolar bose-einstein condensates, J. Chem. Phys. 147, 044304 (2017).
  * Kullback [1959] S. Kullback,  _Statistics and Information Theory_ (Wiley, New York, 1959).
  * Rao [1987] C. Rao,  _Differential Geometry in Statistical Interference, IMS-Lecture Notes_ , Vol. 10 (Wiley, New York, 1987) p. 217.
  * Baratpour and Habibi [2012] S. Baratpour and R. A. Habibi, Testing goodness-of-fit for exponential distribution based on cumulative residual entropy, Commun. Stat. - Theory Methods 41, 1387 (2012).
  * Cover and Thomas [2006] T. M. Cover and J. A. Thomas,  _Elements of Information Theory_ (John Wiley and Sons, Hoboken, NJ, 2006).
  * J.Lin [1991] J.Lin, Divergence measures based on the shannon entropy, IEEE TRANSACTIONS ON INFORMATION THEORY. 37, 145 (1991).
  * Majtey  _et al._ [2004] A. Majtey, P. Lamberti, M. Martin, and A. Plastino, Wootters’ distance revisited: a new distinguishability criterium, arxiv (2004).


# Theoretically Modelling Photoionized Regions with Fractal Geometry in Three
Dimensions

Yifei Jin (金刈非) Research School for Astronomy & Astrophysics, Australian
National University, Canberra, Australia, 2611  ARC Centre of Excellence for
All Sky Astrophysics in 3 Dimensions (ASTRO 3D)  Lisa J. Kewley Research
School for Astronomy & Astrophysics, Australian National University, Canberra,
Australia, 2611  ARC Centre of Excellence for All Sky Astrophysics in 3
Dimensions (ASTRO 3D)  Ralph S. Sutherland Research School for Astronomy &
Astrophysics, Australian National University, Canberra, Australia, 2611

(Received xx, xx, 2022; Revised xx, xx, 2022; Accepted xx, xx, 2022)

###### Abstract

We create a photoionization model embedded in the turbulent ISM by using the
state-of-the-art Messenger Monte-Carlo MAPPINGS V code (M3) in conjunction
with the CMFGEN stellar atmosphere model. We show that the turbulent ISM
causes the inhomogeneity of electron temperature and density within the
nebula. The fluctuation in the turbulent ISM creates complex ionization
structures seen in nearby nebulae. The inhomogeneous density distribution
within the nebula creates a significant scatter on the spatially-resolved
standard optical diagnostic diagrams, which cannot be represented by the
spherical constant density photoionization model. We analyze the dependence of
different optical emission lines on the complexity of nebular geometry,
finding that the emission-lines residing on the nebular boundary are highly
sensitive to the complexity of nebular geometry, while the emission-lines
produced throughout the nebula are sensitive to the density distribution of
the ISM within the nebula. Our fractal photoionization model demonstrates that
a complex nebular geometry is required for accurate modeling of H ii regions
and emission-line galaxies, especially for the high-redshift galaxies, where
the ISM is highly turbulent based on the increasing observational evidence.

galaxies, abundances — galaxies: fundamental parameters — galaxies: high-
redshift

††journal: ApJL

##  1 Introduction

Robust photoionization models are fundamental to understand galaxy evolution
and chemical evolution. Theoretical emission-line models are required to
measure the chemical abundance (McGaugh, 1991; Kewley & Dopita, 2002; Dopita
et al., 2006; Kewley et al., 2019), the electron density (Osterbrock, 1989;
Kewley et al., 2019) and the pressure of interstellar medium (ISM) (Kewley et
al., 2019). With these diagnostic tools, we can analyze the chemical evolution
of galaxies (Kobulnicky & Kewley, 2004; Tremonti et al., 2004; Zahid et al.,
2011; Yuan et al., 2013; Strom et al., 2022), the star formation history of
galaxies (Hopkins & Beacom, 2006; Madau & Dickinson, 2014; Tacchella et al.,
2022), and the power sources of galaxies across cosmic time (Kewley et al.,
2013b; Rigby et al., 2021).

One of the most powerful applications of theoretical photoionization models
has been the determination of power sources of galaxies using strong emission-
lines available in the optical spectra of galaxies. Baldwin et al. (1981)
first proposed usage of the [N ii]/H$\alpha$ and [O iii]/H$\beta$ ratios to
classify the excitation sources in galaxies, which is known as the
Baldwin–Phillips–Terlevich (BPT) diagram. The [N ii]/H$\alpha$ and [O
iii]/H$\beta$ ratios excited by active galactic nuclei (AGN), shocks and
planetary nebulae are larger than the ratios excited by star formation due to
the difference in the hardness of ionizing radiation fields. Veilleux &
Osterbrock (1987) derived semi-empirical classification scheme and extend the
standard optical diagnostic line ratios to include the [S ii]/H$\alpha$ and [O
i]/H$\alpha$ ratios.

Kewley et al. (2001) derived the first theoretical classification scheme for
local galaxies by computing photoionization and shock models with the MAPPINGS
III photoionization code (Binette et al., 1985; Sutherland & Dopita, 1993) in
conjunction with the stellar population synthesis models (Leitherer et al.,
1999). Kewley et al. (2013a, b) further model the changes of the diagnostic
line ratios as a function of cosmic time, finding the [N ii]/H$\alpha$ and [O
iii]/H$\beta$ ratios are likely to be large in high-redshift galaxies due to
the extreme ISM conditions in high-redshift (Steidel et al., 2014; Shapley et
al., 2005).

All previous models for the spectra in star-forming galaxies assume a
spherical or plane parallel geometry. However, spatially-resolved spectra of
local H ii regions show the complex structures within the nebula. Kinematic
and morphology studies reveal the filaments (Gendron-Marsolais et al., 2018;
Zavala et al., 2022) and turbulent structures (Rubin et al., 2011) within the
nebula

These complex structures are not spherical or plane-parallel (Jin et al. in
prep). Snijders et al. (2007) show that a clumpy photoionization model is
required to reproduce the extremely high ionized densities and ionization
parameters as found in the local starburst galaxies M82 and the Antennae.

In M82 and the Antennae, the young massive H ii regions are studied as
candidates of the star forming regions in high-redshift galaxies through their
dense and compact morphology (Smith et al., 2006) and the large ionization
parameters (Snijders et al., 2007; Rigby et al., 2011). In particular, the
large ionization parameters seen in local extreme star forming regions appear
to be widespread in the current set of high-redshift galaxies observed using
near infrared spectroscopy (Shirazi et al., 2014; Sanders et al., 2016;
Papovich et al., 2022).

In this Letter, we present the first geometrical photoionization model using
over 30 chemical elements embedded in the turbulent ISM by using $\rm M^{3}$
(Messenger Monte-Carlo MAPPINGS V Jin et al., 2022), with detailed stellar
atmosphere models as the input ionizing source. We present the complex
structures of the ionizing radiation field, the electron temperature, the
electron density and the ionization states of different elements within the
model. We investigate how the nebular geometry affects the standard optical
diagnostic emission-lines in the scope of both spatially-resolved and
integrated spectroscopy. Our results have important implications for studies
of emission-line ratios in high-redshift galaxies.

##  2 The method

Figure 1: Upper panel: The projected H$\beta$ luminosity map of 3D nebula with
fractal geometry. Lower panel: The projected distribution of log-normally
distributed ISM. The white curve outlines the shape of our modelled nebula.
The ionizing source is at the centre of the domain (indicated by the white
star).

$\rm M^{3}$ (Jin et al., 2022) is a descendant of the MAPPINGS V
photoionization code (Sutherland & Dopita, 1993). $\rm M^{3}$ combines the
Monte-Carlo radiative transfer technique (Lucy, 1999) with MAPPINGS V (Binette
et al., 1985), which comprehensively models the microphysics of over 30
chemical elements in the ISM. $\rm M^{3}$ is designed to create
photoionization models with arbitrary geometry in three dimensions, by
considering both stellar and diffuse ionizing radiation fields.

Our model is set up in a fractal ISM density field (thereafter $the$ $fractal$
$model$) to mimic the hierarchical structures of turbulence in ISM (Larson,
1981). Following Sutherland & Bicknell (2007), we generate the fractal ISM
density field by adopting a standard log-normal density distribution with
$\mu=1.0$ and $\sigma^{2}=5.0$ and a Kolmogorov spatial structure power law.
The average density of the ISM density field is $n_{H}$=100 $\rm cm^{-3}$. The
static log-normal density distribution is mathematically well-constrained and
shows the structure and density fluctuations similar to real turbulent
molecular clouds. The fractal ISM density field reproduces the clumps and
turbulent structures seen in the real nearby HII region, like Orion (Rubin et
al., 2011). We adopt the solar abundance based on the results from Asplund et
al. (2009). In this Letter, we are creating a dust-free photoionization model
without dust depletion.

We select the ionizing spectrum from the CMFGEN stellar library (Hillier &
Miller, 1998, 1999). The CMFGEN library accurately describes the UV stellar
spectra by using non-LTE stellar atmosphere. In this work, we choose an O-star
with the luminosity $L$=$10^{6}$ $L_{\sun}$, the temperature $\rm
T_{eff}$=40000 K and the gravity lg($g$)=4. The selected spectrum well
describes the hard ionizing radiation from massive stars (Simón-Díaz &
Stasińska, 2008).

Our model is set up in a cube with 129$\times$129$\times$129 spaxels (spatial
pixels). The physical size of the entire computational domain is
43$\times$43$\times$43 pc3. The size is selected based on the size-density
relation for the H ii regions in both star-forming galaxies and starburst
galaxies (Hunt & Hirashita, 2009). A spatial resolution of 0.25 pc per spaxel
guarantees that the sub-structures of the nebulae can be resolved.

We also produce an equivalent spherical nebula model as a reference to
illustrate the geometric effects on emission-line distributions, emission-line
fluxes and ratios. The spherical model is created by MAPPINGS V with an
assumption of spherical nebular geometry, in conjunction with the stellar
ionizing spectrum from the CMFGEN stellar library. Except for the geometry,
the input parameters of the spherical model are the same as the input
parameters of the fractal model.

##  3 Results

###  3.1 The Radiation Bounded Nebula

Figure 1 presents the 2D map of the integrated H$\beta$ emission-line
luminosity of the fractal model. The nebula has a fractal boundary consisting
of pixels ionized by both stellar and diffuse photons. The nebula shape is
determined by the surrounding ISM clumps. The photoionized region extends
furthest in low-density regions around the ionizing source. The low density of
the ISM in these regions allows photons to easily pass through. In the dense
regions, the ionizing photons are absorbed by the dense clumps, leading to a
nebular boundary much closer to the central ionizing source. Beyond the dense
nebular boundary, there are regions ionized by diffuse photons, which are
scattered by dense clumps. These regions which are ionized by the diffuse
photons rather than the stellar photons are defined as the diffuse ionized
regions.

Top left panel in Figure 2 shows the local ionization parameter distribution
across the middle slice of the nebula. The local ionization parameter is
defined as:

$U=\frac{f_{\nu>13.6eV}}{cN_{H}}$ (1)

where $f_{\nu>13.6eV}$ is the flux of the ionizing photons above 13.6 eV
through a unit area, $N_{H}$ is the local number density of hydrogen and $c$
is the speed of light. The ionization parameter indicates the capability of
ionizing radiation field to ionize neutral gas. The average dimensionless
ionization parameter is log$U$=-2.8, within the range of the typical
ionization parameter -3.2$<$log$U$$<$-2.7 for local Hii regions (Dopita et
al., 2000) and star-forming galaxies (Moustakas et al., 2010). The diffuse
ionized regions have an ionization parameter as low as log$U$=-5 because these
regions are ionized by weak diffuse ionizing photons.

The electron temperature (top central panel in Figure 2) ranges from 5000 K to
10000 K with average inhomogeneity of $\langle T_{e}\rangle$=700 K within the
main body of the nebula. The electron temperature increases with radius on
average but shows strong azimuthal asymmetry. The positive temperature
gradient is caused by the lack of coolant ions when the ionization parameter
declines at large nebular radius. The temperature is 5000 K at the nebular
center and increases to 10000 K at the edge of the nebula. The diffuse ionized
regions have a low electron temperature of 2000 K because the weak ionization
field can only weakly heat these regions. The electron temperature shows a
spatial variation within the nebula. The high-density clumps have higher
electron temperature and the low-density regions have lower electron
temperature.

The electron density distribution (shown in Figure 2 top right panel) presents
a similar log-normal distribution as the neutral gas density. The electron
density spans four orders of magnitude from 0.1 $\rm cm^{-3}$ to 800 $\rm
cm^{-3}$. The average electron density is $n_{e}$=30 $\rm cm^{-3}$. The
inhomogeneity of electron density, which is described by the standard
deviation, is $\langle n_{e}\rangle$=35 $\rm cm^{-3}$.

###  3.2 The Ionization Structure

Figure 2: Maps of fundamental parameters across the Z=0 pc slice of the
nebula. Grey contours indicate the distribution of the input of ISM density.
We present maps of the dimensionless ionization parameter $U$, the electron
temperature $T_{e}$, the electron density $n_{e}$ and the distribution of
fluxes of the H$\alpha$, H$\beta$, [O iii], [O i], [N ii] and [S ii] emission
lines. The nebula is assumed at the distance of the Large Magellanic Cloud of
48.5kpc.

Figure 2 (lower two panels) presents the internal structure of the six
diagnostic emission-lines, H$\alpha$, H$\beta$, [O iii]$\lambda$5007, [O
i]$\lambda$6300, [N ii]$\lambda$6584 and [S ii]$\lambda\lambda$6717,31. These
emission-lines show diverse structures based on their critical densities and
ionization potentials.

The Balmer H$\alpha$ and H$\beta$ lines are produced uniformly throughout the
main body of the nebula, extending to the diffuse ionized region. The
variation of the H$\alpha$ and H$\beta$ luminosity is smaller than one
magnitude across the entire nebula. The Balmer lines in the diffuse ionized
region are fainter than the Balmer lines from the denser ionized gas. On
average, the H$\alpha$ and H$\beta$ luminosity in the diffuse ionized region
is 0.8 magnitude fainter than the average luminosity of the dense ionized
region.

The [O iii] emission-line is also produced throughout the main body of nebula,
which is ionized by stellar photons. In the diffuse ionized region, the [O
iii] emission is fainter than the [O iii] emission from the nebular main body
by nine magnitudes. The extremely faint diffuse [O iii] emission is caused by
the fact that the hardness of diffuse ionization field is too low to produce
$\rm O^{++}$ ions, which has the ionization potential of 35.1 eV.

The [O i], [N ii] and [S ii] lines are produced at the outer edge of the
nebula. As shown in Figure 3, 99.5% of the [O i] emission, 87.2% of the [N ii]
emission and 87.5% of the [S ii] emission are from the outer 10% radius of
nebula. In the turbulent Hii region model, the [O i], [N ii] and [S ii] lines
highlight the “coastline” of the nebula. The fractal nebular boundary extends
the perimeter of the outer edge of the nebula, increasing the total luminosity
of [O i], [N ii] and [S ii] emission-lines compared with a standard spherical
model.

The inhomogeneous density distribution of the ISM leads to the substructures
of [O i], [N ii] and [S ii] lines in the fractal model. The emission-line
fluxes present a clumpy distribution within the nebula, where bright clumps
trace high-density regions in ISM. We also find an ionization bar at the
center of nebula crossing from the upper-left to the lower-right of the z=0
slice. The substructures seen in the fractal model are similar to the
ionization structures seen in the Orion Nebula, where the complex ionization
structures coexist with the complex density structures (Rubin et al., 2011).

###  3.3 The Spatially-Resolved Standard Optical Diagnostic Diagrams

Figure 3: Spherical MAPPINGS V model. The spherical model has the same inputs
to the fractal model except for the geometry. Here, we present the radial
distribution of the [O iii], [O i], [N ii] and [S ii] lines.  Figure 4: The
comparison of diagnostic emission-line ratios from the fractal model against
the spherical photoionization model. We employ the [N ii]/H$\alpha$ versus [O
iii]/H$\beta$ (left column), [S ii]/H$\alpha$ versus [O iii]/H$\beta$ (middle
column) and [O i]/H$\alpha$ versus [O iii]/H$\beta$ (right column) diagrams.
The background grey scale indicates the number density of spaxels from the
collapsed 2D map of the fractal model. The red stars are the spaxels from the
imitated long-slit data of the spherical model. The solid line in each panel
is the demarcation line given by Kewley et al. (2001) to separate star-forming
galaxies/H ii regions from AGN. The dashed line in the [N ii]/H$\alpha$ versus
[O iii]/H$\beta$ panel is the demarcation line given by Kauffmann et al.
(2003) to separate star-forming galaxies/H ii regions from AGN. Figure 5: The
comparison of integrated emission-lines fluxes and ratios between the fractal
model and the spherical model. The emission-lines are classified into three
categories, which are the Balmer lines (grey shadow), the $Volume~{}species$
(blue shadow) and the $Boundary~{}species$ (red shadow). The black horizontal
dashed line indicates no changes of the line fluxes and ratios between the
fractal model and the spherical model.

We investigate the [N ii]/H$\alpha$, [S ii]/H$\alpha$, [O iii]/H$\beta$, [O
i]/H$\alpha$ and [O iii]/H$\beta$ diagnostics predicted from the fractal model
in comparison with the spherical model used in previous theoretical
classification schemes for the optical diagnostic diagrams. We collapse the 3D
distributions of emission-lines into 2D maps to imitate spatially-resolved
data.

Figure 4 shows the number density distribution of the spaxels from our model
on the optical diagnostic diagrams. On all diagnostic diagrams, the spaxels
lie on a sequence in the locus of H ii regions bounded by the demarcation
between star formation and harder ionizing sources given in Kewley et al.
(2001). The spaxel sequence shows a significant scatter on the diagnostic
diagrams with a scatter of 0.8 dex in the log([O iii]/H$\beta$) ratio, 0.6 dex
in the log([N ii]/H$\alpha$) ratio, 0.8 dex in the log([S ii]/H$\alpha$) ratio
and 1.2 dex in the log([O i]/H$\alpha$) ratio. The scatter is defined as the
region that includes 16% to 84% of the spaxels in the emission-line ratio
distributions. The scatter in the diagnostic line ratios is caused by the
complex internal structures of the ionization parameter and the hardness of
the ionizing radiation field, which vary from spaxel to spaxel, as shown in
Figure 2.

We collapse the spherical model along the x-axis and the y-axis to produce a
1D data to imitate a long-slit observation. Figure 4 shows that the spatially-
resolved spherical model lies on the spaxel sequence of the fractal model and
has agreement with the location on the standard optical diagnostic diagrams
where the spatially-resolved fractal model spaxels are concentrated. The
spherical model pixels agree with 19% of the fractal model spaxels on the [N
ii]/H$\alpha$ versus [O iii]/H$\beta$ diagram, 24% of the fractal model
spaxels on the [S ii]/H$\alpha$ versus [O iii]/H$\beta$ diagram, and 20% of
the fractal model spaxels on the [O ii]/H$\alpha$ versus [O iii]/H$\beta$
diagram. However, the spatially-resolved fractal model covers a larger range
of [N ii]/H$\alpha$, [S ii]/H$\alpha$, [O i]/H$\alpha$ and [O iii]/H$\beta$
ratios than the spatially-resolved spherical model.

###  3.4 Geometric Effects On Total Emission-Lines Fluxes Over A Nebula

Figure 5 gives the deviation between these two models between the spherical
and fractal models for the total (integrated) emission-line fluxes for each
line.

The complex geometry makes only a minor change in the fluxes of the Balmer
lines of 4% in the H$\alpha$ flux and 3% in the H$\beta$ flux. This result is
expected: Kennicutt & Evans (2012) previously found that the fluxes of Case B
recombination Balmer lines are solely correlated with the ionizing luminosity
of the central source with limited dependence on the nebular geometry.

The nebular geometry has a pronounced influence on the forbidden lines. We
classify the forbidden lines into “$Volume~{}species$” and
“$Boundary~{}species$” based on their locations within nebula (Figure 3).

“$Volume~{}species$” lines are the forbidden lines which are produced
throughout the entire nebula and can trace the properties of the overall
volume of the ISM within the nebula, [O iii] is an example of a volume-
sensitive emission-line. In this work, the total flux of the [O iii] line of
the fractal nebula model is 31% lower than the flux of the spherical model.
The [O iii]/H$\beta$ ratio of the fractal model is also 33% lower than the
ratio of the spherical model. The reduced [O iii] flux is caused by the log-
normal density distribution where 93% of the nebular volume are occupied by
the spaxels with density lower than the average density of 100 $\rm cm^{-3}$.

“$Boundary~{}species$” lines are the ionization lines residing on the boundary
of nebula and are sensitive to the perimeter of the nebular boundary. Complex
nebular geometry elongates the nebular boundary, increasing the total flux of
“$Boundary~{}species$” lines. In this work, the [N ii], [S ii] and [O i] lines
are the “$Boundary~{}species$” lines. Compared to the spherical model, the
complex geometry of the fractal model increases the total emission-line flux
by 29% for the [N ii] line, 161% for the [S ii] line and 253% for the [O i]
line. The total flux of the “$Boundary~{}species$” is positively correlated
with the amount of emission-lines concentrated on the edge of nebula. As shown
in Section 3.2, 99.5% of the [O i] emission, 87.5% of the [S ii] emission and
87.2% of the [N ii] emission is concentrated on the boundary. Therefore, the
[O i] line has the largest change in total flux from the spherical to the
fractal model among the three lines while the [N ii] has the least change in
total flux. Correspondingly, the line ratio is 24% larger for the [N
ii]/H$\alpha$ ratio, 151% larger for the [S ii]/H$\alpha$ ratio and 239%
larger for the [O i]/H$\alpha$ ratio for the fractals model, relative to the
1D model.

##  4 conclusions

We compute the first fractal-geometry photoionization model for an HII region
with 30 chemical elements to estimate the electron temperature, density and
emission-lines throughout the nebula in 3D. Our fractal model shows that the
density fluctuation of the turbulent ISM can cause the inhomogeneity of the
electron temperature and electron density, and lead to the complex ionization
structures seen in nebula. The fractal geometry more strongly impacts
emission-line fluxes of boundary elements [O i], [S ii] and [O iii]. The
emission-lines of [O iii] and the Balmer lines are only marginally affected by
nebular geometry.

We find that the spatially-resolved standard optical diagnostic diagrams of
the fractal model cannot be reconstructed by an equivalent spherical model,
but the average emission-line ratios of these two models have agreement within
0.64 dex for log([O iii]/H$\beta$) ratio, 0.16 dex for log([N ii]/H$\alpha$)
ratio, 0.13 dex for log([S ii]/H$\alpha$) ratio and 0.23 dex for log([O
i]/H$\alpha$) ratio. In conjunction with the complex nebular geometry, the
fluctuation of the ISM density within the fractal nebula causes a large
scatter in the optical diagnostic diagrams, corresponding to different
ionization parameter, density and temperature regions across an HII region.
The photoionization model with uniform density distribution and spherical
geometry is insufficient to interpret the spatially-resolved observations of
nearby nebula.

Geometric photoionization models are the future in interpreting emission-lines
from complex ISM conditions, especially for high-redshift galaxies which have
the evidently turbulent ISM. The ISM turbulence produces the complex H ii
region geometry and the density inhomogeneity within nebulae. The
perturbations from gas accretion and minor mergers in high-redshift galaxies
drive the turbulence of the ISM. Therefore, a more complex nebular geometry
driven by turbulence may play a role in the larger [N ii]/H$\alpha$ ratios
seen in high redshift galaxies.

This research was conducted on the traditional lands of the Ngunnawal and
Ngambri people. This research was supported by the Australian Research Council
Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
through project number CE170100013. YJ and LK gratefully acknowledge the
support of Lisa Kewley’s ARC Laureate Fellowship (FL150100113).

## References

  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Binette et al. (1985) Binette, L., Dopita, M. A., & Tuohy, I. R. 1985, ApJ, 297, 476. doi:10.1086/163544
  * Dopita et al. (2000) Dopita, M. A., Kewley, L. J., Heisler, C. A., et al. 2000, ApJ, 542, 224. doi:10.1086/309538
  * Dopita et al. (2006) Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006, ApJS, 167, 177. doi:10.1086/508261
  * Gendron-Marsolais et al. (2018) Gendron-Marsolais, M., Hlavacek-Larrondo, J., Martin, T. B., et al. 2018, MNRAS, 479, L28. doi:10.1093/mnrasl/sly084
  * Hillier & Miller (1998) Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407
  * Hillier & Miller (1999) Hillier, D. J., & Miller, D. L. 1999, ApJ, 519, 354
  * Hopkins & Beacom (2006) Hopkins, A. M. & Beacom, J. F. 2006, ApJ, 651, 142. doi:10.1086/506610
  * Hunt & Hirashita (2009) Hunt, L. K. & Hirashita, H. 2009, A&A, 507, 1327. doi:10.1051/0004-6361/200912020
  * Jin et al. (2022) Jin, Y., Kewley, L. J., & Sutherland, R. 2022, ApJ, 927, 37. doi:10.3847/1538-4357/ac48f3
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055. doi:10.1111/j.1365-2966.2003.07154.x
  * Kennicutt & Evans (2012) Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531. doi:10.1146/annurev-astro-081811-125610
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
  * Kewley & Dopita (2002) Kewley, L. J., & Dopita, M. A. 2002, ApJS, 142, 35
  * Kewley et al. (2013b) Kewley, L. J., Maier, C., Yabe, K., et al. 2013, ApJ, 774, L10
  * Kewley et al. (2013a) Kewley, L. J., Dopita, M. A., Leitherer, C., et al. 2013, ApJ, 774, 100
  * Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., Sutherland, R., et al. 2019, ApJ, 880, 16. doi:10.3847/1538-4357/ab16ed
  * Kobulnicky & Kewley (2004) Kobulnicky, H. A. & Kewley, L. J. 2004, ApJ, 617, 240. doi:10.1086/425299
  * Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3. doi:10.1086/313233
  * Lucy (1999) Lucy, L. B. 1999, A&A, 344, 282
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
  * McGaugh (1991) McGaugh, S. S. 1991, ApJ, 380, 140. doi:10.1086/170569
  * Moustakas et al. (2010) Moustakas, J., Kennicutt, R. C., Tremonti, C. A., et al. 2010, ApJS, 190, 233. doi:10.1088/0067-0049/190/2/233
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, by Donald E. Osterbrock. Published by University Science Books, ISBN 0-935702-22-9, 408pp, 1989.
  * Papovich et al. (2022) Papovich, C., Simons, R. C., Estrada-Carpenter, V., et al. 2022, arXiv:2205.05090
  * Rigby et al. (2011) Rigby, J. R., Wuyts, E., Gladders, M. D., Sharon, K., & Becker, G. D. 2011, ApJ, 732, 59
  * Rigby et al. (2021) Rigby, J. R., Florian, M., Acharyya, A., et al. 2021, ApJ, 908, 154. doi:10.3847/1538-4357/abcfc9
  * Rubin et al. (2011) Rubin, R. H., Simpson, J. P., O’Dell, C. R., et al. 2011, MNRAS, 410, 1320
  * Sanders et al. (2016) Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 816, 23. doi:10.3847/0004-637X/816/1/23
  * Shapley et al. (2005) Shapley, A. E., Coil, A. L., Ma, C.-P., & Bundy, K. 2005, ApJ, 635, 1006
  * Shirazi et al. (2014) Shirazi, M., Brinchmann, J., & Rahmati, A. 2014, ApJ, 787, 120
  * Smith et al. (2006) Smith, L. J., Westmoquette, M. S., Gallagher, J. S., et al. 2006, MNRAS, 370, 513. doi:10.1111/j.1365-2966.2006.10507.x
  * Snijders et al. (2007) Snijders, L., Kewley, L. J., & van der Werf, P. P. 2007, ApJ, 669, 269. doi:10.1086/521522
  * Steidel et al. (2014) Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165
  * Simón-Díaz & Stasińska (2008) Simón-Díaz, S. & Stasińska, G. 2008, MNRAS, 389, 1009. doi:10.1111/j.1365-2966.2008.13444.x
  * Strom et al. (2022) Strom, A. L., Rudie, G. C., Steidel, C. C., et al. 2022, ApJ, 925, 116. doi:10.3847/1538-4357/ac38a3
  * Sutherland & Dopita (1993) Sutherland, R. S. & Dopita, M. A. 1993, ApJS, 88, 253. doi:10.1086/191823
  * Sutherland & Bicknell (2007) Sutherland, R. S. & Bicknell, G. V. 2007, ApJS, 173, 37. doi:10.1086/520640
  * Tacchella et al. (2022) Tacchella, S., Smith, A., Kannan, R., et al. 2022, MNRAS. doi:10.1093/mnras/stac818
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898
  * Veilleux & Osterbrock (1987) Veilleux, S., & Osterbrock, D. E. 1987, ApJS, 63, 295
  * Yuan et al. (2013) Yuan, T.-T., Kewley, L. J., & Richard, J. 2013, ApJ, 763, 9
  * Zahid et al. (2011) Zahid, H. J., Kewley, L. J., & Bresolin, F. 2011, ApJ, 730, 137
  * Zavala et al. (2022) Zavala, S., Toalá, J. A., Santamaría, E., et al. 2022, MNRAS. doi:10.1093/mnras/stac1097