# 4$\alpha$ linear-chain state produced by 9Be+9Be collision

Tomoyuki Baba1, Yasutaka Taniguchi2,6 Masaaki Kimura3,4,5,6 1Kitami Institute
of Technology, 090-8507 Kitami, Japan
2Department of Information Engineering, National Institute of Technology
(KOSEN), Kagawa College, 769-1192 Mitoyo, Japan
3Department of Physics, Hokkaido University, 060-0810 Sapporo, Japan
4RIKEN Nishina Center, 351-0198 Wako, Japan.
5Reaction Nuclear Data Centre (JCPRG), Hokkaido University, 060-0810 Sapporo,
Japan
6Research Center for Nuclear Physics (RCNP), Osaka University, 567-0047
Ibaraki, Japan

###### Abstract

Extreme nuclear deformations provide great insight into the geometric
formation of quantum many-body systems. In this work, the $4\alpha+2n$ linear
chain is assessed in 18O. We predict excitation energies, moment-of-inertia,
$\alpha$-, and 9Be-decay widths by using the antisymmetrized molecular
dynamics. We show that the $K^{\pi}=3^{-}$ linear-chain states may be verified
by the head-on ${}^{9}{\rm Be}+{}^{9}{\rm Be}$ collision experiments.

Introduction—Nuclear clustering has many similarities to molecules. For
instance, extra neutrons surrounding $\alpha$-particles play a glue-like role
analogous to the covalent electrons of atomic molecules as in beryllium
($2\alpha+xn$) seya81 ; oertzen96 ; enyo99 ; 12Beex1 ; itagaki00 ; itoBe ;
ito12Be , carbon ($3\alpha+xn$) itag01 , and in the heavier mass isotopes
furu08 ; kimu07 ; tani14 . In the carbon isotopes ($3\alpha+xn$), it has been
expected that the extremely deformed states which have the intrinsic structure
of linearly aligned three $\alpha$-particles, called as linear-chain (LC)
structure, will be stabilized by the assist of the glue-like role of the
covalent neutrons. Recent experimental free14 ; frit16 ; tian16 ; yama17 ;
li17  and theoretical suha10 ; ebra14 ; baba16 ; baba17 ; enyo20 ; ren20
studies have identified the linear-chain states in 14C($3\alpha+2n$), and the
research of 16C ($3\alpha+4n$) has also been making progress baba18 ; liu20 ;
baba20 .

A fascinating and long-standing open question is how many $\alpha$-particles
can compose the linear-chain structure. Unlike the 3$\alpha$ linear chains,
linear chains containing four or more $\alpha$-particles have not been
identified. In order to produce the $4\alpha$ linear-chain state in oxygen
isotopes, the ABe+ABe resonant scattering may be a natural way. Although the
formation of the $4\alpha$ linear chain in 16O was predicted by several
theoretical works suzu72 ; ichi11 ; suha14 ; yao14 ; suzu17 ; inak18 , 8Be is
the unbound nucleus, so that it is not easy to prove its existence from the
8Be+8Be reaction. Given that 9Be is the only stable beryllium isotope, 9Be+9Be
scattering is the most feasible way to confirm the $4\alpha+2n$ linear chain.
Therefore, we study the excited states of 18O to search for the candidate of
$4\alpha+2n$ linear chain.

Linear chain produced by 9Be+9Be collision—Let us consider the $4\alpha$
linear-chain configurations which will be produced by the head-on 9Be+9Be
collisions. For simplicity, we approximate the ground state of 9Be ($3/2^{-}$)
as the $\alpha+\alpha+n$ system with a valence neutron occupying the $p_{3/2}$
orbit with $j_{z}=\pm 3/2$. As illustrated in Fig. 1, there are two ways to
linearly align two 9Be, which yields different valence neutron configurations;
(a) the anti-parallel and (b) the parallel alignments with respect to the
valence neutron’s $j_{z}$. The anti-parallel alignment yields the intrinsic
state with $K=0$, where $K$ denotes the $z$-component of the intrinsic angular
momentum equal to the sum of valence neutron’s $j_{z}$. Because this
configuration is an admixture of the positive- and negative-parity states, we
expect that it leads to a pair of the rotational bands;
$J^{\pi}=0^{+},2^{+},4^{+},...$ and $1^{-},3^{-},5^{-},...$ bands. Indeed, we
obtain both bands in this work. However, we focus on the only
$J^{\pi}=0^{+},2^{+},4^{+},...$ band because the negative-parity band is
located at higher than the positive-parity band.

The parallel alignment (Fig. 1 (b)) yields the intrinsic state with $K=3$. The
parity of this configuration is uniquely determined as negative due to the
following reason. The spin wave function is symmetric ($S=1$) as both neutron
spins are aligned to $s_{z}=1/2$. Because the isospin wave function is also
symmetric, the spatial wave function must be asymmetric with respect to the
exchange of two valence neutrons (parity transformation). Thus, this
configuration forms a $J^{\pi}=3^{-},4^{-},5^{-},...$ band. In short, the
anti-parallel alignment yields a pair of the positive- and negative-parity
bands with $K=0$, whereas the parallel alignment yields a negative-parity band
with $K=3$.

Figure 1: (color online) Schematic illustration of the $4\alpha$ linear chains
constructed by (a) anti-parallel and (b) parallel alignments with respect to
valence neutron’s $j_{z}$ of 9Be. The black arrows indicate $j_{z}=+3/2$ orbit
of the valence neutron while the red arrow indicates $j_{z}=-3/2$.

Calculated properties of the linear-chain states—In order to describe the
linear-chain states, we use the antisymmetrized molecular dynamics (AMD). We
employ the Hamiltonian with the Gogny D1S nucleon-nucleon interaction gogn91 .
The AMD wave function $\Psi^{\pi}_{\rm AMD}$ is a parity-projected Slater
determinant of single particle wave packets,

$\displaystyle\Psi^{\pi}_{\rm AMD}=\hat{P}^{\pi}\Psi_{\rm
AMD}=\hat{P}^{\pi}{\mathcal{A}}\\{\psi_{1},\psi_{2},...,\psi_{A}\\}.$ (1)

Here, $\hat{P}^{\pi}$ is the parity-projection operator, and $\psi_{i}$ is the
single particle wave packet which is a direct product of the deformed Gaussian
for the spatial part, spin ($\chi_{i}$) and isospin ($\xi_{i}$) parts kimu04 ,

$\displaystyle\phi_{i}({\bm{r}})$
$\displaystyle=\prod_{\sigma=x,y,z}\exp\biggl{\\{}-\nu_{\sigma}\Bigl{(}r_{\sigma}-\frac{Z_{i\sigma}}{\sqrt{\nu_{\sigma}}}\Bigr{)}^{2}\biggr{\\}}\otimes\chi_{i}\otimes\xi_{i},$
(2) $\displaystyle\chi_{i}$
$\displaystyle=a_{i}\chi_{\uparrow}+b_{i}\chi_{\downarrow},\quad\xi_{i}={\rm
proton}\quad{\rm or}\quad{\rm neutron}.$

The centroids of the Gaussian wave packets $\bm{Z}_{i}$, the direction of
nucleon spin $a_{i},b_{i}$, and the width parameter $\nu_{\sigma}$ are the
variables determined by the frictional cooling method enyo95 . In this study,
we impose the constraint on the quadrupole deformation parameter $\beta$ to
describe extremely deformed $4\alpha$ linear chain. After the variational
calculation, the eigenstate of the total angular momentum $J$ is projected
out. We perform the generator coordinate method hill54  by employing the
quadrupole deformation parameter $\beta$ as the generator coordinate.

In our previous work baba19 , it has been shown that the AMD plausibly
describes the low-lying states of 18O. The binding energy of 18O is calculated
as 139.97 MeV whereas the observed value is 139.81 MeV. The low-lying excited
states including the ${}^{14}{\rm C}+{}^{4}{\rm He}$ cluster states are also
reasonably described. Therefore, we expect that the AMD can also precisely
describe higher-lying states of 18O.

Figure 2: (color online) Calculated energies above ${}^{9}{\rm Be}+{}^{9}{\rm
Be}$ threshold of the linear-chain states as a function of angular momenta.
The energy is relative to 23.64 MeV.

Note that the present calculation does not assume a prior linear-chain
configuration. In fact, we obtain many excited states with various cluster and
non-cluster states. Among these excited states, we assign two rational bands
as the linear-chain candidates shown in Fig. 2; a positive-parity band built
on the $0^{+}$ state followed by $2^{+},4^{+},...$ states and a negative-
parity band built on the $3^{-}$ state followed by $4^{-},5^{-},...$ states.
We remark that the spin-parity of these bands is consistent with that expected
from the ${}^{9}{\rm Be}+{}^{9}{\rm Be}$ configurations discussed above. The
reasons for this assignment are as follows. The member states of each band
have the same intrinsic structure. All positive-parity band members have a
large overlap with the intrinsic wave function shown in Fig. 3(a). Its proton
density distribution shows the linear alignment of four $\alpha$ particles,
and the two valence neutrons occupy the negative-parity orbits with $j_{z}=\pm
3/2$. This intrinsic structure is approximately the anti-parallel
configuration shown in Fig. 1(a). However, differently from Fig. 1(a), the two
valence neutrons are localized around two $\alpha$-particles at the center
because of the attraction between valence neutrons and no Pauli exclusion.
Therefore, the total system forms the $\alpha+{}^{10}{\rm Be}+\alpha$ -like
structure.

The intrinsic state of the negative-parity band is shown in Fig. 3 (b1) and
(b2). The proton density distribution indicates that this band also has
4$\alpha$ linear-chain core. A valence neutron occupy the negative-parity
orbit with $j_{z}=3/2$ (Fig. 3(b1)), and the other occupies the positive-
parity orbit with $j_{z}=3/2$ (Fig. 3(b2)). These single-particle orbits can
be understood as a linear combination of the $p_{3/2}$ orbits of two 9Be. Let
us denote the $p_{3/2}$ ($j_{z}=3/2$) orbit of the left (right) side 9Be as
$\varphi_{L}$ ($\varphi_{R}$), which are schematically illustrated as black
arrows in Fig. 1(b). Then, the single-particle orbits are represented as,

$\displaystyle\varphi_{\pm}=\frac{1}{\sqrt{2}}(\varphi_{L}\pm\varphi_{R}).$

They generate an orthogonalized pair of the negative- and positive-parity
orbits with $j_{z}=3/2$, which corresponds to Fig. 3 (b1) and (b2),
respectively. This intrinsic state corresponds to the parallel configuration
in Fig. 1. In contrast to the anti-parallel configuration, the two valence
neutrons separately locate the left and right side Be because of the Pauli
principle which results from their parallel spin. We find that only these two
bands have the structure corresponding to the ${}^{9}{\rm Be}+{}^{9}{\rm Be}$
head-on collision in the vicinity of the energies near the ${}^{9}{\rm
Be}+{}^{9}{\rm Be}$ threshold.

Figure 3: (color online) Density distributions of the intrinsic states of the
linear-chain bands. Contour lines show the proton density distributions and
color plots show the valence neutron orbits. Panel (a) shows the intrinsic
state of the positive-parity band, in which two valence neutrons occupy the
negative-parity orbits with $j_{z}=\pm 3/2$. Panel (b1) and (b2) show the
intrinsic state of the negative-parity band, in which a valence neutron
occupies the negative-parity orbit with $j_{z}=3/2$ (b1), and the other
occupies the positive-parity orbit with $j_{z}=3/2$ (b2).

These two bands have strongly deformed intrinsic shapes compatible with the
$4\alpha$ linear chain. The quadrupole deformation parameters of the positive-
and negative-parity intrinsic states are equally $\beta=1.34$. Consequently,
they have enormous moment-of-inertia as large as $\hbar/2\Im=55$ keV for the
positive-parity band and $58$ keV for the negative-parity band. These are even
larger than a rigid-rotor estimation, $\hbar/2\Im=85$ keV estimated as
follows. The classical moment-of-inertia of the prolate spheroid with 2:1,
3:1, and 4:1 deformation satisfy the following relation,

$\displaystyle\Im_{2:1}:\Im_{3:1}:\Im_{4:1}=5:15:34.$

Since the moment-of-inertia of the $2\alpha$ and $3\alpha$ linear chains are
measured as $\hbar/2\Im=590$ keV ($2\alpha+2n$) till04  and $\hbar/2\Im=190$
keV in ($3\alpha+2n$) yama17 , the moment-of-inertia of $4\alpha$ linear chain
is estimated as $\hbar/2\Im=85$ keV.

