# Throughput Maximization in Multi-Band Optical Networks with Column
Generation

Cao Chen1, Shilin Xiao1, Fen Zhou2, Massimo Tornatore3  1State Key Laboratory
of Advanced Optical Communication Systems and Networks, Shanghai Jiao Tong
University, China  2 CERI-LIA, University of Avignon, France  3Department of
Electronics, Information and Bioengineering, Politecnico di Milano
Emails:{chen.cao, slxiao}{at}sjtu.edu.cn, fen.zhou{at}univ-avignon.fr,
massimo.tornatore{at}polimi.it

###### Abstract

Multi-band transmission is a promising technical direction for spectrum and
capacity expansion of existing optical networks. Due to the increase in the
number of usable wavelengths in multi-band optical networks, the complexity of
resource allocation problems becomes a major concern. Moreover, the
transmission performance, spectrum width, and cost constraint across optical
bands may be heterogeneous. Assuming a worst-case transmission margin in U, L,
and C-bands, this paper investigates the problem of throughput maximization in
multi-band optical networks, including the optimization of route, wavelength,
and band assignment. We propose a low-complexity decomposition approach based
on Column Generation (CG) to address the scalability issue faced by
traditional methodologies. We numerically compare the results obtained by our
CG-based approach to an integer linear programming model, confirming the near-
optimal network throughput. Our results also demonstrate the scalability of
the CG-based approach when the number of wavelengths increases, with the
computation time in the magnitude order of $10\text{\,}\mathrm{s}$ for cases
varying from 75 to 1200 wavelength channels per link in a 14-node network.

###### Index Terms:

Multi-Band, Flexible Optical Networks, Column Generation, Throughput
Maximization.

##  I Introduction

In the past decade, traditional wavelength division multiplexing (WDM) optical
networks characterized by fixed-grid frequency slots of 50 GHz have moved
toward flexible optical networks that benefit from a higher flexibility in
assigning fiber spectrum resources, and which can now allocate variable-size
frequency slots at granularities multiples of 12.5 GHz. More recently, multi-
band (MB) transmission was proposed to exploit fiber spectrum resources beyond
the C-band, e.g., L and U-bands[1]. For these bands, the spectrum width,
transmission performance, and cost of devices are very different. These
differences across different bands shall be carefully considered when
performing the Route and Wavelength Assignment (RWA). Moreover, due to the
increasing number of usable wavelengths in MB optical networks, the RWA
complexity becomes a new matter of concern.

This paper investigates the problem of throughput maximization in MB networks
that requires to solve the Route, Wavelength, and Band Assignment (RWBA). We
formulate RWBA as an integer linear programming (ILP) model and address the
high complexity by proposing an approach based on column generation (CG). In
this approach, a column corresponds to a colorless wavelength configuration,
with multiple link-disjoint routes to be assigned to the same wavelength. A
column may be used multiple times across the whole spectrum. This formulation
had different names in prior studies regarding RWA, such as, among others,
independent set[2] and routing configuration[3]. We choose here wavelength
configuration[4] to highlight its strict relation to wavelength channel
assignments. After the first appearance of the concept of wavelength
configuration in CG in [3], new and improved formulations [5], enhancements to
the algorithm’s scalability [4], and expansion to multicast routing [6] were
studied. In this paper, we extend the concept of wavelength configuration in
CG to incorporate: (i) flexible transceivers that can modulate the
transmission parameters according to the quality of transmission (QoT) of
lightpath, and (ii) MB transmission, considering the heterogeneous
requirements across optical bands, and addressing its inherent scalability
issue.

##  II Network Model and Problem Statement

This paper studies the RWBA problem considering a fixed-grid optical network.
$G=(V,E)$ is used to denote an optical network with node set $V$ and link set
$E$. Length of a link $\ell\in E$ is characterized by $N_{\ell}$ spans. Each
link has a pair of fibers that support bidirectional traffic and an equal
amount of spectrum resources. We assume a spectrum width of
$15\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ per fiber including the U, L, and
C-bands, with a bandwidth of $5\text{\,}\mathrm{T}$Hz per band111This paper
assumes U+L+C bands due to the possible minimum fiber loss around its center
frequency at around 1570 nm [7]. It cannot be excluded that transmission
devices for S-band will appear before those for U-band.. A lightpath $p$
denotes a simple path from source $s$ to destination $d$ on a particular
wavelength. The assigned wavelength should be consistent across all traversed
links to maintain wavelength continuity. The carried transmission capacity is
a function of the path length and the transmission margin requirement on its
band (to be explained later).

Given a normalized traffic demand matrix $\hat{D}$,
$\sum_{s,d}\hat{D}_{s,d}=1$, we aim to maximize the network throughput, i.e.,
the total data assigned for the demands. We need to allocate transmission
capacities for each node pair using adequate lightpaths, subject to the
following constraints:

  * •

Wavelength conflict constraint: no two lightpaths can use the same wavelength
on the same link;

  * •

Wavelength limitation constraint: available wavelengths should not exceed the
spectrum resources. Note that more attention is required when modeling the
RWBA problem, since different bands may hold different spectrum widths;

  * •

Demand-capacity constraint: traffic demand between two nodes should not exceed
the total transmission capacity between them, so that we can assess the
network performance in a non-blocking state;

  * •

Cost constraint: we express the possible cost constraint from network
operators by limiting the number of transceivers in a network, or lightpath
count equivalently.

The transmission capacity of a lightpath varies with the modulation format
based on the signal-to-noise ratio (SNR). We consider the modulation formats
{PM-BPSK, PM-QPSK, PM-8QAM, PM-16QAM, PM-32QAM, PM-64QAM, PM-128QAM,
PM-256QAM}222These formats are with net spectral efficiencies of {1.6, 3.1,
4.7, 6.3, 7.8, 9.4, 10.9, 12.5} in
$\text{\,}\mathrm{bit}\text{\,}{\mathrm{s}}^{-1}\text{\,}\mathrm{Hz}$ and
minimum required SNRs of {3.7, 6.7, 10.8, 13.2, 16.2, 19, 21.8, 24.7} in dB
scale[8].. Once the SNR of a lightpath exceeds the minimum required SNR, any
modulation format from the candidate list can be chosen. The transmission
capacity is given as follows,

$\displaystyle C=S_{m}\cdot R_{s}$ (1)

where $S_{m}$ denotes the spectral efficiency for the $m^{\text{th}}$
modulation format, $R_{s}$ is the transceiver’s baud-rate.

For the QoT evaluation, we selected a physical layer model that can capture
the effect of band-dependent transmission impairments, while being simple
enough to permit the scalable QoT computation required for a networking study.
We adopt the approach in [9] that computes the nonlinear interference based on
the conventional incoherent Gaussian Noise model. This approach first computes
the effective attenuation coefficient $\alpha_{\text{eff},w}$ and effective
length $L_{\text{eff},w}$ for a wavelength channel $w$ on a span with fiber
attenuation coefficient $\alpha_{w}$ [7] and span length $L_{s}$. The optical
signals experience the inter-channel stimulated Raman scattering (ISRS) effect
that can be approximated by a linear function with slope $C_{R}$ and cutoff
frequency $f_{R}$. The power loss and tilt after each span are assumed to be
fully compensated. For a fully-loaded reference spectrum with a center
frequency $f_{c}$, a channel spacing $B_{ch}$, total bandwidth $B_{t}$,
wavelength set $\\{w\in\mathbb{Z}\left|1\leq w\leq W\right.\\}$, a uniform
power level $P_{0}$ per channel, and total power $P_{t}$, the SNR for the
$w^{\text{th}}$ channel can then be calculated as follows,

$\displaystyle SNR_{w}$
$\displaystyle=P_{0}/\left(P_{\text{ASE},w}+\eta_{w}P_{0}^{3}\right)$ (2)
$\displaystyle P_{\text{ASE},w}$
$\displaystyle=\exp\left(\alpha_{\text{eff},w}L_{s}\right)\cdot\text{NF}\cdot
R_{s}\cdot h(f_{c}+f_{w})$ (3) $\displaystyle\eta_{w}$
$\displaystyle\approx\frac{8}{27}\gamma_{w}^{2}\alpha_{\text{eff},w}L_{\text{eff},w}^{2}\frac{\textrm{asinh}\left(\frac{\pi^{2}|\beta_{2}|B_{t}^{2}}{2\alpha_{\text{eff},w}}\right)}{\pi|\beta_{2}|R_{s}^{2}}$
(4)

with noise figure NF, nonlinear coefficient $\gamma_{w}$[7], channel frequency
$f_{c}+f_{w}$, and group velocity dispersion $\beta_{2}$. We assume an ideal
rectangular Nyquist-WDM spectrum for each channel and the channel spacing
equal to the baud-rate, and ignore the impact of transceiver noise and
filtering effect.