Table 1: Excitation energies, energies above ${}^{9}{\rm Be}+{}^{9}{\rm Be}$ threshold, partial $\alpha$-decay widths and dimensionless reduced widths of the $4\alpha$ linear-chain states. 14C is assumed to be the $0^{+}$ and $2^{+}$ states of the linear-chain band. The channel radius is $7.0$ fm. $J^{\pi}$ |  $E_{x}$ [MeV] |  $E_{{\rm Be}}$ [MeV] |  $\Gamma_{\alpha}(0^{+}$; LC) [keV] |  $\theta^{2}_{\alpha}(0^{+}$; LC) $[\times 10^{-2}]$ |  $\Gamma_{\alpha}(2^{+}$; LC) [keV] |  $\theta^{2}_{\alpha}(2^{+}$; LC) $[\times 10^{-2}]$
---|---|---|---|---|---|---
$0^{+}$ | 27.85 | 4.21 | 382 | 8.21 | 1136 | 31.8
$2^{+}$ | 28.32 | 4.68 | 346 | 7.70 | 1044 | 13.4 ($l=4$)
$4^{+}$ | 29.37 | 5.73 | 288 | 3.51 | 918 | 5.37 ($l=2$)
$6^{+}$ | 30.16 | 6.51 | 160 | 5.39 | 572 | 10.2 ($l=4$)
$3^{-}$ | 31.09 | 7.45 | 0 | 0.00 | 0 | 0.00
$4^{-}$ | 31.46 | 7.82 | 0 | 0.00 | 0 | 0.00
$5^{-}$ | 31.86 | 8.22 | 0 | 0.00 | 0 | 0.00
$6^{-}$ | 32.58 | 8.95 | 0 | 0.00 | 0 | 0.00
$7^{-}$ | 33.68 | 10.0 | 0 | 0.00 | 0 | 0.00

Decay properties of the linear-chain states—In order to take a deep dive on
the difference of the two linear-chain bands, we discuss their decay
properties. Calculated excitation energies, partial $\alpha$-decay widths, and
dimensionless reduced widths are listed in Table. 1. A dimensionless reduced
width $\theta^{2}_{l}(a)$ is defined by the ratio of the reduced width to
Wigner limit,

$\displaystyle\theta^{2}_{l}(a)=\frac{a}{3}|ay_{l}(a)|^{2},$ (3)

where $y_{l}(a)$ is the reduced width amplitude,

$\displaystyle y_{l}(r)=\sqrt{\frac{A!}{A_{\rm C1}!A_{\rm
C2}!}}\langle\phi_{\rm C1}[\phi_{\rm
C2}Y_{l0}({\hat{r}})]_{J^{\pi}M}|\Psi^{J\pi}_{Mn}\rangle,$ (4)

In the $\alpha+{}^{14}{\rm C}$ channel, the daughter nucleus 14C is assumed to
be the $0^{+}$ and $2^{+}$ states of the $3\alpha$ linear-chain states of 14C
studied in Ref.baba16 . The reduced width in the $\alpha+{}^{14}{\rm C}$(g.s.)
is almost zero due to the extreme deformation of the $4\alpha$ linear chain.
The $\alpha$-decay properties are different between the $K^{\pi}=0^{+}$ and
$K^{\pi}=3^{-}$ linear chains. The $K^{\pi}=0^{+}$ linear-chain states have
very large $\alpha$-decay widths into the linear-chain states of 14C. In
particular, the widths to the $\alpha+{}^{14}{\rm C}(2^{+}$; LC) channel are
huge because of the strong angular correlation between the linearly aligned
$\alpha$-particles. This character is in contrast to the Hoyle state where
$\alpha$ particles are weakly bound with $l=0$ and hence, the ${}^{8}{\rm
Be}(0^{+}_{1})$ component dominates funa15 . The large reduced width amplitude
for $J^{\pi}=0^{+},2^{+},...$ is a significant linear-chain structure feature.
As the 14C orientation is fixed, it does not become an eigenstate but a mixed
state of angular momentum. The similar character was also discussed for the
$3\alpha$ linear chain of the carbon isotopes suzu72 ; baba16 ; baba18 . The
$K^{\pi}=3^{-}$ linear-chain states have almost zero $\alpha$-decay widths
into the linear chain of 14C. It is clear because the both valence neutrons in
the linear-chain states of 14C are negative parity, as a result they are
orthogonal with the linear chain of 18O with the positive- and negative-parity
valence neutrons shown in Fig. 3 (b). In addition, we did not find the
negative-parity $3\alpha$ linear-chain states of 14C with the positive- and
negative-parity valence neutrons baba16 . In present calculation, therefore,
the $K^{\pi}=3^{-}$ linear-chain states do not decay to any
$\alpha+{}^{14}{\rm C}$ channels.

Table 2: Partial ${}^{9}{\rm Be}$-decay widths and dimensionless reduced widths of the $4\alpha$ linear-chain states. ${}^{9}{\rm Be}(3/2^{-})+{}^{9}{\rm Be}(3/2^{-})$ and ${}^{9}{\rm Be}(3/2^{-})+{}^{9}{\rm Be}(5/2^{-})$ are assumed. The channel radius is $7.0$ fm. $J^{\pi}$ |  $E_{x}$ [MeV] |  $E_{{\rm Be}}$ [MeV] |  $\Gamma_{\rm Be}(3/2^{-})$ [keV] | $\theta^{2}_{\rm Be}(3/2^{-})[\times 10^{-2}]$ |  $\Gamma_{\rm Be}(5/2^{-})$ [keV] | $\theta^{2}_{\rm Be}(5/2^{-})[\times 10^{-2}]$
---|---|---|---|---|---|---
$0^{+}$ | 27.85 | 4.21 | 6 | 0.31 | 9 | 0.52 ($l=2$)
$2^{+}$ | 28.32 | 4.68 | 6 | 0.30 | 8 | 0.23 ($l=4$)
$4^{+}$ | 29.37 | 5.73 | 6 | 0.14 | 8 | 0.10 ($l=2$)
$6^{+}$ | 30.16 | 6.51 | 4 | 0.32 | 7 | 0.25 ($l=4$)
$3^{-}$ | 31.09 | 7.45 | 264 | 5.37 ($l=1$) | 253 | 1.84 ($l=1$)
$4^{-}$ | 31.46 | 7.82 | 151 | 3.02 ($l=3$) | 144 | 0.92 ($l=3$)
$5^{-}$ | 31.86 | 8.22 | 335 | 5.81 ($l=3$) | 278 | 2.21 ($l=3$)
$6^{-}$ | 32.58 | 8.95 | 307 | 6.33 ($l=5$) | 245 | 2.10 ($l=5$)
$7^{-}$ | 33.68 | 10.0 | 281 | 4.87 ($l=5$) | 251 | 1.51 ($l=5$)

Partial ${}^{9}{\rm Be}$-decay widths, and dimensionless reduced widths are
listed in Table. 2. We consider the ${}^{9}{\rm Be}(3/2^{-})+{}^{9}{\rm
Be}(3/2^{-})$, ${}^{9}{\rm Be}(3/2^{-})+{}^{9}{\rm Be}(5/2^{-})$ channels,
where ${}^{9}{\rm Be}(3/2^{-})$ and ${}^{9}{\rm Be}(5/2^{-})$ are the ground
and excited states of ${}^{9}{\rm Be}$. Characters of the decay widths are
different between the $K^{\pi}=0^{+}$ and $K^{\pi}=3^{-}$ linear chains. The
$K^{\pi}=0^{+}$ linear-chain states have very small widths into the
${}^{9}{\rm Be}+{}^{9}{\rm Be}$ channel, which is consistent with the
$\alpha+{}^{10}{\rm Be}+\alpha$ picture shown in Fig. 3 (a). On the other
hand, the $K^{\pi}=3^{-}$ linear-chain states show large ${}^{9}{\rm Be}$
widths. Therefore, the $K^{\pi}=3^{-}$ linear chain manifests so strong
9Be+9Be correlation shown in Fig. 1 (b) and these linear-chain states can be
observed by a 9Be+9Be collision.

In order to clarify the characteristic ${}^{9}{\rm Be}$-decay modes, we
calculate overlaps between the Brink and the AMD wave functions defined as,

$\displaystyle\mathcal{O}(r)=\frac{|\langle\Phi^{K\pi}_{\rm
BB}(r)|\hat{P}^{J}_{KK}|\Psi^{\pi}_{\rm
AMD}\rangle|^{2}}{|\langle\Phi^{K\pi}_{\rm
BB}(r)|\hat{P}^{J}_{KK}|\Phi^{K\pi}_{\rm BB}(r)\rangle||\langle\Psi^{\pi}_{\rm
AMD}|\hat{P}^{J}_{KK}|\Psi^{\pi}_{\rm AMD}\rangle|},$ (5)

where $\hat{P}^{J}_{KK}$ is the angular momentum projection operator. The
Brink wave functions $\Phi^{K\pi}_{\rm BB}(r)$ are constructed by the linearly
alignment of two 9Be shown in Fig. 1; (a) anti-parallel $\Phi^{0+}_{\rm
BB}(r)$ and (b) parallel $\Phi^{3-}_{\rm BB}(r)$ alignments.

$\displaystyle\Phi^{0+}_{\rm BB}(r)$
$\displaystyle=\hat{P}^{\pi}\mathcal{A}\\{\phi^{jz=3/2}_{\rm
Be}(-r/2)\phi^{jz=-3/2}_{\rm Be}(r/2)\\},$ (6) $\displaystyle\Phi^{3-}_{\rm
BB}(r)$ $\displaystyle=\hat{P}^{\pi}\mathcal{A}\\{\phi^{jz=3/2}_{\rm
Be}(-r/2)\phi^{jz=3/2}_{\rm Be}(r/2)\\}.$ (7)

Here, the wave function of 9Be is described as,

$\displaystyle\phi^{jz}_{\rm Be}$
$\displaystyle=\mathcal{A}\\{\phi_{\alpha}\otimes\phi_{\alpha}\otimes(0p_{3/2})\\},$
(8)

where $(0p_{3/2})$ is represented by a infinitesimally shifted Gaussian wave
packet based on the antisymmetrized quasicluster model suhajj . Figure 4 shows
the calculated overlap as a function of the distance $r$. The $K^{\pi}=3^{-}$
linear chain has a large amount of overlap 0.84 with the Brink wave function
$\Phi^{3-}_{\rm BB}$($r=6.5$ fm) shown in Fig. 1 (b). On the other hand, the
$K^{\pi}=0^{+}$ linear chain has a amount of overlap 0.52 with $\Phi^{0+}_{\rm
BB}$($r=5.0$ fm) shown in Fig. 1 (a). Compared to the $K^{\pi}=3^{-}$ linear
chain, the 9Be+9Be correlation is small in the $K^{\pi}=0^{+}$ linear chain.
The $K^{\pi}=3^{-}$ and $K^{\pi}=0^{+}$ states have different features in the
outer region. The $K^{\pi}=3^{-}$ state has more considerable overlap in the
outer region than the $K^{\pi}=0^{+}$ state, with the outer peak position. At
the present channel radius of $7.0$ fm, the overlap of $K^{\pi}=3^{-}$ is one
order of magnitude larger than that of $K^{\pi}=0^{+}$, which leads to the
much larger 9Be+9Be decay width. Indeed, this difference is reflected in the
difference of the 9Be-reduced widths in Table. 2. The $K^{\pi}=0^{+}$ linear
chain has the weak 9Be+9Be correlation so that the 9Be-reduced widths are
small, while the $K^{\pi}=3^{-}$ linear chain has so strong 9Be+9Be
correlation that the member states have large 9Be-reduced widths. Therefore,
we conclude that the $K=3^{-}$ linear-chain states of 18O can be observed by a
9Be+9Be resonant scattering shown in Fig. 1 (b).

Figure 4: (color online) Calculated overlaps between the 9Be+9Be Brink wave
functions and the AMD wave functions. $\Psi^{+}_{LC}$ and $\Psi^{-}_{LC}$
correspond to Fig. 3 (a) and (b), respectively.

Summary—In summary, we have presented the first assessment of the $4\alpha+2n$
linear-chain configuration in 18O using the AMD calculation. There are two
different $4\alpha$ linear-chain bands, $K^{\pi}=0^{+}$ and $K^{\pi}=3^{-}$.
We predict their excitation energies, moment-of-inertia, $\alpha$-, and 9Be-
decay widths. In both bands, the moment-of-inertia is rather large, which is a
strong evidence for the extreme deformation. The two types of linear chains
show different decay properties. The $K^{\pi}=0^{+}$ linear-chain states show
large decay widths into the $\alpha+{}^{14}{\rm C}$(LC) channel while the
$K^{\pi}=3^{-}$ linear-chain states show large those into the ${}^{9}{\rm
Be}+{}^{9}{\rm Be}$ channel. In order to clarify this difference, we calculate
the overlaps with the 9Be+9Be Brink wave functions. As the results, the
$K^{\pi}=0^{+}$ linear chain has small overlap while the $K^{\pi}=3^{-}$
linear chain has large overlap. It means that the $K^{\pi}=3^{-}$ linear chain
shows a strong 9Be+9Be correlation. Therefore, we expect that the linear chain
of 18O can be produced by the head-on 9Be+9Be collision. We believe that these
are promising properties that can be investigated in future experiments and
establish the existence of the exotic $4\alpha$ linear chain.

We thank Prof. M. Ito and Dr. H. Yamaguchi for useful discussions. This
calculation has been done on a supercomputer at Research Center for Nuclear
Physics, Osaka University. The authors acknowledge the support of the
collaborative research program 2021 at Hokkaido University. This work was
supported by the JSPS KAKENHI Grant No. 19K03859.

## References

  * (1) M. Seya, M. Kohno, and S. Nagata, Prog. Theor. Phys. 65, 204 (1981).
  * (2) W. von Oertzen, Z. Phys. A 354, 37 (1996); ibid. 357, 355 (1997).
  * (3) Y. Kanada-En’yo, H. Horiuchi, and A. Dote, Phys. Rev. C 60, 064304 (1999).
  * (4) M. Freer, J. C. Angelique, L. Axelsson, B. Benoit, U. Bergmann, W. N. Catford, S. P. G. Chappell, N. M. Clarke, N. Curtis, A. D’Arrigo, E. de Goes Brennard, O. Dorvaux, B. R. Fulton, G. Giardina, C. Gregori, S. Grevy, F. Hanappe, G. Kelly, M. Labiche, C. Le Brun, S. Leenhardt, M. Lewitowicz, K. Markenroth, F. M. Marques, M. Motta, J. T. Murgatroyd, T. Nilsson, A. Ninane, N. A. Orr, I. Piqueras, M. G. Saint Laurent, S. M. Singer, O. Sorlin, L. Stuttge, and D. L. Watson, Phys. Rev. Lett. 82, 1383 (1999).
  * (5) N. Itagaki and S. Okabe, Phys. Rev. C 61, 044306 (2000).
  * (6) M. Ito, N. Itagaki, H. Sakurai, and K. Ikeda, Phys. Rev. Lett. 100, 182502 (2008).
  * (7) M. Ito, Phys. Rev. C 85, 044308 (2012).
  * (8) N. Itagaki, S. Okabe, K. Ikeda and I. Tanihata, Phys. Rev. C 64, 014301 (2001).
  * (9) N. Furutachi, M. Kimura, A. Dot’e, Y. Kanada-En’yo and S. Oryu, Prog. Theor. Phys. 119, 3 (2008).
  * (10) M. Kimura, Phys. Rev. C 75, 034312 (2007).
  * (11) Y. Taniguchi, Phys. Rev. C 90, 054308 (2014).
  * (12) M. Freer, J. D. Malcolm, N. L. Achouri, N. I. Ashwood, D. W. Bardayan, S. M. Brown, W. N. Catford, K. A. Chipps, J. Cizewski, N. Curtis, K. L. Jones, T. Munoz-Britton, S. D. Pain, N. Soic, C. Wheldon, G. L. Wilson, and V. A. Ziman, Phys. Rev. C 90, 054324 (2014).
  * (13) A. Fritsch, S. Beceiro-Novo, D. Suzuki, W. Mittig, J.J. Kolata, T. Ahn, D. Bazin, F. D. Becchetti, B. Bucher, Z. Chajecki, X. Fang, M. Febbraro, A. M. Howard, Y. Kanada-En’yo, W. G. Lynch, A. J. Mitchell, M. Ojaruega, A. M. Rogers, A. Shore, T. Suhara, X. D. Tang, R. Torres-Isea, and H. Wang, Phys. Rev. C 93, 014321 (2016).
  * (14) Z. Y. Tian, Y. L. Ye, Z. H. Li, C. J. Lin, Q. T. Li, Y. C. Ge, J. L. Lou, W. Jiang, J. Li, Z. H. Yang, J. Feng, P. J. Li, J. Chen, Q. Liu, H. L. Zang, B. Yang, Y. Zhang, Z. Q. Chen, Y. Liu, X. H. Sun, J. Ma, H. M. Jia, X. X. Xu, L. Yan, N. R. Ma, and L. J. Sun, Chinese Phys. C 40, 111001 (2016).
  * (15) H. Yamaguchi, D. Kahl, S. Hayakawa, Y. Sakaguchi, K. Abe, T. Nakao, T. Suhara, N. Iwasa, A. Kim, D. H. Kim, S. M. Cha, M. S. Kwag, J. H. Lee, E. J. Lee, K. Y. Chae, Y. Wakabayashi, N. Imai, N. Kitamura, P. Lee, J. Y. Moon, K. B. Lee, C. Akers, H. S. Jung, N. N. Duy, L. H. Khiem, and C. S. Lee, Phys. Lett. B 766 (2017) 11-16.
  * (16) J. Li, Y. L. Ye, Z. H. Li, C. J. Lin, Q. T. Li, Y. C. Ge, J. L. Lou, Z. Y. Tian, W. Jiang, Z. H. Yang, J. Feng, P. J. Li, J. Chen, Q. Liu, H. L. Zang, B. Yang, Y. Zhang, Z. Q. Chen, Y. Liu, X. H. Sun, J. Ma, H. M. Jia, X. X. Xu, L. Yang, N. R. Ma, and L. J. Sun, Phys. Rev. C 95, 021303(R) (2017).
  * (17) T. Suhara and Y. Kanada-En’yo, Phys. Rev. C 82, 044301 (2010).
  * (18) J. -P. Ebran, E. Khan, T. Niksic, and D. Vretenar, Phys. Rev. C 90, 054329 (2014).
  * (19) T. Baba and M. Kimura Phys. Rev. C 94, 044303 (2016).
  * (20) T. Baba and M. Kimura Phys. Rev. C 95, 064318 (2017).
  * (21) Y. Kanada-En’yo and K. Ogata, Phys. Rev. C 101, 014317 (2020).
  * (22) Z. X. Ren, P. W. Zhao, and J. Meng, Phys. Lett. B 801, 135194 (2020).
  * (23) T. Baba and M. Kimura Phys. Rev. C 97, 054315 (2018).
  * (24) Y. Liu, Y. L. Ye, J. L. Lou, X. F. Yang, T. Baba, M. Kimura, B. Yang, Z. H. Li, Q. T. Li, J. Y. Xu, Y. C. Ge, H. Hua, J. S. Wang, Y. Y. Yang, P. Ma, Z. Bai, Q. Hu, W. Liu, K. Ma, L. C. Tao, Y. Jiang, L. Y. Hu, H. L. Zang, J. Feng, H. Y. Wu, J. X. Han, S. W. Bai, G. Li, H. Z. Yu, S. W. Huang, Z. Q. Chen, X. H. Sun, J. J. Li, Z. W. Tan, Z. H. Gao, F. F. Duan, J. H. Tan, S. Q. Sun, and Y. S. Song, Phys. Rev. Lett. 124, 192501 (2020).
  * (25) T. Baba, Y. Liu, J. X. Han, and Y. L. Ye, Phys. Rev. C 102, 041302(R) (2020).
  * (26) Y. Suzuki, H. Horiuchi, and K. Ikeda, Prog. Theor. Phys. 47, 5 (1972).
  * (27) T. Ichikawa, J. A. Maruhn, N. Itagaki, and S. Ohkubo, Phys. Rev. Lett. 107, 112501 (2011).
  * (28) T. Suhara, Y. Funaki, B. Zhou, H. Horiuchi, and A. Tohsaki, Phys. Rev. Lett. 112, 062501 (2014).
  * (29) J. M. Yao, N. Itagaki, and J. Meng, Phys. Rev. C 90, 054307 (2014).
  * (30) Y. Suzuki and W. Horiuchi, Phys. Rev. C 95, 044320 (2017).
  * (31) T. Inakura and S. Mizutori, Phys. Rev. C 98, 044312 (2018).
  * (32) J. F. Berger, M. Girod, and D. Gogny, Comput. Phys. Comm. 63 (1991) 365.
  * (33) M. Kimura, Phys. Rev. C 69, 044319 (2004).
  * (34) Y. Kanada-En’yo and H. Horiuchi, Prog. Theor. Phys. 93, 115 (1995).
  * (35) D. L. Hill and J. A. Wheeler, Phys. Rev. 89, 1102 (1953).
  * (36) T. Baba and M. Kimura, Phys. Rev. C 100, 064311 (2019).
  * (37) D. R. Tilley, J. H. Kelley, J. L. Godwin, D. J. Millener, J. E. Purcell, C. G. Sheu, H. R. Weller, Nucl. Phys. A 745, 155-362 (2004).
  * (38) Y. Funaki, H. Horiuchi, and A. Tohsaki, Prog. Part. Nucl. Phys. 82 78-132 (2015).
  * (39) T. Suhara, N. Itagaki, J. Cseh, and M. Ploszajczak, Phys. Rev. C 87, 054334 (2013).


# Disorder Enhances the Fracture Toughness of Mechanical Metamaterials

Sage Fulco,1,∗ Michal K. Budzik2, Hongyi Xiao,3,4, Douglas J. Durian3,1,
Kevin T. Turner1,∗
1Department of Mechanical Engineering and Applied Mechanics,
University of Pennsylvania, Philadelphia, PA 19104
2Department of Mechanical Engineering and Production,
Aarhus University, Aarhus, Denmark
3Department of Physics and Astronomy,
University of Pennsylvania, Philadelphia, PA 19104
4Department of Mechanical Engineering,
University of Michigan, Ann Arbor, MI 48109