Next, we describe the key assumption used to set the margins in the three
bands. We transform the wavelength-dependent SNR into band-dependent SNR by
choosing the minimum SNR on that optical band. As such, we evaluate the
transmission performance on each optical band under full-load worst-case
assumptions, regardless of the actual deployments of optical amplifiers and
lightpaths. Given the minimum SNR after the first span
$SNR_{\min}^{(1)}=\min_{b\in\mathcal{B}}SNR_{b}$, the SNR after $N$ spans is
given as follows333By applying the typical parameters of a single mode fiber
and power spectral density $14\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{/}\mathrm{GHz}$, we have $SNR_{\min}\approx 120$ in a linear scale and
$SNR_{w}=\frac{120}{M_{b}N}$, $\forall w\in\mathcal{W}_{b}$, where $M_{U}$,
$M_{L}$, and $M_{C}$ adopts -4.4, -4.1, and 0 dB, respectively.,

$\displaystyle
SNR_{b,$\text{\,}\mathrm{dB}$}^{(N)}=SNR_{\min,$\text{\,}\mathrm{dB}$}^{(1)}-10\log_{10}N-M_{b,$\text{\,}\mathrm{dB}$}$
(5)

where $M_{b}$ is a preserved margin difference between the worst SNR of all
bands and the worst SNR on that band. Based on (2), we plot the SNR
distribution for different channel frequencies in Fig. 1, where we also report
the margin values for the U and L bands in the RWBA. Instead, when solving the
RWA (i.e., when we do not differentiate the performance of the three bands),
the SNR margin is considered to be 0s assuming the lightpaths on all optical
bands work with a flat worst-case SNR of about $20.8\text{\,}\mathrm{dB}$, so
that we can neglect the band allocation while always ensuring successful
transmission.

Figure 1:  Illustration of the SNR on different channel frequencies $f_{w}$
after the first span assuming the full-ISRS (blue) and no ISRS (red),
respectively. The transmission margins $M_{U}$ and $M_{L}$ are obtained by
differentiating the worst-case SNR for that optical band (solid gray line,
RWBA) and for all bands (dashed gray line, RWA). Two states of ISRSs are
considered in order to fully express the network-occupied statuses in worst-
case.

##  III Integer Linear Programming Model

The following subsection III-A formulates the problem of RWBA with throughput
maximization as an ILP model. Note that the ILP is based on path formulation,
where we precalculate $K$ routes via Yen’s algorithm[10]. As such we can
easily pre-estimate the SNR of each lightpath via (5) and choose the maximum
capacity from the highest-order modulation format compatible with that SNR.
Then, in subsection III-B, we use an example on a 4-node network to clarify
some of the issues of the ILP, and point out a promising solution.

###  III-A ILP Model

Parameters

  * •

$A$: maximum number of available transceivers.

  * •

$w\in\\{1,2,\cdots,W\\}$: wavelength index.

  * •

$\hat{D}_{s,d}$: normalized traffic demand profile,
$\sum_{s,d}\hat{D}_{s,d}=1$.

  * •

$p_{s,d,k}$: a path ID of the $k^{\text{th}}$ shortest path between $(s,d)$.

  * •

$\beta_{s,d,k,\ell}$: equals 1 if $p_{s,d,k}$ uses link $\ell$, 0 otherwise.

  * •

$C_{s,d,k,w}$: maximum transmission capacity for $p_{s,d,k}$ on $w$ assuming
the highest order modulation formats.

Variables

  * •

$\delta_{s,d,k,w}$: equals 1 if the $p_{s,d,k}$ adopts the $w^{\text{th}}$
wavelength, 0 otherwise.

  * •

$D_{s,d}$: traffic demand between $(s,d)$.

  * •

$T_{s,d}$: transmission capacity between $(s,d)$.

  * •

$TH$: network throughput.

$\displaystyle\\!{\max\limits_{\delta_{s,d,k,w}\in\mathbb{B}}}$ $\displaystyle
TH\quad\bm{}$ s.t. $\displaystyle D_{s,d}=TH\cdot\hat{D}_{s,d}\quad\quad$
$\displaystyle\forall(s,d)$ (6a) $\displaystyle D_{s,d}\leq
T_{s,d}\quad\quad\quad\quad$ $\displaystyle\forall(s,d)$ (6b) $\displaystyle
T_{s,d}=\sum\nolimits_{k}\sum\nolimits_{w}\delta_{s,d,k,w}C_{s,d,k,w}$
$\displaystyle\forall(s,d)$ (6c)
$\displaystyle\sum\nolimits_{s}\sum\nolimits_{d\neq
s}\sum\nolimits_{k}\delta_{s,d,k,w}\beta_{s,d,k,\ell}\leq 1$
$\displaystyle\forall\ell,w$ (6d)
$\displaystyle\sum\nolimits_{s}\sum\nolimits_{d\neq
s}\sum\nolimits_{k}\sum\nolimits_{w}\delta_{s,d,k,w}\leq A$ (6e)

The objective is to maximize the network throughput. Constraints (6a) ensure
the traffic demands distribution is compliant with the given parameters.
Constraints (6b) are the demand-capacity constraints. Constraints (6c) compute
the transmission capacities of the lightpaths connecting $(s,d)$. Constraints
(6d) avoid wavelength clash for all wavelengths on each link. Constraint (6e)
is the cost constraint used for limiting the transceiver count.

###  III-B A Small Example

Figure 2:  Illustration of a 4-node network with the span length on each fiber
and with a table on the right illustrating the transmission capacities for $9$
paths on $8$ wavelengths, $w\in\\{1,...,8\\}$, or on $2$ bands,
$b\in\\{1,2\\}$. (a), (b), and (c) show the respective assignment based on
path formulation and wavelength configuration formulations,
$c\in\\{1,2,3,4\\}$.

For the above model, scalability becomes a concern in the case of a large
number of wavelengths. Let us consider the 4-node network in Fig. 2 and
discuss the variables required to optimize $3$ demands $(1,4)$, $(2,3)$, and
$(2,4)$ when given $3$ candidate routes and $8$ wavelengths. In the
previously-discussed path-based formulation, at least $3\times 3\times 8$
boolean variables $\delta_{s,d,k,w}$ should be used to compute the maximal
$TH$,

$\displaystyle\max_{\delta_{s,d,k,w}\in\mathbb{B}}\left\\{TH\left|TH\leq\min\left\\{\frac{T_{1,4}}{\hat{D}_{1,4}},\frac{T_{2,3}}{\hat{D}_{2,3}},\frac{T_{2,4}}{\hat{D}_{2,4}}\right\\},\text{(\ref{eq:
conAddi_inthroughput_DT})-(\ref{eq: con3_cost})}\right.\right\\}$ (7)

where $\mathbb{B}$ is {0, 1} set. A possible lightpath assignment can refer to
Fig. 2(a), where
$T_{1,4}=\delta_{1,4,3,1}C_{p_{14},1}+\delta_{1,4,3,2}C_{p_{14},2}+\cdots=8C_{p_{14}}+C_{p_{134}}+C_{p_{124}}$,
$T_{2,3}=7C_{p_{213}}+3C_{p_{243}}$, and $T_{2,4}=5C_{p_{24}}$. For this
formulation, the number of variables increases linearly with the number of
wavelengths per link. As the solution space grows exponentially with the
number of variables, excessive computing time might soon be needed to search
for the optimal solution.

Inspired by the route similarity of the wavelengths 1, 2, 3, and 4, we propose
to apply the concept of “wavelength configuration” to solve the throughput
maximization problem such that a given routing configuration on a wavelength
is captured as one single variable. The first graph in Fig. 2(b) shows this
configuration, marked in gray to express that it is currently colorless. Note
that such a similar idea has been used in prior studies [4, 3], but prior
studies only considered lightpaths with fixed capacity that can always achieve
successful transmission for arbitrary lengths. Instead, we herein study the
wavelength configurations enabled by flexible transceivers with a distance-
adaptive transmission capacity.

As an example, in Fig. 2(b) and Fig. 2(c), we show some possible
configurations and compare the maximal $TH$ between fixed and flexible
transmission capacity and between RWA and RWBA. The transmission capacity
$T_{s,d}$ is now given as follows,
$\left[T_{1,4},T_{2,3},T_{2,4}\right]=z_{1}[C_{p_{14}},C_{p_{213}},C_{p_{24}}]+z_{2}[C_{p_{14}},C_{p_{213}}+C_{p_{243}},0]+z_{3}[C_{p_{14}}+C_{p_{124}}+C_{p_{134}},0,0]$,
where $z_{c}$ denotes the number of times a configuration is used. For the
fixed transmission capacity, the maximal $TH$ of
$2400\text{\,}\mathrm{Gbit}\text{/}\mathrm{s}$ can be achieved with the
optimal solution $z^{*}=(8,0,0)$. For the flexible transmission capacity, the
maximal $TH$ of $3000\text{\,}\mathrm{Gbit}\text{/}\mathrm{s}$ is achieved
with $z^{*}=(4,3,1)$. The network throughput can be further increased if we
account for a lower SNR margin preserved for wavelengths 7 and 8 on the
$2^{\text{nd}}$ band, $b=2$. Fig. 2(c) shows the optimal wavelength
configurations, achieving the maximal $TH$ of
$3150\text{\,}\mathrm{Gbit}\text{/}\mathrm{s}$ with $z^{*}=(5,1,1,1)$. Note
that the optimization variables in the above cases have changed from boolean
$\delta_{s,d,k,w}$ to integer $z_{c}$, reducing the memory consumption from
$72$ to $3$ or $4$.