∗To whom correspondence should be addressed;
E-Mails<EMAIL_ADDRESS>(S. Fulco<EMAIL_ADDRESS>(K.T. Turner)

> Mechanical metamaterials with engineered failure properties typically rely
> on periodic unit cell geometries or bespoke microstructures to achieve their
> unique properties. We demonstrate that intelligent use of disorder in
> metamaterials leads to distributed damage during failure, resulting in
> enhanced fracture toughness with minimal losses of strength. Toughness
> depends on the level of disorder, not a specific geometry, and the confined
> lattices studied exhibit a maximum toughness enhancement at an optimal level
> of disorder. A mechanics model that relates disorder to toughness without
> knowledge of the crack path is presented. The model is verified through
> finite element simulations and experiments utilizing photoelasticity to
> visualize damage during failure. At the optimal level of disorder, the
> toughness is more than 2.6$\times$ of an ordered lattice of equivalent
> density.

Architected mechanical metamaterials leverage geometry to achieve outstanding
mechanical properties, including superior stiffness- and strength-to-weight
ratios (?, ?, ?, ?, ?, ?, ?, ?). Only recently has the fracture behavior of
architected materials been investigated (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?,
?), despite the importance of failure via crack propagation in many materials.
Architecture provides an opportunity to improve mechanical performance by
controlling and enhancing fracture toughness, but this has not been fully
realized, as previous investigations have primarily focused on materials with
simple, periodic microstructures (?, ?, ?).

Seminal works on the fracture of lattice materials (?, ?, ?, ?) showed that
the toughness, $G_{c}$, is primarily a function of the connectivity
(controlling if the lattice is stretch- or bending-dominated), the relative
density $\rho^{n}$ ($n=1$ if stretch-dominated and $n>1$ if bending-
dominated), the unit cell size, $L$ (?), and the orientation of the ligaments.
Thus, for a fixed relative density, stretch-dominated lattices with large unit
cells achieve the highest toughness, with the only other meaningful variable
being the orientation of the ligaments, the effects of which are not well
understood.

Although often inspired by nature (?, ?, ?, ?), the structural complexity of
architected materials is comparatively limited, predominantly utilizing
regular, periodic patterns (?, ?, ?, ?, ?). This simplicity contrasts sharply
with the disordered and non-periodic structures found in natural materials
with high fracture toughness, such as bone (?) and nacre (?). Natural
structures are also often confined or a component in a larger system (?, ?),
unlike many of the freestanding structures that have been explored to date (?,
?). When more complicated or disordered geometries are investigated (?),
typically a small number of representative geometries are considered, leaving
it unclear if the measured properties are a product of the disorder and thus
generalizable, or are unique to the specific geometry.

Disordered structures have not been systematically studied primarily because
traditional fracture mechanics analyses of architected materials assume a
self-similar crack front during propagation (?, ?, ?), limiting the analysis
to periodic structures typically at small length scales. However, these same
analyses show that toughness is enhanced by increasing the unit cell size.
Previous works (?, ?, ?) on fracture of ordered architected materials posited
that the largest structures that are expected to significantly influence the
toughness have a comparable length scale to the process zone, $\lambda_{0}$
(?), the region of inelastic damage resulting from the singular crack-tip
stresses. In an elastic-brittle lattice, the stresses are not truly singular,
as the crack-tip is not sharp, nor is there a traditional process zone arising
from inelastic effects. However, there remains a region of high stress near
the crack tip and a length scale over which these stresses decay that is akin
to an effective process zone size (?). Features larger than the decay length
are unlikely to impact the toughness as the local fracture stresses will fully
develop within the feature itself.

While structures comparable to this length scale have already demonstrated
significant enhancements in toughness (?, ?, ?), they remain periodic and
assume an a priori crack path and a self-similar crack front. Identifying the
crack path in complicated disordered structures is challenging with recent
works using data-driven approaches, such as graph neural networks (?), to
predict crack path and toughness. These approaches, however, do not provide
insight into the mechanisms through which disorder affects toughness.

We investigate the fracture of 2-D plane strain triangular lattices with unit
cells comparable to the effective process zone in size, and quantify how
geometric disorder relates to the distribution of damage during crack
propagation and the toughness. We introduce a framework utilizing an effective
crack length that demonstrates how toughness may be estimated by identifying
the number of ligament failures per crack advance and the strengths of the
lattice unit cells. Results are verified through a combination of finite
element simulations and fracture experiments, which reveal that the
introduction of disorder changes the damage from a straight, continuous crack
path to discontinuous distributed damage, resulting in higher toughness. For a
confined lattice, we identify an optimal level of disorder where toughness is
maximized. These enhancements in toughness are achieved with minimal loss of
strength and negligible change in stiffness.

Figure 1: (A) A single unit cell showing the introduction of disorder by
perturbing node locations by $\delta$. Representative networks with varying
average levels of disorder (B) $\bar{\delta}=0$, (C) $\bar{\delta}=0.15$, and
(D) $\bar{\delta}=0.27$.

We consider a lattice of finite height (e.g., a lattice at an interface), as
shown in Fig. 1(B-D). Geometric disorder is incorporated while maintaining
consistent connectivity by perturbing the node locations in 2-D by distances
$\Delta x$ and $\Delta y$ (Fig. 1(A)). Each perturbation is selected from a
normal distribution. For a unit cell with length $L$, the disorder magnitude
is defined as $\delta=\sqrt{\Delta x^{2}+\Delta y^{2}}/L$. The network’s
average disorder, $\bar{\delta}$, is calculated as the mean of all the
perturbations, with larger $\bar{\delta}$ values corresponding to greater
disorder across the network, as shown in Figs. 1(C) and (D). The impacts of
different definitions of disorder are discussed below. Lattices with
perturbations that result in overlapping nodes are not considered, as they
would change the connectivity of the lattice, which is outside the scope of
this analysis. This places a geometric limit on the perturbation of one node
to be less than $\delta=0.433$. Additional practical constraints arising from
the finite thickness of the beams and neighboring nodes limit the average
level of disorder in a lattice to less than $\bar{\delta}\approx 0.375$.

## Effect of Disorder on Local Lattice Failure

The impact of disorder on the local strength of the lattice is assessed by a
unit cell analysis consisting of a single perturbed node and its connecting
ligaments, as shown in Fig. 1(A). As the lattice is stretch-dominated, only
the axial stresses in the ligaments are significant (?). Bending stresses and
stress concentrations at the nodes are neglected. The axial stress in each
ligament is $\sigma_{a_{i}}=F_{i}/A$, where $i$ denotes the ligament, $F_{i}$
is the force in the ligament, and $A=tb$ is the ligament’s cross-sectional
area, with $t$ being the in-plane width and $b$ the out-of-plane thickness. An
average normal stress $\bar{\sigma}_{yy}$ is applied to the upper and lower
boundaries of the unit cell, and equilibrium at each node requires that
$\sum\limits_{i\epsilon N_{j}}F_{i_{x}}=0$ and $\sum\limits_{i\epsilon
N_{j}}F_{i_{y}}=0$, at each node $j$, where $N_{j}$ are the set of ligaments
that connect to each node. The solutions of these equations determine all
$\sigma_{a_{i}}$. Assuming the ligaments are elastic-brittle with a failure
stress, $\sigma_{f}$, this model predicts the unit cell’s strength,
$\sigma_{m}$, for any node perturbation. This approach applies to stretch-
dominated lattices with narrow ligaments ($t=L/10$ here), but can be adapted
to other cases through more sophisticated failure criteria.

Figure 2: Representative finite element stress distributions showing the
maximum principal stresses $\sigma_{mps}$ relative to the average applied
stress, $\sigma_{appl}$, for unit cells with (A) $\delta=0$, (B)
$\delta=0.10$, (C) $\delta=0.20$, (D) $\delta=0.30$. (E) Relative unit cell
strength as a function of disorder. Result from the analytical model is also
indicated by the solid curve.

The predictions from the analytical model are verified using 2-D finite
element simulations of the unit cell with varying perturbations to the central
node (Fig. 2)(A-D) under uniaxial tension ($\bar{\sigma}_{yy}=\sigma_{appl}$).
The node is perturbed at disorder $\delta$ from 0 to 0.4 in increments of
0.05, with seven simulations at each level of disorder corresponding to
perturbing the node an angle from 0$\degree$ to 180$\degree$ in increments of
30$\degree$. Periodic boundary conditions were applied to the left and right
sides of the unit cell. The average maximum principal stress, $\sigma_{mps}$,
in the middle of each ligament was used to predict the failure.

The strengths of the unit cells as a function of disorder, relative to an
unperturbed unit cell are given in Fig. 2(E) for the simulations and the
analytical model. The agreement between the two is close, although the
analytical model slightly underestimates the strength. This discrepancy arises
primarily from the omission of periodic boundary conditions in the analytical
model. The strength of the unit cell monotonically decreases with $\delta$,
but the loss of strength remains $\leq 10\%$ on average at the highest level
of disorder. The reduction in strength is due to the loss of symmetry, with
the stress being transferred to at least one ligament that is oriented more
along the line of loading than the others, causing premature failure.

## Effect of Disorder on Toughness

The fracture behaviors of lattices with varying geometries were calculated
using 2-D linear elastic finite element simulations. Crack growth was
simulated by running the model multiple times, as described below. For a given
network geometry, the model was run once to calculate stresses in the network.
The ligament with the highest stress was identified and deleted. The updated
model was then run again to identify the next ligament that would fail. This
process was repeated for up to 50 ligament failures. As the applied load and
stresses are linearly related, the applied load at each ligament failure was
determined by scaling the ligament stress to the failure stress. Additional
details about the model geometry and simulations are provided in the Materials
and Methods section.

The load-displacement curves for the ordered lattice and ten distinct
disordered lattices, each with $\bar{\delta}\approx 0.15$, are shown in Fig.
3(A). While the damage curves for the disordered configurations exhibit
considerable variability, they predominantly maintain loads that exceed those
of the ordered lattice, indicating higher damage resistance. Crack propagation
and toughness were investigated for 420 unique geometries of varying levels of
disorder between $0\leq\bar{\delta}\leq 0.375$. This allows for the effect of
disorder to be quantified independent of the specific lattice geometry.

The crack path through the network, as shown for the representative geometries
in Fig. 3(B), is identified by tracking the location of ligament failures. In
the case of the ordered lattice, crack propagation follows a straight path,
akin to a crack in a homogeneous material. Conversely, in most disordered
lattices, the damage path is neither straight nor always continuous. We define
an effective crack tip by considering the stress distribution in the network.
As shown in Fig. 3(B), averaging the stress across the height of the network
at every $x$-location, results in a homogenized stress profile that is similar
to the stress distribution around a crack tip in a continuum material with a
process zone (?). The effective crack length, $a^{\ast}$, is determined by
identifying the location of the maximum tensile stress arising from the
singular-like region of the stress distribution (Fig. 3(B)). This methodology
for defining $a^{\ast}$ not only aligns with the physical behavior of cracks
but also simplifies the subsequent toughness calculation, since the effective
crack length monotonically increases during failure but does not require an
assumption of self-similarity of the crack front. The stress distribution also
provides a measurement of the effective process zone size. We define the
effective process zone size, $\lambda^{\ast}$, as the full-width at half-
maximum of the stress distribution around the effective crack length. For the
ordered lattice, $\lambda^{\ast}\approx 2L$, which confirms that the unit cell
is comparable in size to the process zone.

Figure 3: (A) Representative load-displacement behavior for the ordered
lattice and ten geometries with $\bar{\delta}=0.15$. (B) Representative
ordered and disordered geometries, with ligament failure indicated. The
average stress profiles along the crack plane are shown with the effective
crack tip at the peak stress indicated along with the process zone,
$\lambda^{\ast}$. (C) Toughness (ave. $\pm$ std. dev. for all ligament
ruptures) of each simulated geometry, relative to the ordered lattice. Results
plotted versus eq. (3). A linear fit ($R^{2}=0.97$) is also indicated

The toughness is predicted by considering the strain energy change during
crack propagation. Most of the strain energy change is localized to the unit
cells around the ligaments that rupture, giving the total change in strain
energy, $\Delta U_{el}$, released during crack propagation as

$\Delta U_{el}=\frac{\eta}{2}\int_{V}\epsilon\sigma_{m}dV,$ (1)

where $V$ is the volume of the ligaments in the unit cell, $\sigma_{m}$ is the
stress in the unit cell at failure, $\epsilon$ is the strain, and $\eta$ is
the number of ligaments that fail. Since the network is stretch-dominated,
stress and strain are taken as approximately uniform through the volume of
each ligament, with $\epsilon_{m}=\sigma_{m}/E$, where $E$ is the Young’s
modulus. Thus, the critical energy release rate of the 2-D lattice is

$G_{c}=-\frac{1}{b}\frac{\Delta U}{\Delta
a^{\ast}}\approx-\frac{1}{b}\left(\frac{\bar{\sigma}_{m}^{2}}{E}\right)V\eta^{{}^{\prime}},$
(2)

where $\eta^{{}^{\prime}}=\Delta\eta/\Delta a^{\ast}$, and $\bar{\sigma}_{m}$
is the average stress in the unit cell at failure. Thus, the relative
toughness of a disordered lattice, compared to the ordered lattice toughness,
$G_{o}$, is

$\frac{G_{c}}{G_{o}}=\left(\frac{\bar{\sigma}_{m}}{\bar{\sigma}_{0}}\right)^{2}\left(\frac{\eta^{{}^{\prime}}}{\eta_{o}^{{}^{\prime}}}\right),$
(3)

where $\bar{\sigma}_{o}$ is the strength of an unperturbed unit cell,
$\eta_{o}^{{}^{\prime}}$ is the number of ligament failures per effective
crack propagation for the ordered lattice, and any small volume changes of the
ligaments between the ordered and disordered lattices are neglected.

The first term on the right-hand side of eq. (3) is $\leq 1$ for all $\delta$
(see Fig. 2(E)). Thus, any enhancements in toughness must come from sufficient
distributed damage to account for the loss in local strength. The toughness of
all simulated geometries (mean $\pm$ std. dev.), relative to the ordered
lattice, are shown in Fig. 3(C) as a function of eq. (3). The results are
well-described by the relationship and significant enhancements in toughness
(up to more than $3\times$) are achieved as a result of the distributed
damage. A linear fit to the data in Fig. 3(C) ($R^{2}=0.97$) has a slope of
1.1, slightly higher than the predicted slope of 1 due to a number of outlier
cases that are high and none that are low. This is expected based on eq. (3),
as there is no mechanism to generate significantly lower toughness than the
ordered lattice. The number of ligament failures cannot decrease relative to
the ordered lattice, and the losses in local strength are small ($\leq 10\%$),
resulting in a toughness that is typically equal to, or greater than, the
ordered lattice.

## Experimental Results & Discussion

Fracture experiments were performed on laser-cut polymethylmethacrylate (PMMA)
specimens as detailed in the Materials and Methods section (?). Three sets of
specimens were prepared: an ordered lattice, and two sets of disordered
structures. The disordered structures were generated by selecting two
geometries randomly from the simulated lattices, and then scaling the average
perturbation of each geometry to achieve average levels of disorder in the
range $0\leq\bar{\delta}\leq 0.34$. This allowed the effect of the disorder
magnitude to be investigated without introducing any bias arising from
specific geometries. A total of 33 specimens were tested, which included 3
ordered specimens, and 15 specimens in each of the disordered sets.

During fracture testing, specimens were analyzed using photoelasticity. For a
photoelastic material such as PMMA (?), the local stresses in the material
alter the propagation of light through the specimen. When imaged through cross
polarizers, this effect is visible such that the intensity of the observed
light is directly related to the difference in principal stresses.

Figure 4: Photoelastic images of (A) an ordered lattice, and (B) a disordered
lattice ($\bar{\delta}=0.27$), with ligament failures indicated. (C)
Representative load-displacement behavior for experimental specimens at
varying levels of disorder.

This technique allows for the identification of failing ligaments by
monitoring changes in local light intensity during fracture. Instances of
ligament failure are shown in Fig. 4(A) for an ordered lattice and in (B) for
a disordered lattice, where each specimen has the same length of crack
propagation. The observations reveal that disordered lattices exhibit more
distributed damage and a higher number of ligament failures during fracture
(33 for a disordered lattice with $\bar{\delta}=0.27$ vs. 27 for the ordered
lattice for cracks propagating over ${\sim}100$ mm).

Representative load-displacement curves for ordered and disordered lattices
are given in Fig. 4(C). Specimens with disordered structures sustained higher
loads during failure compared to the ordered lattice. However, specimens at
the highest levels of disorder ($\delta\geq 15\%$) showed no additional
enhancements in load capacity.

Figure 5: (A) Toughness, relative to the ordered lattice, as a function of
disorder percentage, for all finite element simulations (mean $\pm$ std. dev.
of $\sim 60$ geometries for each point), and experiments (mean $\pm$ std. dev.
for 3 identical specimens at each point).

Experimental toughness values as a function of the disorder are shown in Fig.
5, with each point being the mean and standard deviations for three specimens.
The toughness of the lattices was calculated from the load-displacement data
using a modified double cantilever beam analysis (?), as described in the
Materials and Methods section. Toughness predictions from the finite element
simulations are also shown, where each data point corresponds to ${\sim}60$
geometries (mean $\pm$ std. dev.). By considering many geometries at each
level of disorder, the toughness represents the average toughness for a system
where the number of unit cells is large, rather than simply a measurement of
any one unique geometry. The qualitative and quantitative agreement between
the simulations and experiments further verifies this result.

Consistent with the predictions of eq. (3), distributed damage enhances the
toughness of the structure compared to the ordered lattice, since
$\eta^{{}^{\prime}}/\eta_{o}^{{}^{\prime}}>1$. The effects of confinement
place an inherent constraint on the extent of the distributed damage and,
consequently, the maximum possible toughness enhancement. When the peak value
of $\eta^{{}^{\prime}}$ is reached but disorder increases further, a reduction
in toughness is expected due to diminished local strength. This is apparent in
both the experimental and simulation data, where the introduction of disorder
initially leads to a rapid increase in toughness as a result of the
distributed damage. The enhancement plateaus at $\bar{\delta}=0.15$ due to the
confined nature of the lattice studied here, and beyond $\bar{\delta}=0.2$ the
toughness begins to decrease as the local strength is reduced.

Experimentally, the maximum toughness achieved was 2.6$\times$ that of the
ordered lattice at $\bar{\delta}=0.15$. The value of the optimal level of
disorder is a function of the definition of disorder, which is not
standardized. Other methods for characterizing disorder, such as Voronoi
tessellation (?), will provide qualitatively identical results to those shown
here, but at varying quantitative levels of disorder.

## Conclusions

Disorder increases the fracture toughness of architected lattice materials by
generating distributed damage. Through finite element simulations of varying
lattice geometries, enhancements are shown to be common to the set of
geometries characterized by the level of disorder, not for any specific
architecture. These enhancements are achieved with minimal losses in strength
of $\leq 10\%$, relative to the ordered lattice, and negligible changes in
stiffness. For a confined lattice, a maximum enhancement in toughness is
identified due to the inherent limit of the extent of the distributed damage,
leading to an optimal level of disorder. The enhancements in toughness are
verified through experiments using photoelasticity to visualize the damage
during failure. The maximum increase in toughness was 2.6$\times$ at
$\bar{\delta}=0.15$. These enhancements in toughness demonstrate how
constraining the consideration of architected materials to periodic structures
significantly limits the possible mechanical performance.

A triangular lattice was considered in this work as a representative geometry,
but the approach can be extended to a broad set of lattices with varying
topology and connectivity. The approach is material-agnostic and is expected
to apply to most elastic-brittle materials. The framework presented in this
work can be applied to inform the failure analysis of a variety of disordered
structures, significantly expanding the design spectrum of architected
materials and enabling new enhancements in toughness.

## Materials and Methods

### Specimen Fabrication

Fracture specimens were fabricated from 6 mm-thick, cast, transparent
polymethylmethacrylate (PMMA) sheets, which were laser-cut using an x1250
laser cutter (Eduard, Inc., Denmark) with a 150-Watt laser operating at 80%
power and a cutting speed of 15 mm/min. A dual-pass laser cutting approach was
used to mitigate the effects of residual heat. Furthermore, to avoid
unintended failure of the solid beams and to isolate failure to the lattice
ligaments, rather than the nodes, an additional raster cut was introduced in
the middle of each ligament at a laser speed of 30 mm/min, resulting in a
consistent cut depth of ${\sim}2$ mm.

Figure 6: Optical image of laser-cut PMMA lattice (without raster cuts) with
ligament width and unit cell length indicated

Ligament thicknesses were measured optically and found to be $t=$ 0.55 mm
$\pm$ 0.05 mm, as shown in Fig. 6(A), and the unit cell length was measured to
be $L=$ 4.95 mm $\pm$ 0.05, both of which closely match the nominal value used
in the simulations. Specimens were annealed for 4 hours in a 90°C oven to
relieve residual stresses from fabrication. This temperature was chosen as it
is close to but below the glass transition temperature of PMMA (37).

Two sets of disordered lattices were fabricated using the seeds shown in Fig.
7. To generate disordered lattices, the nodes of the ordered lattice are
perturbed as shown. Lattices of varying magnitudes of disorder are generated
by changing the average magnitude of the perturbations, with the directions
and relative levels of perturbation being conserved. Structures at five levels
of disorder between $\bar{\delta}=0.07–0.34$ are generated, along with an
ordered lattice with $\bar{\delta}=0$ for a total of 11 unique geometries.
Three identical specimens of each geometry level were fabricated and tested
(33 total specimens).

Figure 7: Representative disordered lattices used for the experiments, with
$\bar{\delta}\approx 0.15$ in both cases. To generate the disordered lattices,
an ordered lattice is perturbed following the directions and relative
magnitudes shown for each seed. To generate lattices of varying levels of
disorder, the average magnitude of the perturbations is changed.

### Fracture Testing & Photoelasticity

Specimens were tested on a custom-built tensile testing instrument with dual-
head control. Specimens were tested under displacement control, with each head
moving at a rate of 2 mm/min, for a total displacement rate of 4 mm/min. Load
and displacement measurements were made at a rate of 100 Hz. Tests were
terminated after a crack had propagated through the entire specimen.

For each seed geometry and disorder level, one specimen was selected for
imaging via photoelasticity. The uniformity across specimens enabled the use
of a single specimen per group for representative analysis. The transmission
photoelastic setup was implemented as detailed in (38), incorporating circular
polarizers (Edmund Optics, NJ) and employing white light illumination. Imaging
was executed with a Point Grey camera, capturing data at a frequency of 2 Hz,
to visualize the stress distribution and crack propagation behavior. The
intensity of the observed light $I$, is known to be

$I\propto\sin{\left(\frac{\pi C_{\Lambda}b\Delta\sigma}{\Lambda}\right)}^{2},$
(4)

where $\Lambda$ is the wavelength of the light, $b$ is the specimen thickness,
$\Delta\sigma$ is the difference in principal stresses, and $C_{\Lambda}$ is a
material constant that relates the wavelength of the light to the
photoelasticity of the material (38).For a stretch-dominated lattice, the
axial stress in the ligament is dominant, thus
$\Delta\sigma\approx\sigma_{a}$. Thus, as ligaments fail, the intensity of the
observed light in each ligament will change according to the resulting stress
distribution. This change in light intensity is used to identify ligament
failures and the path of damage through the lattice during fracture.

Toughness was calculated from the load-displacement data using a modified
analysis for a double cantilever beam, with corrections for shear and a
compliant foundation using the Winkler foundation model. The crack length, $a$
was determined using the theoretical compliance (39) of the system, $C$, given
by

$C=\frac{8}{Eb}\left[\frac{3}{2}\alpha^{0.75}+\left(\frac{E}{4K\mu}+3\sqrt{\alpha}\right)\left(\frac{a}{h}\right)+\left(3\alpha^{0.25}+\frac{3}{2\pi}\sqrt{\frac{E}{\mu}}\right)\left(\frac{a}{h}\right)^{2}+\left(\frac{a}{h}\right)^{3}\right],$
(5)

with $\alpha=E/6E_{T}$, where $E$ is the Young’s modulus of the PMMA, $\mu$ is
the shear modulus, here taken as $\mu=E/2(1+\nu)$ for an isotropic material,
where $\nu\approx 0.3$ is the Poisson’s ratio, $E_{T}$ is the effective
tensile modulus of the lattice, $K$ is a constant equal to $5/6$, and $h$ is
the thickness of the beams.

In order to predict the effective crack length, the effective tensile modulus
of the foundation, $E_{T}$ must be known. For a stretch-dominated lattice, the
effective modulus is predicted (36) to be $E_{T}\approx\frac{1}{3}(\rho)E$,
where $\rho$, is the relative density of the lattice. For the lattices used in
this work where $t=L/10$, the effective modulus is predicted to be
$E_{T}\approx 0.129E$. This was found to be in good agreement with the
experimentally measured compliance.

The energy release rate of the lattice can therefore be calculated as a
function of crack length as

$G=\frac{P^{2}}{2b}\frac{\partial C}{\partial a},$ (6)

where $P$ is the load. For a homogeneous linear elastic material, the critical
energy release rate, $G_{c}$, is expected to be a constant with continuous
crack propagation. However, in the disordered lattices this is not guaranteed.
Variations in toughness result in cracks that arrest and propagate, as shown
by the jumps in the load-displacement curves in Fig. 4(C). The toughness of
the lattice corresponds to the critical energy release rate such that crack
advance occurs. Thus, the average toughness of the lattice was calculated from
all $G$ values during which the crack length increased as indicated by the
change in compliance of the specimen.

### Model Generation & Finite Element Simulations

Models were generated for finite element simulations using the open-source
computer-aided design software, FreeCAD, integrated with Python. Lattice
geometries were generated by creating a triangular lattice of nodes and
connecting them with rectangular ligaments. To introduce disorder to the
lattice, nodes were perturbed in $x$ and $y$ by selecting two random numbers
from a Gaussian distribution with zero mean. As the magnitude of the variance
in the distribution is increased, larger node perturbations occur. Using the
“random” command from SciPy, pseudo-random numbers are generated instead of
purely random numbers, with the exact number being reproducible if the same
“seed” value is input into the number generator. For this work, ${\sim}60$
seeds were generated, and scaled to 7 different levels of disorder (as shown
in Fig. 7), creating 420 unique lattices. The lattices consist of 30
horizontal nodes and 5 vertical nodes, with a unit cell size of $L=5$ mm and
ligament thickness $t=0.5$ mm. The overall length of the lattice is 150 mm,
and the thickness is 20 mm, with a portion of the top and bottom of the
lattice being embedded in the solid beams. The beams are 13 mm thick and 225
mm long, with the loading points 175 mm from the far edge of the lattice. A
representative model geometry with an ordered lattice is shown in Fig. 8, with
the boundary conditions indicated.

Figure 8: Double cantilever beam specimen geometry used for simulations and
experiments, with an ordered lattice. The lower loading point is fixed
horizontally and vertically, and the top loading point is displaced
vertically, causing a vertical reaction force, $F$.

Finite element simulations were performed in ABAQUS (2020, Providence, RI) by
importing the FreeCAD model file of the lattice. The lattice was modeled under
2-D plane strain conditions, with one mounting point in the beams being fixed
in $x$ and $y$, and the other under a 1 mm displacement boundary condition in
$y$ and fixed in $x$. The elastic modulus was taken as $E=1$ GPa, and the
failure stress was taken as $\sigma_{f}=20$ MPa. The linearity of the system
allows for the results to be easily scaled to materials of differing stiffness
and/or strength. The models were meshed with bilinear plane-strain
quadrilateral elements (CPE4), with a free mesh in the lattice and a uniform
mesh in the beams, with a mesh size of 0.05 mm. Approximately 110,000-130,000
elements were used per simulation, depending on the exact geometry. The mesh
size was refined until there was a less than 0.5% difference in the failure
force of an ordered and representative disordered lattice with a 10% increase
in mesh density.

To simulate crack propagation, the ligament under highest stress is identified
in the finite element model. The applied displacement and corresponding load
are scaled until the average von Mises stress in the ligament equals the
failure stress, $\sigma_{f}$. The von Mises stress was used to account for
possible nonuniformity in the stress distribution; however, as shown by the
unit cell analysis, the results are expected to be very similar to a maximum
principal stress criterion. The stresses are averaged in each ligament inside
an area centered in the ligament with a radius of $2t$. This ligament is then
deleted from the FreeCAD model and a new finite element simulation is run.
This process is run multiple times and up to 50 times per geometry. Since each
simulation in the sequence is treated as independent of the others, it is
possible that the predicted displacement that a ligament fails at is smaller
than the displacement at which the previous ligament failed. During a
continuous fracture test under displacement control, this is not possible.
Instead, these ligaments are expected to fail unstably during testing, and
these data points are not considered in the analysis.

## Acknowledgments

This research was funded primarily by the National Science Foundation (NSF)
MRSEC program [awards DMR-1720530 and DMR-2309043]. SF acknowledges support
from the US Department of Defense (DoD) through the National Defense Science &
Engineering Graduate (NDSEG) Fellowship Program. MKB acknowledges the Villum
Foundations for support under the Villum Experiment programme (VIL50302).

## References

  * 1. B. Compton, J. Lewis, Advanced Materials 26, 5930 (2014).
  * 2. T. Schaedler, et al., Science 334, 962 (2011).
  * 3. J. Mueller, J. Raney, K. Shea, J. Lewis, Advanced Materials 30 (2018).
  * 4. X. Zheng, et al., Science 344, 1373 (2014).
  * 5. J. Bauer, A. Schroer, R. Schwaiger, O. Kraft, Nature Materials 15, 438 (2016).
  * 6. L. Meza, S. Das, J. Greer, Science 345, 1322 (2014).
  * 7. S. Park, et al., ACS Nano 14, 8383 (2020).
  * 8. R. Ramachandramoorthy, et al., Applied Materials Today 20 (2020).
  * 9. S. Heide-Jorgensen, M. Budzik, K. Turner, Journal of the Mechanics and Physics of Solids 137 (2020).
  * 10. A. Shaikeea, H. Cui, M. O’Masta, X. Zheng, V. Deshpande, Nature Materials 21, 297 (2022).
  * 11. A. Mateos, W. Huang, Y.-W. Zhang, J. Greer, Advanced Functional Materials 29 (2019).
  * 12. A. Athanasiadis, M. Dias, M. Budzik, Extreme Mechanics Letters 48 (2021).
  * 13. M. Hossain, C.-J. Hsueh, B. Bourdin, K. Bhattacharya, Journal of the Mechanics and Physics of Solids 71, 15 (2014).
  * 14. C.-J. Hsueh, L. Avellar, B. Bourdin, G. Ravichandran, K. Bhattacharya, Journal of the Mechanics and Physics of Solids 120, 68 (2018).
  * 15. S. Luan, E. Chen, S. Gaitanaros, Journal of the Mechanics and Physics of Solids 169 (2022).
  * 16. M. O’Masta, L. Dong, L. St-Pierre, H. Wadley, V. Deshpande, Journal of the Mechanics and Physics of Solids 98, 271 (2017).
  * 17. L. Zhang, X. Mao, New Journal of Physics 20 (2018).
  * 18. Y. Cui, Y. Gao, H. Chew, International Journal of Solids and Structures 200-201 (2020).
  * 19. R. Muro-Barrios, Y. Cui, J. Lambros, H. Chew, Journal of the Mechanics and Physics of Solids 159 (2022).
  * 20. S. Fulco, M. Budzik, E. Bain, K. Turner, Extreme Mechanics Letters 57 (2022).
  * 21. S. Fulco, M. Budzik, K. Turner, Journal of the Mechanics and Physics of Solids 184 (2024).
  * 22. R. Hutchinson, N. Fleck, Journal of the Mechanics and Physics of Solids 54, 756 (2006).
  * 23. S. Maiti, M. Ashby, L. Gibson, Scripta Metallurgica 18, 213 (1984).
  * 24. N. Fleck, X. Qui, Journal of the Mechanics and Physics of Solids 55, 562 (2007).
  * 25. S. Luan, E. Chen, J. J., S. Gaitanaros, Science Advances 9 (2023).
  * 26. F. Barthelat, Z. Yin, M. Buehler, Nature Reviews: Materials 1 (2016).
  * 27. N. Glassmaker, A. Jagota, C.-Y. Hui, W. Noderer, M. Chaudhury, Proceedings of the National Academy of Science 104, 10786 (2007).
  * 28. M. Bacca, J. Booth, K. Turner, Journal of the Mechanics and Physics of Solids 96, 428 (2016).
  * 29. J. Chung, M. Chaudhury, Journal of The Royal Society Interface 2, 55 (2005).
  * 30. Q. Yang, B. Cox, R. Nalla, R. Ritchie, Biomaterials 27, 2095 (2006).
  * 31. Z. Huang, X. Li, Scientific Reports 3 (2013).
  * 32. S. Choukir, D. Egmond, B. Hatton, G. Hibbard, C. Sing, International Journal of Engineering Science 188, 1 (2023).
  * 33. M. Hedvard, M. Dias, M. Budzik, International Journal of Solids and Structures 288 (2024).
  * 34. G. Irwin, Journal of Applied Mechanics 24, 361 (1957).
  * 35. K. Karapiperis, D. Kochmann, Communications Engineering 2, 1 (2023).
  * 36. M. Ashby, Philosophical Transactions of The Royal Society: A 364, 15 (2005).
  * 37. K. Ute, N. Miyatake, K. Hatada, Polymer 36, 1415 (1995).
  * 38. K. Daniels, J. Kollmer, J. Puckett, Review of Scientific Instruments 88 (2017).
  * 39. R. Olsson, Composites Science and Technology 43, 329 (1992).
  * 40. H. Zhu, J. Hobdell, A. Windle, Journal of the Mechanics and Physics of Solids 49, 857 (2001).


11institutetext: Observatoire de Genève, Université de Genève, Chemin Pegasi
51, 1290 Versoix, Switzerland 22institutetext: CNRS, IRAP, 14 Avenue E. Belin,
31400 Toulouse, France 33institutetext: Cosmic Dawn Center (DAWN), Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200,
Denmark 44institutetext: Center for Astrophysics, Harvard & Smithsonian, 60
Garden Street, Cambridge, MA 02138, USA 55institutetext: Bogolyubov Institute
for Theoretical Physics, National Academy of Sciences of Ukraine, 14-b
Metrolohichna str., Kyiv, 03143, Ukraine

# First look with JWST spectroscopy: $z\sim 8$ galaxies resemble local
analogues

D. Schaerer 1122 R. Marques-Chaves 11 L. Barrufet 11 P. Oesch 1133 Y. I.
Izotov4 R. Naidu5 N. G. Guseva4 G. Brammer 33

(Accepted for publication in A&A Letters)

Deep images and near-IR spectra of galaxies in the field of the lensing
cluster SMACS J0723.3-7327 were recently taken in the Early Release
Observations program of JWST. Among these, two NIRSpec spectra of galaxies at
$z=7.7$ and one at $z=8.5$ were obtained, revealing for the first time rest-
frame optical emission line spectra of galaxies in the epoch of reionization,
including the detection of the important [O iii]$\lambda 4363$ auroral line
(see JWST PR 2022-035). We present an analysis of the emission line properties
of these galaxies, finding that these galaxies have a high excitation (as
indicated by high ratios of [O iii]$\lambda 5007$/[O ii] $\lambda$3727, [Ne
iii] $\lambda$3869/[O ii] $\lambda$3727), strong [O iii]$\lambda 4363$, high
equivalent widths, and other properties which are typical of low-metallicity
star-forming galaxies. Using the direct method we determine oxygen abundances
of $12+\log({\rm O/H})\approx 7.9$ in two $z=7.7$ galaxies, and a lower
metallicity of $12+\log({\rm O/H})\approx 7.4-7.5$ ($\sim 5$ % solar) in the
$z=8.5$ galaxy using different strong line methods. More accurate metallicity
determinations will require better data. With stellar masses estimated from
SED fits, we find that the three galaxies lie close to or below the $z\sim 2$
mass-metallicity relation. Overall, these first galaxy spectra at $z\sim 8$
show a strong resemblance of the emission lines properties of galaxies in the
epoch of reionization with those of relatively rare local analogues previously
studied from the SDSS. Clearly, the first JWST observations demonstrate
already the incredible power of spectroscopy to reveal properties of galaxies
in the early Universe.

###### Key Words.:

Galaxies: high-redshift – Galaxies: ISM – Cosmology: dark ages, reionization,
first stars

##  1 Introduction

Optical emission line spectroscopy has long provided important insights on the
physical composition, properties of the interstellar medium (ISM), and the
nature of the ionizing power of galaxies, yielding thus key information to
understand many key aspects of galaxy evolution (see review of Kewley et al.
2019). The well-known emission lines of H, He, O, N, S, Ne, and other elements
detected in optical galaxy spectra, emitted in the ionized ISM (Hii regions
primarily) have been detected in nearly one million of galaxy spectra, out to
redshifts of $z\sim 0.5-1$ with the Sloan Digital Sky Survey (SDSS, Ahumada et
al. 2020). Ground-based near-IR spectroscopy has recently pushed these limits
to $z\sim 1.5-3$, where $\sim 1500$ measurements of the strongest rest-optical
lines have been possible, e.g. with the MOSDEF survey (Kriek et al. 2015),
revealing thus ISM properties at cosmic noon (e.g. Förster Schreiber & Wuyts
2020).

The recent launch of the JWST and the spectroscopic capabilities of its multi-
object near-infrared spectrograph (NIRSpec) in the opens now a completely new
window into the early Universe, where, for the first time, all the “classical”
optical diagnostics developed at low-$z$ can be used to study the properties
of galaxies over a wide redshift range, from $z\sim 3$ out to the epoch of
reionization ($z>6.5$).

The first public NIRSpec observations, part of the Early Release Observations
(ERO) of JWST, have covered the SMACS J0723.3-7327 galaxy cluster, providing
$1.8-5.2$ $\mu$m spectra of 35 objects in the field. Among those, three
galaxies in the epoch of reionization were observed, showing spectacular, rich
rest-frame optical emission line spectra of objects at $z=7.7$ and $z=8.5$
(see JWST Press release 2022-035111https://webbtelescope.org/contents/news-
releases/2022/news-2022-035). We here report the first determination of the
metallicity (O/H) of these galaxies, a detailed analysis of their emission
line properties, and a comparison with observed properties of low-$z$ emission
line galaxies. At the time of revision, the JWST targets have been the object
of other studies analysing their emission line spectra and derived properties
(see Arellano-Córdova et al. 2022; Brinchmann 2022; Carnall et al. 2022; Curti
et al. 2022; Katz et al. 2022; Rhoads et al. 2022; Taylor et al. 2022; Trump
et al. 2022) .

##  2 Observations

Figure 1: 2D (top) and 1D (bottom) NIRSpec/JWST spectrum of 06355 at $z=7.664$
(black) and 1$\sigma$ uncertainty (grey). Vertical dashed lines mark the
position of well-detected nebular emission lines. The continuum emission is
also detected as seen in the 2D spectrum. X-axis in the bottom and top panels
refer to the observed ($\mu$m) and rest-frame wavelengths (Å), respectively.
Figure 2: Relative line intensities of the measured lines with respect to
H$\beta$. Full and open symbols show the ratios prior to and after applying
our flux correction (see text).

###  2.1 JWST NIRSpec observations

Rest-frame UV and optical spectra were obtained on 30 June 2022 using NIRSpec
with the micro-shutter assembly (MSA). Observations consist of two different
pointings (s007 and s008), each of them using two grating/filter combinations:
G235M/F170LP and G395M/F290LP. The total exposure time is 2$\times 8754$
seconds for each grating/filter. This provides a spectral resolution $R\simeq
1000$ and a continuous spectral coverage of $\simeq 1.75-5.20\mu$m. Fully
calibrated spectra (calibration level 3) were retrieved from the Mikulski
Archive for Space Telescopes222https://mast.stsci.edu/, which were previously
processed with the JWST Science Calibration Pipeline (cal_ver: 1.5.3 and
crds_ctx: jwst-0916.pmap). For each source, two individual 1D spectra (s007
and s008) are combined using the average flux, and after masking spectral
regions affected by cosmic rays and other artifacts, after inspection of the
2D spectra. For 04590 we exclude data above $\lambda_{\rm obs}\geq 4.5\mu$m
(i.e., after the detector gap) for the observation s008 as the 2D spectrum in
this region shows several artifacts that might compromise the flux
calibration, including a considerable flux gradient along the spatial
direction. An example of the combined spectrum of 06355 at $z\simeq 7.7$ is
shown in Fig. 1.

The spectra of the three sources show nebular emission lines, many of them
detected with high significance ($>3\sigma$). These include Balmer lines
(H$\beta$, H$\gamma$, H$\delta$), [O ii] and [O iii], [Ne iii], and He i. In
particular, the auroral [O iii]$\lambda 4363$ line – key for accurate
determinations of the O/H abundance using the direct method – is detected in
all galaxies with a significance of $3.9-5.7\sigma$. Gaussian profiles are
fitted to each line using the Python nonlinear least-squares function curve-
fit, and assuming a constant level for the continuum (in $f_{\lambda}$). We
derive the redshifts for these sources using the brightest lines, namely
H$\gamma$, H$\beta$, and [O iii] $\lambda\lambda$4959,5007 lines (Table 1).
The continuum is clearly detected in the two $z\simeq 7.7$ sources (10612 and
06355, see Fig.1) allowing the determination of the equivalent widths of the
lines (Table 1). However, these do not include any correction for slit-losses
and other possible effects due to different morphologies of the continuum
emission and nebular lines.

After careful inspection of the flux measurements we noticed that some line
ratios have nonphysical values, suggesting systematics on the flux calibration
and throughput in the current JWST pipeline. For example, the observed Balmer
line ratios are found to be larger than the intrinsic H$\gamma$/H$\beta$ and
H$\delta$/H$\beta$ ratios assuming case B recombination. Before more accurate
calibration reference files are available, we overcame this issue by applying
an ad hoc correction to the flux calibration. More specifically, we fit a
power law ($\propto\lambda^{\alpha}$) to i) the intrinsic
H$\gamma$/H$\beta=0.47$, H$\delta$/H$\beta=0.26$ (case B), and (He
i+H8)/H$\beta=0.2$ and to ii) the observed line ratios. The correction factor
as a function of wavelength is thus given by the division between i) and ii)
for each source, and it is applied to the full spectral range covered by the
G235M/F170LP grating/filter configuration. This empirical correction results
in a maximum increase of $\sim 20-50$% (depending on galaxy) for the O32=[O
iii] $\lambda$5007/[O ii] $\lambda$3727 ratio, and smaller corrections for
other line ratios. By doing this, we are also applying a first-order
correction to the Balmer decrement, thus losing the information on the dust
attenuation. Given the uncertainties on the flux calibration and on our
empirical correction, we conservatively add $20$% of uncertainties on our flux
measurements, which should also account for other sources of uncertainty, such
as those arising from the estimation of the continuum level. This depends on
the different assumptions (e.g., size of the spectral windows used in the fit,
constant level or straight line, noise fluctuations, etc.) and may affect the
flux measurements of different lines. The measured line ratios, prior to and
after correction, are shown in Fig. 2. Where possible, our analysis will focus
mostly on lines ratios between nearby lines (see Section 3), thus minimizing
possible issues on the flux calibration, and dust attenuation.