This example shows the advantage of wavelength configuration in reducing the
optimization variables, how flexible transmission capacities can be
incorporated within the wavelength configuration concept, and how wavelength
configuration effectively lends itself to application for design of MB optical
networks. However, the number of potential configurations can also grow
exponentially in real network scenarios. So, in the next section, to more
efficiently identify and generate high-quality configurations, we present our
proposed CG approach.

##  IV Column Generation Approach

We begin by introducing CG decomposition in RWA to understand the iterative
process for generating a new column, including the restricted master problem
(RMP) and pricing problem. We then turn our attention to the RWBA to deal with
different optical bands. The CG approaches for RWA and RWBA are presented in
subsections IV-A and IV-B, respectively. Algorithm 1 shows the outline of CG.
Note that in line 1, we initialize the wavelength configuration set by a
sequential loading algorithm based on $k$-Shortest-Path (SP) and First-Fit
(FF). Up to $|V|(|V|-1)$ configurations, or $W$ configurations if
$W\leq|V|(|V|-1)$, are filled. As $k$SP-FF is also used for benchmarking, we
provide the details in Sec. V-A.

Algorithm 1 Column Generation Approach

1: Initialize set $\Omega_{\text{init}}=\\{c_{1},c_{2},...\\}$

2:
$c_{\text{new}}\longleftarrow\emptyset,\Omega\longleftarrow\Omega_{\text{init}}$

3: Do

4: $\Omega\longleftarrow\Omega\cup\\{c_{\text{new}}\\}$

5: Solve RMP$(\Omega)$ with relaxing variables, $z_{c}\in\mathbb{R}_{\geq 0}$

6: Choose a column $c_{\text{new}}$ via pricing problem(s)

7: While Reduced cost in line 6 is positive

8: Solve RMP($\Omega$) with integer variables, $z_{c}\in\mathbb{Z}_{\geq 0}$

9: Assign a wavelength to each configuration

###  IV-A CG for RWA

Parameters

  * •

$c\in\Omega$ : wavelength configuration ID.

  * •

$C_{s,d,k}$ : transmission capacity of path $p_{s,d,k}$.

  * •

$\delta_{s,d,k,c}$ : boolean parameter, equals 1 if path $p_{s,d,k}$ is used
by the configuration $c$, 0 otherwise.

  * •

$a_{s,d,c}$ : number of transceivers for a node pair $(s,d)$ in the
configuration $c$, $a_{s,d,c}=\sum_{k}\delta_{s,d,k,c}$.

  * •

$T_{s,d,c}$ : total transmission capacities for node pair $(s,d)$ in
configuration $c$, $T_{s,d,c}=\sum_{k}C_{s,d,k}\cdot\delta_{s,d,k,c}$.

  * •

$a_{c}$ : number of transceivers in a configuration,
$a_{c}=\sum_{s}\sum_{d\neq s}\sum_{k}\delta_{s,d,k,c}$.

Variables

  * •

$z_{c}$: number of times configuration $c$ is used.

####  IV-A1 RMP

$\displaystyle\\!\\!\\!\\!\max_{z_{c}\geq 0}$ $\displaystyle~{}TH\quad\bm{}$
s.t. $\displaystyle\hat{D}_{s,d}\cdot
TH+\sum\nolimits_{c\in\Omega}(-T_{s,d,c})z_{c}\leq 0\,\forall(s,d)$
$\displaystyle[\sigma_{s,d}]$ (8a)
$\displaystyle\sum\nolimits_{c\in\Omega}a_{c}z_{c}\leq A$
$\displaystyle[\sigma_{A}]$ (8b)
$\displaystyle\sum\nolimits_{c\in\Omega}z_{c}\leq W$
$\displaystyle[\sigma_{W}]$ (8c)

where $\sigma_{s,d}$, $\sigma_{A}$, and $\sigma_{W}$ are the dual variables
for demand-capacity constraint (8a), cost constraint (8b), and wavelength
limitation constraint (8c).

####  IV-A2 Pricing Problem

According to the dual of RMP, we can derive the pricing problem as follows,

$\displaystyle\max\limits_{\delta_{s,d,k,c^{\prime}}\in\mathbb{B}}$
$\displaystyle\sum\nolimits_{s}\sum\nolimits_{d\neq
s}\sigma_{s,d}\left(+T_{s,d,c^{\prime}}\right)-a_{c^{\prime}}\sigma_{A}-\sigma_{W}$
s.t. $\displaystyle\sum\nolimits_{s}\sum\nolimits_{d\neq
s}\sum\nolimits_{k}\delta_{s,d,k,c^{\prime}}\beta_{s,d,k,l}\leq 1\quad\forall
l$ (9a) $\displaystyle
T_{s,d,c^{\prime}}=\sum\nolimits_{k}C_{s,d,k}\delta_{s,d,k,c^{\prime}}\quad\forall(s,d)$
(9b) $\displaystyle a_{c^{\prime}}=\sum\nolimits_{s}\sum\nolimits_{d\neq
s}\sum\nolimits_{k}\delta_{s,d,k,c^{\prime}}$ (9c)

where $c^{\prime}$ is a configuration under design, with multiple routes
denoted by $\delta_{s,d,k,c^{\prime}}$. Constraints (9a) ensure no wavelength
clash on this configuration. Constraints (9b) and (9c) calculate the variables
in the objective function. The configuration is deemed promising if the
objective can be greater than 0.

Regarding the first two terms of the objective, this pricing problem can also
be expressed as a maximum-weight independent set problem, where each node
denotes a path associated with weight $\sigma_{s,d}C_{s,d,k}-\sigma_{A}$.
Rather than dealing with the time-consuming ILP model, we propose an
alternative algorithm, which constructs the node set by sequentially selecting
a node based on the descending order of node weight, until the weight of the
next node is 0 or negative. Meanwhile, we skip those nodes that create
wavelength clash with the already chosen nodes.

Next, we describe the utilization of the columns. Recalling that the variable
denoting the lightpath usage must be an integer, we compute such solution by
rerunning the RMP and restricting the variable $z_{c}$ to be an integer.
Finally, we take a simple wavelength packing strategy that groups the same
wavelength configuration together, sorts them by the configuration times, and
assigns contiguous wavelengths.

###  IV-B CG for RWBA

For RWBA, a wavelength configuration has to specify an optical band so that
the specific requirements of different bands can be considered. This is
fulfilled by introducing two parameters $\psi_{c,b}$ and $C_{s,d,k,b}$ to note
whether the configuration $c$ belongs to optical band $b$ and the transmission
capacity matrix of each band.

Additional Parameters

  * •

$b\in\mathcal{B}$ : optical band ID, $\mathcal{B}$={U, L, C}.

  * •

$W_{b}$ : available wavelengths on the $b^{\text{th}}$ optical band.

  * •

$\psi_{c,b}$ : equals 1 if wavelength configuration $c$ is used for the
$b^{\text{th}}$ optical band, 0 otherwise.

  * •

$C_{s,d,k,b}$ : transmission capacity of $p_{s,d,k}$ on the $b^{\text{th}}$
band.

####  IV-B1 RMP

$\displaystyle\max_{z_{c}\geq 0}\quad$ $\displaystyle TH\quad\bm{\texttt{}}$
s.t. $\displaystyle(\text{\ref{eq: demand-capacity
constraint}}),(\text{\ref{eq: limited lightpath constraints}})$
$\displaystyle\sum\nolimits_{c}\psi_{c,b}z_{c}\leq W_{b}\quad\quad\forall
b\in\mathcal{B}$ $\displaystyle[\sigma_{W,b}]$ (10a)

Constraints (10a) restrict the wavelength usage on respective optical bands.

####  IV-B2 Pricing Problem

$\displaystyle\max_{\begin{smallmatrix}\delta_{s,d,k,c^{\prime}}&\in\mathbb{B}\\\
\psi_{c^{\prime},b}&\in\mathbb{B}\end{smallmatrix}}$
$\displaystyle\sum_{s}\sum_{d\neq
s}\sigma_{s,d}\left(+T_{s,d,c^{\prime}}\right)-a_{c^{\prime}}\sigma_{A}-\sum_{b}\psi_{c^{\prime},b}\sigma_{W,b}$
s.t. $\displaystyle\text{(\ref{eq: pp_nonoverlapping})},\text{(\ref{eq:
pp_def2})}$ $\displaystyle
T_{s,d,c^{\prime}}=\sum\nolimits_{k,b}C_{s,d,k,b}\,\psi_{c^{\prime},b}\,\delta_{s,d,k,c^{\prime}}\forall(s,d)$
(11a) $\displaystyle\sum\nolimits_{b}\psi_{c^{\prime},b}=1$ (11b)