It must be noted that presently the data reduction is still in a preliminary
stage, which uses in particular pre-launch calibration files (Rigby et al.
2022). To tackle with these limitations, different groups have used different
approaches, including observations of a standard star, different extractions
of the 2D spectra, different combinations of the 2 visits and other
corrections (cf. Curti et al. 2022; Rhoads et al. 2022; Trump et al. 2022). A
comparison shows differences in several line ratios, by amounts which are
sometimes larger than the quoted uncertainties. To illustrate this, we
subsequently compare our measurements with those from Curti et al. (2022) who
have improved the calibration using a standard star. For other comparisons see
also Taylor et al. (2022) and Brinchmann (2022).

###  2.2 JWST NIRCam observations

We make use of a single NIRCam pointing in the six wide filters F090W, F150W,
F200W, F277W, F356W, and F444W, with a uniform exposure time of 2.1hr in each,
and shallower NIRISS imaging in F115W. We obtained the calibrated and
distortion-corrected NIRCam and NIRISS images from the publicly available
reduction by Brammer (2022). The images were processed with the standard JWST
pipeline up to stage 2b, before they were WCS aligned and combined with the
grizli333https://github.com/gbrammer/grizli/ pipeline (see Brammer et al., in
prep). The images were pixel aligned at 40mas pixel scales, before producing a
multi-wavelength catalog with SExtractor (Bertin & Arnouts 1996). Fluxes are
measured in small circular apertures of 0$\aas@@fstack{\prime\prime}$32
diameter and corrected to total fluxes using the AUTO fluxes from the F200W
detection image.

To determine stellar masses, SED fits to the 7 bands were done using the
latest versions of the CIGALE and Prospector codes (Boquien et al. 2019;
Johnson et al. 2021), exploring a relatively wide range of priors. Prior to
correction for magnification, we find masses of $\log(M_{\star}/M_{\odot})\sim
8.9-9.2$ and significant differences between the codes. We also note that the
masses derived by Carnall et al. (2022) are systematically lower than ours,
which should mainly be due to consistently younger ages found by these
authors. Before more detailed analyses of the SEDs and spectra of these
galaxies become available, we adopt the stellar masses from our preferred
CIGALE fits and a conservative uncertainty of $\pm 0.3$ (see Table 1). To
correct the masses for gravitational magnification from the cluster, we use
the magnification factors $\mu$ from the lens models of Caminha et al. (2022).
Other lens models, e.g. those used by Carnall et al. (2022), yield similar
magnifications (differences of 10-50% at most) for the sources studied here.

###  2.3 Comparison samples

For comparison with low-$z$ galaxies we use a sample of 5607 star-forming
galaxies from the SDSS Data Release 14 compiled by Y. Izotov and
collaborators, analysed in earlier publications (e.g. Guseva et al. 2019;
Ramambason et al. 2020). The selection criteria used for the extraction of
star-forming galaxies are presented in Izotov et al. (2014). Then we require a
detection of the [O iii]$\lambda 4363$ line with an accuracy better than
$4\sigma$, allowing thus direct abundance determinations using the
$T_{e}$-method. Subsequently we refer to this sample as Izotov-DR14.