We consider the following algorithm for the above model. Given that an optical
band is unique to a wavelength configuration, we decompose this pricing
problem into $|\mathcal{B}|$ subproblems, each assuming the configuration is
aligned to a specific optical band. We then apply the same algorithm used for
the model (9) and choose the column with the most positive reduced cost.

##  V Illustrative Numerical Results

We now evaluate the computational efficiency and network performance of ILP,
CG, and two benchmark algorithms. We choose RWA to evaluate the impact of a
large number of wavelengths in MB optical networks. Then, we compare the RWBA
against RWA in order to evaluate the impact of the flexible capacity and MB
transmission.

###  V-A Simulation Setup

We consider two real networks for Deutsche Telekom AG (DT), 9-node and 17-link
and 14-node and 23-link noted as DT9 and DT14, respectively[11] with the
physical parameter settings as shown in Sec. II and uniform traffic demand
distribution among the node pairs i.e., $\hat{D}_{s,d}=\frac{1}{|V|(|V|-1)}$.

####  V-A1 Number of Wavelengths

As our aim is to compare the computational efficiency for different numbers of
wavelengths, we vary the number of wavelength channels by choosing different
channel spacings depending on varying the baud of the deployed
transceivers444Although we could also achieve the same goal by multiplying
wavelengths on a broader spectrum beyond U+L+C bands, this approach requires
tuning the physical layer settings, including the margin and transmission
capacity matrix. This adaptation involves a more complex physical layer
modeling that can be considered as future work.. Hence, the number of
wavelength channels varies from 30 to 75, 100, 150, 300, 600, 1000, 1200, and
1500, based on transceivers operating at 500, 200, 150, 100, 50, 25, 15, 12.5,
and 10 GBaud, respectively, when considering a total spectrum of
$15\text{\,}\mathrm{THz}$.

####  V-A2 Benchmark Algorithms

We consider two sequential loading algorithms $k$SP-FF and FF-$k$SP with $10$
candidate routes. $k$SP-FF establishes a lightpath by searching all viable
wavelengths for a given route, and turns to the next possible route, whereas
FF-$k$SP searches all possible routes on a given wavelength and turns to the
next possible wavelength[11, 12]. For the two algorithms, we uniformly load
the network with a fixed-rate demand that equals the minimum positive
transmission capacity of all lightpaths, so that at most one lightpath is
established for a new demand. The network throughput and lightpaths are
recorded every time a demand is loaded. Blocking occurs if a demand cannot be
fulfilled by any lightpath. When the first blocking occurs, the algorithm
stops and returns the throughput.

For the ILP in (6), the relative integer solution gap is set to 5% and the
maximum simulation time is $1\text{\,}\mathrm{h}$. For the RMP with integer
constraint, these parameters change to 1% and $10\text{\,}\mathrm{s}$ in order
to attain optimal performance within a reasonable computation time. Note that,
the typical solving time for the relaxed RMP is around
$2\text{\,}\mathrm{m}\mathrm{s}\mathrm{/}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}$
for DT9 and
$50\text{\,}\mathrm{m}\mathrm{s}\mathrm{/}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}$
for DT14 by using Gurobi 9.0 API and Matlab R2017b on a computer with 8 GB RAM
and a 4-core Intel i5-6300HQ.

###  V-B Performance Comparison in RWA

Fig. 3 shows the computational time and network throughput as the number of
wavelength channels increases. The results of the two heuristic algorithms and
the CG are averaged over $1000$ independent simulations. The points at which
100, 25, and 12.5 GBaud transceiver is assumed are marked by diamond,
pentagon, and circle symbols, respectively.

####  V-B1 Computational Efficiency

In Figs. 3(a) and (b), the computation time for the two benchmark algorithms
and the ILP increases as the number of wavelengths increases. This is caused
by the algorithms’ strategy that carefully searches the exact route and
wavelength for each lightpath. This process requires more time as the number
of usable wavelengths increases. In addition, the ILP can not be established
on the DT14 network with more than $1000$ wavelengths or with a 15 GBaud
transceiver due to memory limitations.

The CG approach, unlike traditional methodologies, shows strong scalability as
the number of wavelengths increases, maintaining a stable time consumption in
the magnitude order of $10\text{\,}\mathrm{s}$. This approach focuses on the
parallel allocation of each configuration rather than on a specific lightpath,
thus not sensitive to the number of wavelengths. Counterintuitively, the
computation time slightly decreases as the number of wavelengths increases.
For CG, the time-consuming module mainly comes from the slow convergence in
lines 3-7 of Algorithm 1, and finding the integer solution in line 8. The
observed time saving comes from the latter. As the wavelength increases, the
integrality gap of the whole problem decreases. The small integrality gap
means the last RMP with integer constraint could be approximated as a linear
programming (LP) model that has a fractional number of wavelengths within
$15\text{\,}\mathrm{THz}$ so that the advantage of high computation efficiency
inherent to LP is obtained.

Figure 3:  Computation time (top row) and network throughput (bottom row) for
ILP, CG, FF-$k$SP, and $k$SP-FF algorithms, for DT9 and DT14 networks, and for
variable numbers of wavelengths in $15\text{\,}\mathrm{THz}$.

####  V-B2 Max. Throughput

Regarding the Figs. 3(c) and (d), all algorithms exhibit an increase in the
network throughput as the number of wavelengths increases. This is well
understood because a more fine-grained spectrum is allowed in the network.
Note that the CG approach easily outperforms the benchmark algorithms in most
cases and is even comparable to the ILP model. A small solution gap between
the ILP model and CG exists, meaning that the current CG approach is near-
optimal instead of optimal. We also note an exception in Fig. 3(d) at the
point with 30 wavelengths where the CG is worse than the others. Most columns
in such a case are used only once, meaning that the RMP in finding the integer
solution is back into a complex 0-1 knapsack problem, which needs more time
instead of the $10\text{\,}\mathrm{s}$ in our setting. Increasing the time,
altering the initial configurations, and introducing a branch-and-bound
algorithm can improve the performance. Nonetheless, this exception does not
represent a big issue for CG, as the typical number of wavelengths exceeds 80,
e.g., CBand transmission on a 50 GHz grid.

###  V-C Performance Comparison in RWBA

Next, we turn to RWBA and compare its network performance with RWA. We choose
the DT9 network and 25 GBaud transceiver as an example, and show the results
in Fig. 4. They are marked with filled circles. We observe that solving RWBA
allows higher throughput, as the actual worst-case SNR in the respective
optical bands is typically greater than the value assumed in the RWA case.
Besides, allowing a larger number of modulation formats increases the network
throughput by 67% for RWA and 125% for RWBA.

As network operators are always concerned about the cost, we also consider the
achievable throughput under a cost constraint represented by a limited number
of transceivers, 8000 as an example. These results are also reported with
unfilled squares in Fig. 4. Adding a cost constraint decreases the achievable
throughput, but the results also confirm the benefit of the adapted modulation
format and the heterogeneous MB transmission requirements.

Figure 4: Maximum throughput of the DT9 using a 25 GBaud transceiver with
variable modulation formats compared for the RWA (red) and RWBA (blue), and
for the unlimited (filled circles) and limited transceivers (unfilled
squares).

##  VI Conclusion

We have proposed and studied the network throughput maximization problem in MB
fixed-grid optical networks with distance-adaptive modulation formats,
including the optimization of routing, wavelength, and band assignment. An ILP
model and a low-complexity CG-based approach were presented. Beginning with
the traditional RWA problem, we have validated the scalable performance of CG
compared to the ILP model as the number of wavelengths increases. Illustrative
numerical results, tested on realistic network topologies with numbers of
wavelengths per link varying from 75 to 1200, show that the CG is able to find
the near-optimal solution in the magnitude of $10\text{\,}\mathrm{s}$. This
enables the CG approach to scale to a greater network, whereas the ILP model
is not always scalable. When turning to the novel RWBA problem, we observe a
larger network throughput compared to the results obtained by solving the RWA
problem, since the modulation formats can be adapted to suit the transmission
margin requirements of these optical bands. Results with adapted modulation
format show a network throughput increase of 67% and 125% in the RWA and the
RWBA problems for the Germany DT9 network, respectively. This highlights the
significant potential of combining the CG approach with flexible transceivers.