We also use the spectra of 89 galaxies at $z\sim 0.3$ from the Low-Z Lyman
Continuum Survey (LzLCS), the first large sample of galaxies with UV
spectroscopy covering both the Lyman continuum and non-ionizing UV (see Flury
et al. 2022a, b). Approximately half of the sample has [O iii]$\lambda 4363$
detections.

Figure 3: Observed Ne3O2 ratio versus O32 (left panel) and [O iii]$\lambda
4363$/H$\gamma$ versus [O iii]$\lambda 5007$/[O iii]$\lambda 4363$ (right
panel) for the $z\sim 8$ galaxies (black stars: our measurements; green stars:
data from Curti et al. (2022)), the Izotov-DR14 sample of star-forming
galaxies (grey and red points, according to EW(H$\beta$)$<100$ Å and $>100$ Å,
respectively), and the LzLCS sample (blue points).

Figure 4: [Ne iii] $\lambda$3869/[O ii] $\lambda$3727 ratio (left panel) and
intensity of the auroral [O iii]$\lambda 4363$ line relative to [O
iii]$\lambda 5007$ (right panel), as a function of O/H. Same symbols as in
Fig. 3. Figure 5: Relation between the observed Ne3O2 ratio (left y-axis),
$\xi_{\rm ion}$ (right y-axis), and EW([O iii]$\lambda 5007$). The $z\sim 8$
galaxies are marked with by black stars for our measurements and green stars
using the line ratios and equivalent widths taken from Curti et al. (2022) and
Trussler et al. (2022) respectively. The shaded stripe shows the mean Ne3O2–EW
relation and scatter in the large low-$z$ sample of Izotov et al. (2021a). The
dotted lines show the relations between $\xi_{\rm ion}$ (right axis) and EW([O
iii]$\lambda 5007$) found at $z\sim 0-1$ (Chevallard et al. 2018; Izotov et
al. 2021a).  Figure 6: Mass-metallicity relation of the three $z\sim 8$
galaxies (black stars) compared to mean relations observed at $z=0$ and 2.2
(Andrews & Martini 2013; Sanders et al. 2020), and results from numerical
simulations at $z=2.2$ and 8 (dotted lines) from Ma et al. (2016). Green stars
show the results from adopting stellar masses determined by Curti et al.
(2022), who quote very small uncertaintes for their stellar masses.

##  3 Observed and derived properties of the $z\sim 8$ galaxies

###  3.1 Emission line properties

First we examine the observed emission line ratios in the three high-$z$
galaxies and compare them to those of the low-$z$ samples (cf. above). Figure
3 shows the main line ratios including the lines of [Ne iii] $\lambda$3869, [O
ii] $\lambda$3727, [O iii]$\lambda 4363$, [O iii]$\lambda 5007$, and
H$\gamma$, which are detected in the JWST spectra. The [Ne iii]/[O ii] (Ne3O2)
ratio is well known to closely trace O32, since both high ionization lines of
[Ne iii] and [O iii] originate in the same zone of the Hii region. Our
measurements in the high-$z$ galaxies (and those from Curti et al. 2022, shown
for comparison) are compatible with the observed correlation, providing
confidence for our empirical flux calibration. By these line ratios, the three
high-$z$ galaxies are found at relatively high excitation, corresponding to
the tail of the distribution in low-$z$ compact star-forming galaxies (cf.
Izotov et al. 2021a). Compared to the low-$z$ samples, the intensity of [O
iii]$\lambda 4363$ is also found to be relatively high, with [O iii]$\lambda
4363$$/{\rm H}\gamma\sim 0.3-0.5$, as illustrated in Fig. 3. The highest-$z$
source (04590) is a bit offset towards low [O iii]$\lambda 5007$/[O
iii]$\lambda 4363$.

The detection of the auroral [O iii]$\lambda 4363$ line and [O iii]$\lambda
5007$ provides access to the electron temperature $T_{e}$, and thus allows
abundance determinations using the so-called “direct method” (see e.g. Kewley
& Dopita 2002). To do this we follow the prescriptions of Izotov et al. (2006)
assuming low densities. The results are listed in Table 1. We find electron
temperatures $T_{e}({\rm O}^{2+})\sim 16000-18000$ K for the two $z=7.7$
galaxies, and nebular O/H abundances of $12+\log({\rm O/H})\approx 7.86-7.88$
for the two objects. This includes both ionic abundances of O+ and O2+,
determined from the optical lines; in both cases O2+ dominates. For 04590 the
unusually low [O iii]$\lambda 5007$/[O iii]$\lambda 4363$ ratio leads to
unphysically high electron temperatures ($T_{e}\sim 35$ kK) and in excess of
the highest accurate electron temperatures measured for star-forming galaxies
at low-$z$ ($T_{e}=24800\pm 900$ K, Izotov et al. 2021b). This result is
consistently found by all other papers who have analysed the JWST, and some
authors have derived electron temperatures as high as $T_{e}=37200\pm 9900$ K
(Rhoads et al. 2022). Although the detection of the [O iii]$\lambda 4363$ line
is beyond doubt, we consider that better data will be needed before the
electron temperatures can accurately be determined. Therefore the direct
metallicity determinations for these high-$z$ galaxies should presently be
taken with caution.

Different metallicity estimates can be obtained for the $z=8.495$ galaxy 04590
(see Table 1): From their finding of $T_{e}\ga 20$ kK, Trump et al. (2022)
conclude that $12+\log({\rm O/H})<7.69$ ($<1/10$ solar), from the empirical
relation between $T_{e}$ and $12+\log({\rm O/H})$ (cf. Pérez-Montero et al.
2021). Since standard strong line methods are not applicable at such low
metallicities and for objects with “extreme” line ratios, we can, e.g., use
the empirical calibration of Izotov et al. (2021b) based on measurements of
R23444$R23=($[O iii] $\lambda\lambda$4959,5007$/$H$\beta$$+$[O ii]
$\lambda$3727$/$H$\beta$$)$ and O32, and which has been established at
$12+\log({\rm O/H})<7.5$. For 04590 we find $12+\log({\rm O/H})=7.39$ with our
measurements, and $12+\log({\rm O/H})=7.24$ using the line ratios reported by
Curti et al. (2022), who derived $12+\log({\rm O/H})=6.99\pm 0.11$ by the
direct method. With the measurements from Rhoads et al. (2022) we obtain an
intermediate metallicity, $12+\log({\rm O/H})=7.30$, using the same
calibration. All studies so far agree that 04590 shows the lowest metallicity
among the three $z\sim 8$ galaxies.

In Fig. 4 we show the Ne3O2 and [O iii]$\lambda 4363$$/{\rm H}\gamma$ ratios
as a function of metallicity for the low-$z$ samples and the three $z\sim 8$
galaxies. We use these line ratios which are close in wavelengths to minimize
possible uncertainties of the flux calibration, differential slit losses and
others. Empirically, these line ratios can also provide a simple estimate of
the metallicity O/H, as discussed by earlier studies (see e.g. Nagao et al.
2006; Sanders et al. 2020). In any case, we see that our metallicity estimates
of the high-$z$ galaxies lead to fairly compatible locations in these
diagrams. We conclude that the three $z\sim 8$ star-forming galaxies have low
metallicities, in the range of $12+\log({\rm O/H})\sim 7.4-8.0$. More accurate
determinations and a proper evaluation of the uncertainties await overall
better data, including higher S/N data, proper calibrations, a more
sophisticated data reduction, and continuum and line flux extractions.

Since the continuum is also detected in the NIRSpec spectra of the two $z=7.7$
galaxies, and very weakly so also in the third object, we have also measured
the [O iii] equivalent width (see Table 1). We compare our measurements with
those from the low-$z$ galaxies in Fig. 5, where strong correlations between
EW([O iii]$\lambda 5007$) and properties such as O32, Ne3O2, and others have
been found (Tang et al. 2019; Izotov et al. 2021a). Possibly, the high-$z$
sources are somewhat offset. In any case, the EWs are high in the galaxies
with significant continuum detections, with EW([O iii]$\lambda 5007$)$\sim
700$ Å in the brightest source.

Empirically, the ionizing photon production efficiency, $\xi_{\rm ion}$, is
found to increase with the [O iii]$\lambda 5007$ equivalent width, as also
shown in Fig. 5. Using the relations found at low-$z$ and $z\sim 1-2$ (see
Tang et al. 2019; Izotov et al. 2021a), we estimate $\log(\xi_{\rm
ion})\approx 25.1-25.5$ erg-1 Hz, up to a factor $\sim 2$ higher than the
“canocical” value often assumed in ionizing photon budget calculations
(Robertson et al. 2013). This is comparable to other estimates of $\xi_{\rm
ion}$ at high redshift (e.g. Stefanon et al. 2022).

Table 1: Observed and derived quantities for the three high-$z$ galaxies