## References

  * [1] N. Sambo _et al._ , “Provisioning in multi-band optical networks,” _IEEE J. Lightw. Technol._ , vol. 38, no. 9, pp. 2598–2605, May 2020.
  * [2] R. Ramaswami _et al._ , “Routing and wavelength assignment in all-optical networks,” _IEEE/ACM Trans. Netw._ , vol. 3, no. 5, pp. 489–500, 1995.
  * [3] T. Lee, K. Lee, and S. Park, “Optimal routing and wavelength assignment in wdm ring networks,” _IEEE J. Sel. Area. Comm._ , vol. 18, no. 10, pp. 2146–2154, Oct. 2000.
  * [4] B. Jaumard _et al._ , “Efficient spectrum utilization in large scale RWA problems,” _IEEE/ACM Trans. Netw._ , vol. 25, no. 2, pp. 1263–1278, 2017.
  * [5] B. Jaumard, C. Meyer, and B. Thiongane, “On column generation formulations for the RWA problem,” _Discrete Applied Mathematics_ , vol. 157, no. 6, pp. 1291–1308, 2009.
  * [6] F. Zhou _et al._ , “Joint optimization for multicast provisioning in mixed-line-rate optical networks with a column generation approach,” _IEEE J. Lightw. Technol._ , vol. 36, no. 3, pp. 637–649, Feb. 2018.
  * [7] N. A. Shevchenko, S. Nallaperuma, and S. J. Savory, “Maximizing the information throughput of ultra-wideband fiber-optic communication systems,” _Opt. Express_ , vol. 30, no. 11, pp. 19 320–19 331, May 2022.
  * [8] S. J. Savory, R. J. Vincent, and D. J. Ives, “Design considerations for low-margin elastic optical networks in the nonlinear regime,” _IEEE/OSA J. Opt. Commun. Netw._ , vol. 11, no. 10, pp. C76–C85, Oct. 2019.
  * [9] D. Semrau _et al._ , “Achievable rate degradation of ultra-wideband coherent fiber communication systems due to stimulated raman scattering,” _Opt. Express_ , vol. 25, no. 12, pp. 13 024–13 034, Jun. 2017.
  * [10] J. Y. Yen, “Finding the k shortest loopless paths in a network,” _Manag. Sci._ , vol. 17, no. 11, pp. 712–716, Jul. 1971.
  * [11] R. J. Vincent, D. J. Ives, and S. J. Savory, “Scalable capacity estimation for nonlinear elastic all-optical core networks,” _IEEE J. Lightw. Technol._ , vol. 37, no. 21, pp. 5380–5391, Sep. 2019.
  * [12] C. Chen, F. Zhou, Y. Liu, and S. Xiao, “Throughput maximization leveraging just-enough SNR margin and channel spacing optimization,” _IEEE J. Lightw. Technol._ , vol. 40, no. 13, pp. 4078–4093, Jul. 2022.


# Spin–gap formation due to spin–Peierls instability in $\pi$–orbital–ordered
NaO2

Mizuki Miyajima1,§,∗, Fahmi Astuti2,3,¶,∗, Takahito Fukuda1, Masashi Kodani1,
Shinsuke Iida4, Shinichiro Asai4, Akira Matsuo4, Takatsugu Masuda4,6,7, Koichi
Kindo4,
Takumi Hasegawa5, Tatsuo C Kobayashi1, Takehito Nakano8, Isao Watanabe2,3, and
Takashi Kambe1,†<EMAIL_ADDRESS>1Department of Physics, Okayama
University, Okayama 700–8530, Japan  2Advanced Meson Science Laboratory, RIKEN
Nishina Center, Wako, Saitama 351–0198, Japan  3Department of Physics,
Hokkaido University, Sapporo 060–0808, Japan  4Institute for Solid State
Physics, The University of Tokyo, Kashiwa, Chiba 277–8581, Japan  5Graduate
School of Advanced Science and Engineering, Hiroshima University,
Higashi–Hiroshima, 739–8521, Japan  6 Institute of Materials Structure
Science, High Energy Accelerator Research Organization, Tsukuba, Ibaraki
305–0801, Japan  7 Trans–scale Quantum Science Institute, The University of
Tokyo, Tokyo 113–0033, Japan  8Institute of Quantum Beam Science, Ibaraki
University, Mito, Ibaraki 310–8512, Japan

###### Abstract

We have investigated the low–temperature magnetism of sodium superoxide
(NaO2), in which spin, orbital, and lattice degrees of freedom are closely
entangled. The magnetic susceptibility shows anomalies at $T_{1}=220$ K and
$T_{2}=190$ K, which correspond well to the structural phase transition
temperatures, and a sudden decrease below $T_{3}=34$ K. At 4.2 K, the
magnetization shows a clear stepwise anomaly around 30 T with a large
hysteresis. In addition, the muon spin relaxation experiments indicate no
magnetic phase transition down to $T=0.3$ K. The inelastic neutron scattering
spectrum exhibits magnetic excitation with a finite energy gap. These results
confirm that the ground state of NaO2 is a spin–singlet state. To understand
this ground state in NaO2, we performed Raman scattering experiments. All the
Raman–active libration modes expected for the marcasite phase below $T_{2}$
are observed. Furthermore, we find that several new peaks appear below
$T_{3}$. This directly evidences the low crystal symmetry, namely, the
presence of the phase transition at $T_{3}$. We conclude the singlet–ground
state of NaO2 due to the spin–Peierls instability.

###### pacs:

Valid PACS appear here

The entanglement of spins, orbitals, charge, and lattice degrees of freedom is
a fundamental issue in solid–state physics Tokura ; Goodenough ; Kugel ;
Kanamori . They yield a variety of phenomena such as superconductivity,
quantum liquids, multi–ferroics, and orbital liquids BaCuSbO ; Kimura ;
Ishihara . While these physics have been discussed mainly in $d$– and
$f$–electron systems, there has been very little discussion in $p$– or
$\pi$–electron systems fullerene ; tdae_kambe . The magnetism of alkali–metal
superoxide, $A$O2, originates from the unpaired $p$–electrons on O2 molecule.
In O2, as two unpaired electrons occupy two $\pi_{g}^{*}$–orbitals,
$\pi_{x}^{*}$ and $\pi_{y}^{*}$, whose spins are in parallel to each other, O2
is a magnetic molecule with the spin quantum number, $S$, of 1. In $A$O2, as
$A$–cation is usually fully ionized to form $+1$ state, the additional
electron on O2 occupies one of two half–filled $\pi_{g}^{*}$–orbitals. This
allows the O${}_{2}^{-}$ molecule to gain the orbital degrees of freedom. As
three electrons exist on two $\pi_{g}^{*}$s, the energy band for the
O${}_{2}^{-}$ state is quarter–filled. The ab initio band calculations for
$A$O2 show that the Fermi energy locates within the $\pi$–band KO2_kim ;
RbO2_Mott ; KO2_Solovyev . On the contrary, the experimental magnetic
properties suggest electron localization, which implies the importance of an
electronic correlation on the O2 molecule. Thus, $A$O2 should be considered as
a Mott insulator KO2_kim ; RbO2_Mott . Spontaneous Jahn–Teller distortion
should occur to lift the degeneracy of $\pi_{g}^{*}$–orbital, which leads to
the selection of the $\pi^{*}$–orbitals. The coherent arrangement of the O2
orbitals is expected to lead to a three–dimensional magnetic exchange
interaction between the spins. Therefore, in $A$O2, spins, orbitals, charge
and lattice degrees of freedom are strongly coupled KO2_Solovyev . $A$O2 is
expected to be a candidate material to exhibit such fascinating phenomena.

Recently, $A$O2 has attracted much attention for its magnetic quantum
phenomena at low temperatures CsO2_Riyadi ; CsO2_NMR ; CsO2_ESR ;
CsO2_Miyajima ; RbO2_Fahmi ; Fahmi_thesis ; Miyajima_thesis . CsO2 shows
one–dimensional (1D) antiferromagnetic (AF) behavior in the magnetic
susceptibility and the high–magnetic field magnetization CsO2_Riyadi ;
CsO2_Miyajima . It was suggested that a 1D chain should form as a result of
the $\pi$–orbital ordering, but the detailed low–temperature structure has not
been determined CsO2_Riyadi . NMR experiments showed a power–law dependence of
the spin–lattice relaxation function, suggesting the emergence of a
Tomonaga–Luttinger Liquid state in the 1D short range ordered phase CsO2_NMR .