Quantity | 04590 | 06355 | 10612
---|---|---|---
redshift $z$ | 8.495 | 7.664 | 7.660
$T_{e}$(OIII) |  | 15584 | 18146
$12+\log({\rm O/H})$a |  | 7.86 | 7.88
$12+\log({\rm O/H})$b | 7.50 | 7.85 | 8.0
$12+\log({\rm O/H})$c | 7.39 | 7.76 | 7.91
O32 | 7.0 | 6.3 | 10.6
EW([O iii]$\lambda 5007$) | $172\pm 150$ | $723\pm 78$ | $638\pm 176$
$m_{F277W}$ | 27.8 | 26.6 | 25.6
$M_{1500}$d | $-20.29$ | $-21.09$ | $-20.38$
$M_{1500}$e | $-18.06$ | $-20.51$ | $-19.81$
$\beta_{1500}$ | $-2.20\pm 0.15$ | $-1.96\pm 0.22$ | $-2.31\pm 0.11$
$\log(M_{\star})\times\mu$ | $9.0\pm 0.3$ | $9.2\pm 0.3$ | $8.9\pm 0.3$
Magnification $\mu$ | 7.9 | 1.7 | 1.7
a direct method, b assuming $T_{e}=16000$ K
c using strong line methods (Izotov+2019,2021)
d observed, including lensing
e corrected for lensing

###  3.2 Physical properties of the $z\sim 8$ galaxies and the mass-
metallicity relation

As discussed, the three $z\sim 8$ galaxies show emission line properties
comparable to compact star-forming galaxies at low-$z$ with strong emission
lines. Izotov et al. (2021a) have shown that the low-$z$ galaxies with strong
lines (EW(H$\beta$)$>100$ Å) are good analogues of many of the $z\sim 1-3$
star-forming galaxies (Lyman alpha emitters and Lyman break galaxies) studied
so far. By inference the emission line properties of the three $z\sim 8$
galaxies studied here therefore also resemble those at intermediate redshifts.

By construction, the galaxies selected here cannot be claimed to be “typical”,
and larger, systematic studies will be needed. From our measurements, one
object (10612) shows a very high ratio [O iii]$\lambda 5007$$/{\rm
H}\beta\approx 10$ and a strong [O iii]$\lambda 4363$ line, which could
indicate an active galaxy (Seyfert 2, see also Brinchmann 2022; Curti et al.
2022). On the other hand, we clearly find evidence for one star-forming galaxy
with a fairly low metallicity (04590 with $12+\log({\rm O/H})\approx
7.2-7.4$).

If we combine our nebular metallicity estimates with the stellar masses
described earlier, we obtain the mass-metallicity relation shown in Fig. 6.
Our objects are found close to or below the mass-metallicity relation observed
at $z\sim 2$ and possibly offset by $\sim 0.2-0.3$ below this relation, in
good agreement with the relation derived by Ma et al. (2016) from simulations.
Similar results were also obtained by Jones et al. (2020) using an indirect
method based on ALMA emission line detections. However, at the present stage
we consider the stellar masses uncertain, since these depend significantly on
assumptions of the star formation history and on the age of stellar
populations (see also Tang et al. 2022), plus uncertainties with the
zeropoints of NIRCam (see Rigby et al. 2022). For example, comparing the
masses derived from SEDs by Carnall et al. (2022), Curti et al. (2022), and
our work, we can see significant differences which are generally larger than
the uncertainties cited. Adopting the lower masses from Carnall et al. (2022),
e.g., which would imply that our $z\sim 8$ galaxies lie close to the $z\sim 2$
mass-metallicity relation. The lower masses are mostly due to the very young
ages ($\sim 1-2$ Myr) inferred by these authors, whereas our SED fits favour
less extreme populations. Although Carnall et al. (2022) claim that the SEDs
show evidence for Balmer jumps (i.e. Balmer breaks in emission due to strong
nebular continuum), we do not see such a behaviour in the NIRSpec spectrum of
the brightest source, 06355, shown in Fig. 1.

Having shown that the $z\sim 8$ galaxies closely resemble strong emission line
galaxies from our low-$z$ sample it is tempting to infer indirectly other
properties using correlations established at low-$z$. Certainly, the high O32
ratios, low metallicity, and blue UV slopes ($\beta$) suggest that these
galaxies could contribute to cosmic reionization, i.e. have escaping ionizing
photons. For example, for the observed values of O32, $12+\log({\rm O/H})$,
and $\beta$, the LzLCS results suggest a 30-60% detection fraction of the
Lyman continuum. Adopting the mean relation between the LyC escape fraction,
$f_{\rm esc}$, and the UV slope, we estimate $f_{\rm esc}=0.03-0.08$ for the
three $z\sim 8$ galaxies, although the LyC escape could also be significantly
higher (see Chisholm et al. 2022).

Future improvements in the calibration and data reduction and additional
observations will allow us to determine more accurately equivalent widths and
total line fluxes of emission lines, measure continuum shapes, combine
photometry and spectra etc. and hence improve our knowledge of the physical
properties of galaxies at high-redshift.

##  4 Conclusion

We have analysed the rest-frame optical spectra of two galaxies at $z=7.7$ and
one at $z=8.5$ from the JWST Early Release Observations. The spectra, exhibit
numerous emission lines of H, He i, [O ii], [O iii], and [Ne iii], as commonly
seen in metal-poor star-forming galaxies at low redshift. They provide, for
the first time in the epoch of reionization, detailed information on the
chemical composition and interstellar medium of these galaxies. Our main
results are summarized as follows:

  * •

The auroral [O iii]$\lambda 4363$ line is significantly detected in all
galaxies with $3.9-5.7\sigma$, allowing the determination of the O/H abundance
(metallicity) using the direct method. With this method and different strong-
line methods we estimate metallicities between $12+\log({\rm O/H})=7.4-8.0$,
i.e. $\sim 5-20$ % of solar.

  * •

All three galaxies show a high excitation, as measured by their line ratios of
O32$=6-11$ and Ne3O2$=0.4-0.7$, or even higher, according to Curti et al.
(2022) and Rhoads et al. (2022). The observed emission line ratios are similar
to those of rare low-$z$ star-forming galaxies, which are considered analogues
of high-redshift ($z\sim 1-3$) galaxies. One of the $z=7.6$ galaxies shows
unusually high [O iii]$\lambda 5007$/H$\beta$$\approx 10$, possibly indicative
of nuclear activity (Seyfert 2).

  * •

The $z\sim 8$ galaxies show quite high equivalent widths, e.g. EW([O
iii]$\lambda 5007$) up to $\sim 700$ Å, as expected from low-$z$ galaxies with
high excitation. Such galaxies are known to be efficient producers of ionizing
photons. We conservatively estimate $\log(\xi_{\rm ion})\sim 25.2-25.5$ erg-1
Hz for our galaxies.

  * •

Using stellar mass estimates from SED fits, we find the $z\sim 8$ galaxies to
lie close to or below the mass-metallicity relation (MZR) at $z\sim 2$. To
assess if the MZR continues to evolve from $z=2$ to 8 will need future studies
of larger galaxy samples and accurate metallicity determinations.

Overall, the first analysis of the rest-frame optical spectra of galaxies at
$z\sim 8$ indicates that the emission lines properties of galaxies in the
epoch of reionization resemble those of relatively rare “local analogues”
previously studied from the SDSS, and for which numerous physical properties
have already been determined. These low-$z$ samples will soon be rivaled by
numerous measurements with NIRSpec onboard JWST. Clearly, the first data
release reveals already an extremely promising “preview” of upcoming science
with JWST in the early Universe. More robust inferences will require better
data, including higher S/N spectra and improvements in the calibration, data
reduction, and continuum and line flux extractions.

###### Acknowledgements.

We are very grateful to the ERO and other teams who made these JWST
observations possible. The Early Release Observations and associated materials
were developed, executed, and compiled by the ERO production team: Hannah
Braun, Claire Blome, Matthew Brown, Margaret Carruthers, Dan Coe, Joseph
DePasquale, Nestor Espinoza, Macarena Garcia Marin, Karl Gordon, Alaina Henry,
Leah Hustak, Andi James, Ann Jenkins, Anton Koekemoer, Stephanie LaMassa,
David Law, Alexandra Lockwood, Amaya Moro-Martin, Susan Mullally, Alyssa
Pagan, Dani Player, Klaus Pontoppidan, Charles Proffitt, Christine Pulliam,
Leah Ramsay, Swara Ravindranath, Neill Reid, Massimo Robberto, Elena Sabbi,
Leonardo Ubeda. The EROs were also made possible by the foundational efforts
and support from the JWST instruments, STScI planning and scheduling, and Data
Management teams. We acknowledge support from the Swiss National Science
Foundation through project grant 200020_207349 (PAO, LB). The Cosmic Dawn
Center (DAWN) is funded by the Danish National Research Foundation under grant
No. 140. Cloud-based data processing and file storage for this work is
provided by the AWS Cloud Credits for Research program. Y.I. and N.G.
acknowledge support from the National Academy of Sciences of Ukraine by its
priority project No. 0122U002259 “Fundamental properties of the matter and its
manifestation in micro world, astrophysics and cosmology”.

## References

  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3
  * Andrews & Martini (2013) Andrews, B. H. & Martini, P. 2013, ApJ, 765, 140
  * Arellano-Córdova et al. (2022) Arellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, arXiv e-prints, arXiv:2208.02562
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103
  * Brammer (2022) Brammer, G. 2022, Images and catalogs of HST and JWST images in the SMACS-0723 field
  * Brinchmann (2022) Brinchmann, J. 2022, arXiv e-prints, arXiv:2208.07467
  * Caminha et al. (2022) Caminha, G. B., Suyu, S. H., Mercurio, A., et al. 2022, arXiv e-prints, arXiv:2207.07567
  * Carnall et al. (2022) Carnall, A. C., Begley, R., McLeod, D. J., et al. 2022, arXiv e-prints, arXiv:2207.08778
  * Chevallard et al. (2018) Chevallard, J., Charlot, S., Senchyna, P., et al. 2018, MNRAS, 479, 3264
  * Chisholm et al. (2022) Chisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, arXiv e-prints, arXiv:2207.05771
  * Curti et al. (2022) Curti, M., D’Eugenio, F., Carniani, S., et al. 2022, arXiv e-prints, arXiv:2207.12375
  * Flury et al. (2022a) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022a, ApJS, 260, 1
  * Flury et al. (2022b) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022b, ApJ, 930, 126
  * Förster Schreiber & Wuyts (2020) Förster Schreiber, N. M. & Wuyts, S. 2020, ARA&A, 58, 661
  * Guseva et al. (2019) Guseva, N. G., Izotov, Y. I., Fricke, K. J., & Henkel, C. 2019, A&A, 624, A21
  * Izotov et al. (2014) Izotov, Y. I., Guseva, N. G., Fricke, K. J., & Henkel, C. 2014, A&A, 561, A33
  * Izotov et al. (2021a) Izotov, Y. I., Guseva, N. G., Fricke, K. J., et al. 2021a, A&A, 646, A138
  * Izotov et al. (2006) Izotov, Y. I., Stasińska, G., Meynet, G., Guseva, N. G., & Thuan, T. X. 2006, A&A, 448, 955
  * Izotov et al. (2021b) Izotov, Y. I., Thuan, T. X., & Guseva, N. G. 2021b, MNRAS, 504, 3996
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22
  * Jones et al. (2020) Jones, T., Sanders, R., Roberts-Borsani, G., et al. 2020, ApJ, 903, 150
  * Katz et al. (2022) Katz, H., Saxena, A., Cameron, A. J., et al. 2022, arXiv e-prints, arXiv:2207.13693
  * Kewley & Dopita (2002) Kewley, L. J. & Dopita, M. A. 2002, ApJS, 142, 35
  * Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511
  * Kriek et al. (2015) Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15
  * Ma et al. (2016) Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140
  * Nagao et al. (2006) Nagao, T., Maiolino, R., & Marconi, A. 2006, A&A, 459, 85
  * Pérez-Montero et al. (2021) Pérez-Montero, E., Amorín, R., Sánchez Almeida, J., et al. 2021, MNRAS, 504, 1237
  * Ramambason et al. (2020) Ramambason, L., Schaerer, D., Stasińska, G., et al. 2020, A&A, 644, A21
  * Rhoads et al. (2022) Rhoads, J. E., Wold, I. G. B., Harish, S., et al. 2022, arXiv e-prints, arXiv:2207.13020
  * Rigby et al. (2022) Rigby, J., Perrin, M., McElwain, M., et al. 2022, arXiv e-prints, arXiv:2207.05632
  * Robertson et al. (2013) Robertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, ApJ, 768, 71
  * Sanders et al. (2020) Sanders, R. L., Shapley, A. E., Reddy, N. A., et al. 2020, MNRAS, 491, 1427
  * Stefanon et al. (2022) Stefanon, M., Bouwens, R. J., Illingworth, G. D., et al. 2022, ApJ, 935, 94
  * Tang et al. (2019) Tang, M., Stark, D. P., Chevallard, J., & Charlot, S. 2019, MNRAS, 489, 2572
  * Tang et al. (2022) Tang, M., Stark, D. P., & Ellis, R. S. 2022, MNRAS, 513, 5211
  * Taylor et al. (2022) Taylor, A. J., Barger, A. J., & Cowie, L. L. 2022, arXiv e-prints, arXiv:2208.06418
  * Trump et al. (2022) Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2022, arXiv e-prints, arXiv:2207.12388
  * Trussler et al. (2022) Trussler, J. A. A., Adams, N. J., Conselice, C. J., et al. 2022, arXiv e-prints, arXiv:2207.14265