Contrary to RbO2 and CsO2, NaO2 has a cubic (Space group; $Fm$3̄$m$) symmetry
at room–temperature, in which O2 has an orientational disorder AO2_XRD . With
decreasing temperature, NaO2 shows successive structural phase transitions at
$T_{1}=220$ and $T_{2}\sim 196$ K AO2_XRD . In the marcasite phase below
$T_{2}$, the degeneracy of the $\pi_{g}^{*}$ orbitals is considered to be
lifted due to the low local symmetry around O2. The magnetic susceptibility
shows a weak decrease below $T_{2}$ and, then, a sudden drop below
$T_{3}=30\sim 40$ K. These experimental findings remind us of a low
dimensionality of the spin system and of a magnetic phase change below
$T_{3}$. Theoretical calculations for the marcasite phase pointed to the
quasi–1D AF character along the $c$–axis and frustration of exchange
interactions between different sublattice spins NaO2_Solovyev . However,
because of no experimental inspection on an existence of magnetic phase
transition and a change in crystal symmetry, the magnetic ground state of NaO2
has not yet been clarified. In this Letter, we investigate low–temperature
magnetic and structural properties of the $\pi$–orbital system NaO2 in detail.
For this purpose, we performed magnetic susceptibility, high–field
magnetization, muon spin relaxation ($\mu$SR), x–ray diffraction (xrd),
inelastic neutron scattering (ins) and Raman scattering experiments using high
quality samples.

First, we define three phases AO2_XRD ; the phase I above $T_{1}$, the phase
II between $T_{1}$ and $T_{2}$, the phase III between $T_{2}$ and $T_{3}$.
NaO2 has a remarkable temperature dependence of magnetic susceptibility,
$\chi=M/B$, where $M$ and $B$ denote the magnetization and the magnetic field,
respectively. Figure 1 (b) shows the temperature dependence of $\chi(T)$ for a
powder sample using a cooling and heating protocol with a magnetic field of
0.1 T. $\chi(T)$ shows anomalies around $T_{1}$, $T_{2}$ and $T_{3}$, which is
consistent with the literature AO2_magnetism . The temperatures of $T_{1}$ and
$T_{2}$ correspond well to the structural phase transition temperatures
supple_1 . In the phase I, the $\chi(T)$ follows the Curie–Weiss law with a
negative Weiss constant of $\theta=-9.4$ K. Small $\theta$ indicates a weak AF
interaction, which should be due to the orientational disorder of O2. In the
phase II, the $\theta$ changes to $41.1$ K, suggesting ferromagnetic
correlations due to the orientational ordering of O2. The effective magnetic
moments above and below $T_{1}$ are estimated to be about 1.82 and 1.68
$\mu_{\rm{B}}$, respectively, which increase slightly from the value expected
from the spin only. This may be due to the orbital effect. Following a clear
hysteresis around $T_{2}$, which indicates a first order phase transition,
$\chi(T)$ shows a weak decrease with decreasing temperature. As the localized
spins on O2 are responsible for the magnetism of NaO2, the remarkable decrease
of $\chi(T)$ in the phase III is an important key to consider the
low–temperature magnetism NaO2_Mahanti ; NaO2_Solovyev . In other words, it
may be reasonable to think that it comes from low–dimensionality. To evaluate
$J$ from the $\chi(T)$, we used the Bonner–Fisher model BF  and a
two–dimensional model with weak inter–chain interaction Keith , as shown in
Fig. S3 in Supplemental Material. However, we could not reproduce the
experiments by these models using the temperature–independent $J$. The xrd
experiments showed that the thermal shrinkage coefficient of the $c$–axis was
$-5.16\times 10^{-4}$ (Å/K), which was the largest value among the principal
axes. Thus, the shrinkage of the nearest neighbor (NN) length between O2s,
which corresponds to the $c$–axis length, may be sufficient to make AF $J$
stronger as the temperature decreases.

Below $T_{3}$, $\chi(T)$ decreases rapidly with decreasing temperature with no
temperature hysteresis and, then, increases. The $\mu$SR experiment showed
that no magnetic ordering was found down to 0.3 K as described later.
Therefore, the decrease of $\chi(T)$ implies the appearance of a spin–gap in
the spin excitation spectrum. Note that the low–temperature Curie–tail is
strongly dependent on the sample batch, indicating that it is due to the
extrinsic spins. To evaluate the intrinsic temperature dependence of
$\chi(T)$, we use the equation: $\chi(T)=C_{0}/T+C/T\exp\left({-2\Delta/k_{\rm
B}T}\right)$ , where the $\Delta$ denotes the spin–gap and $k_{\rm B}$ is the
Boltzmann constant muSR_NaTiSiO ; SP_NMR . The first term is the Curie–tail
contribution, which can be subtracted from the data, and the second term is
responsible for the spins excited from the singlet state to the magnetic
excited states. As shown in Fig. 1 (c), the fitting by this equation with
$\Delta/k_{\rm B}=51.2$ K is good. The Curie tail allows us to estimate the
number of extrinsic spins with $S=1/2$ to be 0.014 mol in this sample.

Figure 1:  (a) Unit cell of the marcasite phase of NaO2. The right figure
shows the unit cell in the ($ac$)–plane. $a_{\rm III}$ and $c_{\rm III}$
denote the axes in the phase III while $a^{\prime}$ and $c^{\prime}$
correspond to that in the phase II. (b) Temperature dependence of magnetic
susceptibility, $\chi(T)$, at $B=0.1$ T in powder NaO2. $\chi(T)$ shows clear
anomaly at $T_{1}$, $T_{2}$ and $T_{3}$. The arrows indicate the hysteresis in
the cooling and heating protocols. (c) Enlarged figure around $T_{3}$. The
experimental data $\chi(T)$ (blue dot), the low–temperature Curie–tail
$C_{0}/T$ (blue line) and the subtracted data $\chi(T)-C_{0}/T$ (red dot) are
shown. $\chi(T)$ is fitted by the equation in the text. (d) High–magnetic
field magnetization at 4.2 K, where the vertical axis indicates the
magnetization per O2. At low–magnetic field region, the $M(B)$ curve can be
fitted by the Brillouin function, $B_{s}$. The experimental data $M(B)$ (blue
dot), the $B_{s}$ (orange dotted line) and the subtracted data $M(B)-B_{s}$
(red dot) are shown.

Figure 1 (d) shows the magnetization curve, $M(B)$, as a function of $B$ up to
60 T at $T=4.2$ K. The saturation magnetization is equivalent to $\sim 1$
$\mu_{\rm B}$. Note that the $M(B)$ experiment was performed using the same
sample as the $\chi(T)$ experiments. At $T=4.2$ K, the $M(B)$ shows a
nonlinear increase at low fields and an anomaly with a large hysteresis around
$\sim 30$ T. To evaluate the low–magnetic field part of the $M(B)$ curve, we
use a Brillouin function as shown in Fig. 1 (d). The fitting was good and the
number of paramagnetic spins was obtained to be 0.011 mol. This value is
consistent with that obtained from the Curie–tail in the $\chi(T)$ experiment,
allowing us to subtract the low–field Brillouin contribution from the $M(B)$
curve. Accordingly, $M(B)$ shows the magnetic field induced transition from
non–magnetic to magnetic state with the large hysteresis around $\sim 30$ T.

To find out any signatures of a magnetic ordering in NaO2, we performed
$\mu$SR experiments down to 0.3 K. Figure 2(a) shows corrected time spectra
measured at 0.3 K ($\ll T_{3}$), 5 K and 100 K ($\gg T_{3}$) in the zero–field
(ZF) condition. Neither the loss of the initial asymmetry at $t$ = 0 nor the
muon–spin precession was observed down to 0.3 K. These findings exclude the
presence of a long–range magnetic order down to 0.3 K.

Figure 2:  (a) Zero–field (ZF) $\mu$SR spectra of NaO2 at 5 K and 100 K. The
inset shows ZF $\mu$SR spectrum at 0.3 K. The solid lines are the fitted
curves using the equation in the text. (b) Temperature dependence of
$\lambda$. The solid line is the fitted curve using the gap–related analysis
function described in the text muSR_NaTiSiO ; SP_NMR . The inset shows the
temperature dependence of $\sigma$. (c) ins spectra after background
subtraction. Details are described in Supplemental Material. (d) ins profiles,
where the spectra between 0.5 Å-1 and 1.5 Å-1 are integrated. Temperature
dependence of the intensity at 9.4 meV is shown in the inset.

To focus on collecting some data points around $T_{3}$, we analyzed the
depolarization rate obtained in ZF by using the function: $A(t)=Ae^{-\lambda
t}e^{-(\sigma t)^{2}}+{\rm B.G.}$ Kubo_Toyabe . In here, the Gaussian term
arises from the muon–spin relaxation caused by randomly distributed internal
fields coming from surrounding nuclear dipoles. The exponential term ascribes
the effect of the fluctuating electronic moments around the muonMuon_Relax_1 ;
Muon_Relax_2 . The B.G. indicates background signals from muons which do not
stop in the sample but in a sample mounting plate. These B.G. signals were
subtracted from the raw signal as a constant term in order to achieve
corrected $\mu$SR time spectra.

The temperature dependence of $\lambda$ and $\sigma$ are shown in Fig. 2 (b).
Above $T_{3}$, we obtained $\lambda\sim 0.1$ and $\sigma\sim 0.22$ $\mu$s-1.
The result of $\lambda$ obtained above $T_{3}$ indicates that electronic spins
likely fluctuate beyond the $\mu$SR characteristic time window ($10^{-6}\sim
10^{-11}$ sec), resulting in the motional narrowing limit. The $\lambda$
increases below $T_{3}$ suggests that the muon is expected to sense the
formation of the spin–gap state in NaO2. One possible scenario to explain this
result is that the muon spin relaxes its polarization by the thermally
activated electronic spins across the spin–gap. The same behavior was observed
in other spin–gap systems muSR_NaTiSiO . Following this scenario, the spin–gap
is estimated from the temperature dependence of $\lambda$ below $T_{3}$ by
applying the following function:
$\lambda(T)=\lambda_{0}\left\\{1+C^{\prime}\exp(-2\Delta/k_{\rm{B}}T)\right\\}^{-1}$
muSR_NaTiSiO ; SP_NMR . Using this equation, we estimated the
$\Delta/k_{\rm{B}}$ to be $\sim$ 44.6 K, which was consistent with that
obtained from the magnetic susceptibility measurement.

Figure 2(c) shows the ins spectra of powder sample at 2.70 K, where the
background contribution was subtracted (see Supplemental Material). The
excitation around $Q\sim 1$ Å-1 is shown to have a finite energy gap. The
intensity decreases with $Q$, which is typical behavior of magnetic
scattering, in the range of $Q\lesssim 2$ Å-1. The enhanced intensity at
$Q\gtrsim 2$ Å-1 is from remnant phonon scattering of the sample cell made of
Aluminum. The first momentum of the dynamical structure factor in 1D AF spin
chain is proportional to $1-\sin Qd/Qd$, where $d$ is the distance between
spins Igor , leading to the pronounced intensity at $Q\sim 1$ Å-1.

To reveal the change in the intensity as a function of temperature, the
spectra integrated between 0.5 Å-1 and 1.5 Å-1 are shown in Fig. 2(d). At 2.70
K, the intensity starts to increase at $E\sim 4$ meV and has a maximum at
$E\sim 9$ meV. The temperature dependence of the intensity at the maximum
energy is shown in the inset. It shows no temperature dependence above
$T_{3}$, and gradually increases with decreasing temperature below $T_{3}$.
This result directly indicates that NaO2 has the magnetic excitation with an
excitation gap energy of 9 meV below $T_{3}$. As the $\mu$SR experiments
indicated no magnetic long–range ordering down to 0.3 K, this peak results not
from magnetic excitation in the magnetic long–range ordered phase but from
singlet–triplet excitation in the non–magnetic ground state.

We searched for structural dimerization of the O2 molecules as the cause of
the non–magnetic state by the xrd and neutron diffraction measurements, but
could not experience the direct evidence on the structural change below
$T_{3}$ supple_2 . Then, we perform Raman scattering experiments because of
high sensitivity to changes in crystal symmetry and/or molecular charge.
Figures 3 (a) and (b) show the temperature dependence of the Raman scattering
spectra in the stretching and the libration mode region of O2, respectively.
We will focus on the change of the Raman–active modes around $T_{3}$.

As the crystal symmetry of the phase III is determined as $Pnnm$ ($D_{2h}$),
the sets of Raman active stretching and libration modes are given by
$\Gamma_{S}=A_{g}+B_{1g}$ and $\Gamma_{L}=B_{1g}+B_{3g}+A_{g}+B_{2g}$,
respectively AO2_Raman ; AO2_opt . As shown in Fig. 3 (a), two peaks can be
clearly observed in the stretching mode region. We can assign that the peaks
at 1163 cm-1 and 1140 cm-1 originate from the in–phase and the out–of–phase
stretching modes, namely, $A_{g}$ and $B_{1g}$, respectively. The $B_{1g}$
mode is observed for the first time. In the libration mode region (see Fig. 3
(b)), two major peaks are observed around 150 cm-1 ($L_{3}$) and 240 cm-1
($L_{1}$). Moreover, very weak peaks around 130 cm-1 ($L_{4}$) and 190 cm-1
($L_{2}$) are observed for the first time. These libration modes can be
assigned from the intensity and energy, but a more detailed consideration will
be needed. Anyway, these results demonstrate that all Raman active modes for
the phase with the $D_{2h}$ symmetry can be successfully detected AO2_Raman .
Below $T_{3}$, while no change was observed in the stretching mode region, in
the libration mode region, new peaks at 173, 86 and 56 cm-1 gradually
appeared, which are represented as $P_{1}$, $P_{2}$ and $P_{3}$ in Fig. 3 (b),
respectively. Figure 3 (c) shows the temperature dependence of the $P_{1}$
peak and Fig. 3 (d) summarizes the temperature dependence of the peak
intensity observed in the libration mode region. While the $\Gamma_{L}$ modes
depended weakly on the temperatures, the $P_{1}$, $P_{2}$ and $P_{3}$ peak
intensities increased markedly below 30 K, following order parameter like
behavior. Note that no splittings of the stretching modes are found,
suggesting an absence of charge ordering on O2. Thus, this result clearly
indicates an existence of a phase transition around $T_{3}$. Because all
Raman–active $\Gamma_{L}$ modes are confirmed in the phase III, the
observation of the new peaks is direct evidence for the low crystal symmetry
below $T_{3}$.

Figure 3: Raman scattering results of NaO2. (a) and (b) Temperature
dependences of the stretching modes and the libration modes below 50 K,
respectively. The spectra are shifted along the vertical axis for clarity. The
asterisk peak is a line included in the laser source. (c) Enlarged figure
around 150 cm-1 for the libration modes. (d) Temperature dependences of peak
intensity for the libration mode region, where the peak intensities are
normalized by the intensity at 5 K. Solid lines are guides to the eye.

To discuss the cause of the non–magnetic state, it is necessary to understand
the crystal structure of the precursor phase, i.e., the phase III. O2 is
octahedrally surrounded by Na atoms in all phases. In the phase II, the
molecular axis is aligned along one of the four equivalent [111]–directions of
the octahedron (see Supplemental Material). The NN molecules arrange their
molecular axes to avoid each other, i.e., coherent antiferro–like arrangement
of the molecular axes. Because the symmetry of the octahedron is still cubic,
the degeneracy of the $\pi_{g}^{*}$ orbital should be conserved. On the
contrary, in the phase III, the three–fold symmetry of the octahedron is lost
and the molecular axis is slightly tilted from the [111]–direction of the
octahedron. Na–Na bond lengths are changed to be not equivalent and the Na–O
bond length is increased along the $b$–axis (In Fig.1 (a), different Na–Na
bond lengths are displayed by colors). The NN molecules within the $ac$–plane
arrange their molecular axes to be parallel to each other, namely, the
ferro–like arrangement of the molecular axes is realized. Thus, the NN
molecular axis is parallel to each other along the $c$–axis while those is
canted along the $a$– and $b$–axis. The two–fold distortion of the octahedron
should break the orbital degeneracy, and stabilizes the $\pi_{g}^{*}$ orbital
perpendicular to the $c$–axis as unoccupied orbital. Namely, a ferro–orbital
ordering realizes in the $ac$–plane. The $c$–axis length, i.e., the length
between the NN molecules, is obtained to be 3.39 Å at 100 K, which is close to
the length between molecules in the $\alpha$–phase of the solid–O2 ($\sim 3.2$
Å) O2_Defoits ; O2_Date . Consequently, we can expect the strong AF exchange
interaction along the $c$–axis. This structural peculiarity should be
manifested in the $\chi(T)$ below $T_{2}$, namely, the low–dimensional nature.
Moreover, because there was no structural dimerization of O2 in any directions
in all phases, we were able to deny both isolated dimerization of O2 and 1D
alternating AF chain as the cause of the spin–gap. Accordingly, we can
conclude that the magnetism in the phase III is based on the uniform AF
spin–chain and, then, the magnetic ground state is the spin–Peierls (SP)
state.

Finally, we consider the SP state in NaO2. As the spin–gap value and the SP
transition temperature are obtained to be $\Delta/k_{\rm B}=51.2$ K and
$T_{3}=T_{\rm SP}=34$ K, respectively, from the magnetic measurements, the
value of $2\Delta/k_{\rm B}T_{\rm{SP}}$ is calculated to be 3.01. This is
comparable with the BCS weak coupling result (3.54). The alternating exchange
interaction constants, $J_{1}$ and $J_{2}$, in the SP state can be written as
$J_{1,2}=J(1\pm\delta)$ using the alternating parameter $\delta$ SP_Pytte .
The spin–gap is also related to $\Delta=2pJ\delta$ with $p\sim 1+2/\pi$. Using
the Bulaevskii’s formula for the $\chi(T)$ below the SP transition Bulaevskii
, we estimate the $\delta$ by fitting to the experiment. When we use the
$J/k_{\rm B}=140$ K in the 1D uniform AF phase, which was obtained around 50
K, we obtain the $\delta$ of 0.11. For the organic SP compounds TTF–CuBDT
TTF_CuBDT_Bray  and the inorganic CuGeO3 CuGeO3_Hase , the $\delta$ was
estimated to be 0.167 and 0.167, respectively. The $\delta$ in NaO2 is
comparable to these values. Moreover, the $\delta$ of 0.11 roughly leads to
the $\Delta/k_{\rm B}$ of 51 K, which is almost identical to the spin–gap
value obtained in the experiment. It is known that the AF $J$ of the solid–O2
magnet depends exponentially on the inter–molecular length O2_Hermert ;
O2_Wormer ; O2_Bussery ; O2_Date . As we expect that the direct magnetic
interaction between NN O2s along the $c$–axis is dominant, the same dependence
can be applied. If so, even though the lattice dimerization in NaO2 would be
extremely so small not to be observed experimentally, the $J$–alternation may
occur in the SP phase. More detailed structural study in the SP phase of NaO2
is a future task.

In summary, we have investigated the low–temperature magnetism of NaO2. We
found no magnetic phase transition down to $T=0.3$ K and confirmed the
spin–singlet ground state below $T_{3}$. Raman scattering experiments clearly
indicated the presence of the phase transition at $T_{3}$. Consequently, we
conclude that the singlet–ground state of NaO2 is due to the SP instability.

The authors acknowledge fruitful discussions with H. O. Jescheke, J. Otsuki,
M. Naka, K. Okada, R. Kondo, T. Goto, H. Sagayama, R. Kumai. The x–ray
diffraction study was performed under the approval of the Photon Factory
Program Advisory Committee (Proposal No. 2017G636, 2019T003, 2020G666). The
neutron scattering experiment at the HRC was approved by the Neutron
Scattering Program Advisory Committee of IMSS, KEK (proposals no. 2019S01),
and ISSP. This work was partly supported by JSPS KAKENHI (15H03529, 20K20896,
21H04441), MEXT, Japan.

† Corresponding author<EMAIL_ADDRESS>
§ Present address; Institute for Molecular Science, Okazaki, Japan
¶ Present address; Department of Physics, Institut Teknologi Sepuluh Nopember,
Indonesia
∗ These authors contributed equally to this work.

## References

  * (1) Y. Tokura and N. Nagaosa,  Science, 288, 462 (2000).
  * (2) J. B. Goodenough.  Magnetism and Chemical Bond, (Interscience, New York, (1963).
  * (3) K. I. Kugel and D. I. Khomskii,  Sov. Phys. Usp. , 25(4), 231(1982).
  * (4) J. Kanamori,  J. Phys. Chem. Solids, 10(2), 87 (1959).
  * (5) Y. Ishiguro, K. Kimura, S. Nakatsuji, S. Tsutsui, A. Q R Baron, T. Kimura, Y. Wakabayashi,  Nat. Commun., 4, 2022, (2013).
  * (6) T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y. Tokura.  Nature, 426, 55–58, (2003).
  * (7) S. Ishihara, M. Yamanaka, and N. Nagaosa.  Phys. Rev. B, 56(2), 686–692, 07 1997.
  * (8) S. Margadonna, K. Prassides, H. Shimoda, T. Takenobu, and Y. Iwasa,  Phys. Rev. B, 64, 132414 (2001).
  * (9) T. Kambe, K. Kajiyoshi, M. Fujiwara and K. Oshima,  Phys. Rev. Lett., 99, 177205, (2007).
  * (10) M. Kim and B. I. Min.  Phys. Rev. B, 89(12), 121106(R), (2014).
  * (11) R. Kováčik, P. Werner, K. Dymkowski, and C. Ederer.  Phys. Rev. B, 86(7), 075130, (2012).
  * (12) I V Solovyev.  New J. Phys., 10(1), 013035, (2008).
  * (13) S. Riyadi, B. Zhang, R. A. de Groot, A. Caretta, P. H. M. van Loosdrecht, T. T. M. Palstra, and G. R. Blake.  Phys. Rev. Lett., 108, 217206, (2012).
  * (14) M. Klanjšek, D. Arčon, A. Sans, P. Adler, M. Jansen, and C. Felser.  Phys. Rev. Lett., 115, 057205, (2015).
  * (15) T. Knaflič, M. Klanjšek, A. Sans, P. Adler, M. Jansen, C. Felser, and D. Arčon.  Phys. Rev. B, 91, 174419, (2015).
  * (16) M. Miyajima, F. Astuti, T. Kakuto, A. Matsuo, D. P. Sari, R. Asih, K. Okunishi, T. Nakano, Y. Nozue, K. Kindo, I. Watanabe, and T. Kambe.  J. Phys. Soc. Jpn., 87(6), 063704, (2018).
  * (17) F. Astuti, M. Miyajima, T. Fukuda, M. Kodani, T. Nakano, T. Kambe, and I. Watanabe.  J. Phys. Soc. Jpn., 88(4), 043701, (2019).
  * (18) Fahmi Astuti.  PhD thesis, Hokkaido University, (2019).
  * (19) Mizuki Miyajima.  PhD thesis, Okayama University, (2021).
  * (20) M. Ziegler, M. Rosenfeld, and W. Känzig,  Helv. Phys. Acta, 49, 57, (1976).
  * (21) I. V. Solovyev, Z. V. Pchelkina, and V. V. Mazurenko.  Cryst. Eng. Comm, 16, 522, (2014).
  * (22) A. Zumsteg, M. Ziegler, W. Känzig, and M. Bösch.  J. Phys. Condens. Matter, 17(4), 267, (1974).
  * (23) The crystal structures for these phases were analyzed by the Rietveld refinement and are shown in Supplemental Material.
  * (24) S. D. Mahanti and G. Kemeny.  Phys. Rev. B, 20, 2105, (1979).
  * (25) J. C. Bonner and M. E. Fisher.  Phys. Rev., 135, A640, (1964).
  * (26) B. C. Keith, C. P. Landee, T. Valleau, M. M. Turnbull, and N. Harrison.  Phys. Rev. B, 84, 104442, (2011).
  * (27) P. J. Baker, S. J. Blundell, F. L. Pratt, T. Lancaster, M. L. Brooks, W. Hayes, M. Isobe, Y. Ueda, M. Hoinkis, M. Sing, M. Klemm, S. Horn, and R. Claessen.  Phys. Rev. B, 75, 094404, (2007).
  * (28) E. Ehrenfreund and L. S. Smith.  Phys. Rev. B, 16(5), 1870, (1977).
  * (29) R. S. Hayano, Y. J. Uemura, J. Imazato, N. Nishida, T. Yamazaki, and R. Kubo.  Phys. Rev. B, 20, 850, (1979).
  * (30) I. Watanabe, M. Akoshima, Y. Koike, S. Ohira, and K. Nagamine.  Phys. Rev. B, 62, 14524, (2000).
  * (31) T. Adachi, N. Oki, Risdiana, S. Yairi, Y. Koike, and I. Watanabe.  Phys. Rev. B, 78, 134515, (2008).
  * (32) I. Zaliznyak and S. Lee, Magnetic Neutron Scattering in Modern Techniques for Characterizing Magnetic Materials, edited by Y. Zhu (Springer, Heidelberg, 2005).
  * (33) In Supplemental Material, we estimated the upper limit of a distortion along the $c$–axis, where the strong AF interaction should be present between the NN O2s, by taking into account the systematic error of the experimental background in the neutron elastic diffraction experiments.
  * (34) J. B. Bates, M. H. Brooker, and G. E. Boyd.  Chem. Phys. Lett., 16(2), 391, (1972).
  * (35) M. Bösch and W. Känzig.  Helv. Phys. Acta, 48(5-6), 743, (1975).
  * (36) G. C. DeFotis.  Phys. Rev. B, 23, 4714, (1981).
  * (37) C. Uyeda, K. Sugiyama, and M. Date.  J. Phys. Soc. Jpn., 54(3), 1107, (1985).
  * (38) E. Pytte.  Phys. Rev. B, 10, 4637, (1974).
  * (39) N. Bulaevskii.  Sov. Phys. Solid state, 11(5), 921, (1969).
  * (40) J. W. Bray, H. R. Hart, L. V. Interrante, I. S. Jacobs, J. S. Kasper, G. D. Watkins, S. H. Wee, and J. C. Bonner.  Phys. Rev. Lett., 35(11), 744, (1975).
  * (41) M. Hase, I. Terasaki, and K. Uchinokura.  Phys. Rev. Lett., 70, 3651, (1993).
  * (42) M. C. van Hemert, P. E. S. Wormer, and A. van der Avoird.  Phys. Rev. Lett., 51, 1167, (1983).
  * (43) Paul E. S. Wormer and Ad van der Avoird.  J. Chem. Phys., 81(4), 1929, (1984).
  * (44) B. Bussery, S. Ya. Umanskii, M. Aubert–Frécon, and O. Bouty.  J. Chem. Phys., 101(1), 416, (1994).