\section{Introduction}
\label{sec:intro}
The X-ray astronomy satellite XRISM is planned to be launched in the Japan fiscal year
2022\cite{Tashiro20}. It is expected to revolutionize high-energy astrophysics with its
high-resolution, high-throughput, and non-dispersive X-ray spectrometer based on the
X-ray microcalorimetry. The instrument called \textit{Resolve} is built on the heritage
of Hitomi SXS\cite{Kelley16}, which demonstrated its capability in orbit but was
discontinued a month after launch due to unexpected failure of the spacecraft attitude
control system in 2016 \cite{Takahashi18}.

\textit{Resolve} is a cryogenic instrument with an X-ray microcalorimeter detector in a
Dewar \cite{Ishisaki18, Fujimoto17}. To maintain the Dewar vacuum on ground and to
protect the detector from the initial spacecraft outgassing in orbit, a partially
transparent lid called the gate valve (GV) is installed at the top of the Dewar along
the X-ray optical path.  GV is composed of a Be window and a stainless steel mesh (see
details in Section~\ref{sec:gv}).  On ground, the GV can be opened and closed
repeatedly. In orbit, it can be opened only once using a non-explosive actuator system
and cannot be closed once opened. It is planned to be opened during the commissioning
phase when the outgassing settles. Until then, all observations will be made through the
GV. Therefore, calibrating the energy dependence of the X-ray transmission of this
apparatus is crucial for the initial scientific outcomes.

Be windows are routinely used as a standard entrance window in X-ray detectors requiring
a leak-tight structure. In the past, the modeling based only on the Be photo-electric
absorption was sufficient (e.g. ASCA GIS\cite{Ohashi96}). For X-ray spectrometers with a
high spectral resolution and a high throughput such as Hitomi SXS\cite{Kelley16}, XRISM
\textit{Resolve}\cite{Ishisaki18}, and Athena X-IFU\cite{Barret18}, this is
insufficient. This was recognized in SXS \cite{Yoshida17}, based on which we develop the
model for \textit{Resolve} in this work.

In the SXS, the GV was not calibrated using the flight model on ground unlike the other
components \cite{Eckart18}. This choice was motivated by the lower priority due to its
not being used under nominal observations and some scheduling constraints. Due to the
unexpected early loss of the spacecraft, all observations were made with the GV
closed. After the loss, the transmission measurements were performed using a spare unit
\cite{Eckart18, Yoshida17, Hoshino17}. As a result, the GV turned out to be the cause
for the largest systematic uncertainty in the effective area calibration in the Crab
observation with the SXS \cite{Tsujimoto18}.

Based on this lesson, we planned to perform the ground calibration measurements for all
parts of the GV for \textit{Resolve} with a clear set of calibration requirements that
are traced back to science goals of the mission.  We aim to understand causes of the
large systematic uncertainty in the SXS GV and reduce it.  The allocated budgets for the
GV are shown in Table~\ref{tab:requirements}.  The calibration results will be delivered
to users through the calibration database (CALDB).

\begingroup
\renewcommand{\arraystretch}{1.3}
\begin{table}[H]
\caption{Calibration requirements of the Be window and the stainless steel mesh\cite{Eckart19req}}\label{tab:requirements}
\scalebox{0.72}{
\begin{tabular}{cccccl}\hline
Component & Items & Accuracy & Range (req) & Range (goal) &\multicolumn{1}{c}{Comment}\\\hline\hline
Be window& \begin{tabular}{c}Broadband\\\relax (absolute transmission) \end{tabular} & $\pm$1\% & 2.5--12.0~keV & 1.8--25.0~keV & \\\cdashline{2-6}[3pt/1.4pt]
& \begin{tabular}{c}Fine structures\\\relax (relative transmission) \end{tabular} & $\pm$5\% & & & \begin{tabular}{l}BDFs and contaminants\\ Knowledge of structure to 2~eV required\end{tabular}\\\cdashline{2-6}[3pt/1.4pt]
& \begin{tabular}{c}Spatial uniformity\\\relax (relative transmission) \end{tabular} & $\pm$3\% & \begin{tabular}{c}single energy\\ ($\sim$3.6~keV)\end{tabular} && \begin{tabular}{l}With $>$100 spatial grids\\ \end{tabular}\\\hline
Mesh & Broadband & $\pm$1\% & 2.5--12.0~keV & 1.8--25.0~keV & \\
\hline
\end{tabular}
}
\end{table}
\endgroup

We performed a ground calibration campaign using various facilities in 2018--2019 and
achieved the required calibration accuracy before the GV was installed in the flight
model (FM) of \textit{Resolve} at the end of 2019. In this paper, we describe the
details of this campaign in the following structure. First, we give a brief description
of the GV in Section~\ref{sec:gv}. We then present the results of the measurements in
Section~\ref{sec:results}. Based on the results, we construct a physical model of the
transmission curve and apply it to the data showing that the best-fit model is
consistent with the data within the required calibration accuracy. Finally, in
Section~\ref{sec:discussion}, we show that we can obtain a better fit to the observed
Crab spectrum with the SXS by applying our model with improved parameters.
We summarize the results in Section~\ref{sec:summary}. More details can be found in Midooka (2020)\cite{Midooka20}.

\section{GATE VALVE}
\label{sec:gv}
The GV has two parts in the optical path when it is closed (Figure
\ref{fig:GV_top_dewar}): (1) Beryllium window of about 270~\textmu m\xspace thickness and (2)
stainless steel mesh for protecting the Be window from accidental damage by any
infalling materials during the ground testing. A picture of these components is shown
in Figure~\ref{fig:Be_mesh}. They have an inner diameter of 29~mm and are placed about
231~mm above the detector.

\begin{figure}[H]
\begin{center}
	\includegraphics[keepaspectratio, width=0.5\columnwidth]{GVclose.pdf}
	    \caption[\textit{Resolve} FM GV at the top of the Dewar]{\textit{Resolve} FM GV in a closed configuration at the top of the Dewar. The photo was provided by Sumitomo Heavy Industries, Ltd.}
	    \label{fig:GV_top_dewar}
\end{center}
\end{figure}

\begin{figure}[H]
\begin{center}
	\includegraphics[width=0.7\columnwidth]{Be_mesh.pdf}
	    \caption[Picture of the Be window and the stainless steel mesh]{Picture of the Be window (calibration model placed on the holder) and the stainless steel mesh (flight model attached to the flange).}
	    \label{fig:Be_mesh}
\end{center}
\end{figure}

\paragraph{Be window}
The Beryllium window was fabricated by Materion Corp. from a high purity (IF-1) rolled
foil with a reported thickness of 270$_{-5}^{+8}$~\textmu m\xspace. Three models were provided from
the same lot, which we use as the flight model (FM), the flight spare (FS), and the
calibration model (CM). Although the geometrical thickness was measured by the provider,
the X-ray measurements are mandatory because (i) the Be window in a thin film form may
have a lower density than the values in the literature, (ii) the thickness has some
gradient over the diameter in the manufacturing process, and (iii) many spectroscopic
features (absorption edges by atoms of impurities and the features caused by Bragg
diffraction) may affect the high-resolution science by the microcalorimeter. They are
known to vary in different lots from the SXS results \cite{Hoshino17,Yoshida17}.

\paragraph{Stainless steel mesh}
The stainless steel mesh was installed by Sumitomo Heavy Industries, Ltd. The mesh is a
plain weave using the SUS~304 stainless steel wires with a circular cross section. In
the drawing, the wire thickness is 0.20~mm, the pitch is 1.27~mm, and the wire gap is
1.07~mm. The design values are the same for SXS and \textit{Resolve}. From these
values, the opening fraction is estimated to be $(1.07)^2/(1.27)^2 \sim 71\%$. The
remaining 29\% is opaque to the X-rays except for the highest end of the X-ray
bandpass. The mesh for SXS had an optically-thick cross structure at the center for
alignment purposes, which was removed for \textit{Resolve}. The transmission by the
presence of this structure is included in the ray-tracing simulation, which is treated
separately from the average transmission of the mesh part to be presented in this
manuscript. The mesh part was not measured in the SXS.

\section{EXPERIMENTS AND RESULTS}
\label{sec:results}

\subsection{Energy Dependence of the X-ray Transmission of the Stainless Steel Mesh}
We first describe the X-ray transmission measurements of the stainless steel mesh
performed using the 30~m X-ray beamline at JAXA ISAS \cite{Hayashi15} from February 4 to
19, 2019. At that time, the FM stainless steel mesh was not manufactured, thus we used a
spare mesh left from the SXS. Figure~\ref{fig:ISAS_BL} shows the schematic view of the
beamline. X-rays from an X-ray generator illuminate one of the five metal targets: Ti,
Cu, Pt, Mo, and Ag. Their fluorescent and scattered X-rays are further monochromatized
using a Ge (220) double crystal monochromator (DCM). The stainless steel mesh was
installed on the sample stage. The parallel beam was formed by a slit of 2.0 $\times$
2.0~mm, which is comparable in size with the mesh pitch.

\begin{figure}[H]
\begin{center}
	\includegraphics[width=1.0\columnwidth]{ISAS_BL.pdf}
	    \caption[Schematic of the ISAS beamline]{Schematic of the ISAS beamline \cite{Hayashi15}. The distance between X-ray generator and a slit is approximately 27~m. The slit, sample stage, and detector stage are installed in the measuring chamber.}
	    \label{fig:ISAS_BL}
\end{center}
\end{figure}

\begin{figure}[H]
\begin{center}
	\includegraphics[width=0.7\columnwidth]{ISASchamber.pdf}
	    \caption[Interior of the measuring chamber in the beamline]{Interior of the measuring chamber in the beamline.}
	    \label{fig:ISAS_chamber}
\end{center}
\end{figure}

X-ray spectrum with and without the mesh on the sample stage was measured with the
proportional counter (PC) filled with a P10 gas (Ar: 90\%, CH$_4$: 10\%) on the detector
stage. We call them transmitted and direct X-rays, respectively. For the transmitted
X-rays, we performed raster scans over the entire mesh to average for the mesh
pattern. A raster scan is composed of several zigzag lines over entire surface of the mesh by moving the
sample stage continuously at a slow speed of 0.1~mm~s$^{-1}$. The entire surface is
scanned many times to make the systematic error small enough. We performed direct X-ray
measurements before and after the transmitted X-ray measurements without moving the
stages. This set was repeated several times to accumulate sufficient photon statistics. We
then calculated the transmission of the mesh by taking the ratio between the transmitted
and direct X-ray photons in an energy band of interest for each target line. The result
is shown in Table~\ref{tab:meshtable} and the red points with errors in
Figure~\ref{fig:mesh_curve}.  Measurement errors include systematic and statistical
errors, though the latter is negligible. The former was estimated as variations among
individual scans and the latter as Poissonian of the photon counts.

\begingrou
\renewcommand{\arraystretch}{1.2
\begin{table}[H]
\begin{center}
\caption{Transmission measurement results of the stainless steel mesh
\scalebox{0.87}{
\begin{threeparttable}
\begin{tabular}{ccccccc}\hline
Target & \begin{tabular}{c}Energy\\\relax[keV] \end{tabular}  & Sample\tnote{a} & \begin{tabular}{c}Total exp.\\\relax [sec] \end{tabular} & \begin{tabular}{c}Total photons\\\relax [counts] \end{tabular}  & \begin{tabular}{c}Transmission\\ \relax[\%] \end{tabular} & \begin{tabular}{c}Error\tnote{b}\\\relax [\%] \end{tabular}\\\hline
\multirow{2}{*}{Ti K{\rm $\alpha$}} & \multirow{2}{*}{4.51} & ON & 1866.3 & 2562596 & \multirow{2}{*}{72.7} & \multirow{2}{*}{$\pm$0.6}\\
 & & OFF & 3627.8 & 6850815 & & \\  \hline
\multirow{2}{*}{Cu K{\rm $\alpha$}} & \multirow{2}{*}{8.05} & ON & 1856.7 & 200030 & \multirow{2}{*}{71.7} &  \multirow{2}{*}{$\pm$0.9}\\
 & & OFF & 3601.7 & 540877 & & \\ \hline
\multirow{2}{*}{Pt L{\rm $\beta$}} & \multirow{2}{*}{11.07} & ON & 1858.0 & 422722 & \multirow{2}{*}{71.5} &  \multirow{2}{*}{$\pm$1.5}\\
 & & OFF & 3605.6 & 1149485 & & \\ \hline
\multirow{2}{*}{Mo K{\rm $\alpha$}} & \multirow{2}{*}{17.48} & ON & 3712.8 & 211733 & \multirow{2}{*}{72.3} & \multirow{2}{*}{$\pm$2.3}\\
 & & OFF & 3601.6 & 285442 & & \\ \hline
 \multirow{2}{*}{Mo K{\rm $\beta$}} & \multirow{2}{*}{19.61} & ON & 11136.2 & 122811 & \multirow{2}{*}{73.6} &  \multirow{2}{*}{$\pm$2.8}\\
 & & OFF & 7200.1 & 109824 & & \\ \hline
\multirow{2}{*}{Ag K{\rm $\alpha$}} & \multirow{2}{*}{22.17} & ON & 1861.2 & 270694 & \multirow{2}{*}{76.4} & \multirow{2}{*}{$\pm$1.3}\\
 & & OFF & 3613.6 & 686769 & & \\ \hline
\end{tabular}
	\begin{tablenotes}
		\item[a] \small{Sample ON/OFF means the measurements of the transmitted/direct X-rays.}
		\item[b] \small{The error is for a 1$\sigma$ level.}
	\end{tablenotes}
\end{threeparttable}
}
\label{tab:meshtable}
\end{center}
\end{table}
\endgroup

\subsection{Energy Dependence of the X-ray Transmission of the Be Window}
\label{sec:E_Be}

X-ray transmission measurements of the FM Be window were performed at KEK Photon Factory
(PF) from June 13 to 23, 2019. We used the beamlines BL-7C\cite{7c} and
BL-11B\cite{11b}. The operating energy ranges are 4.0--12.0~keV and 2.1--5.0~keV and we
used a beam size of 1.0$\times$1.0~mm and 2.5$\times$2.5~mm as the largest possible
square respectively for BL-7C and BL-11B. The intense synchrotron X-rays are
monochromatized using a Si (111) DCM. We measured three positions of the Be window; they
are at the center of the window and at two positions 6.5~mm off-center (top and bottom).

Figure~\ref{fig:PF_sche} shows the schematic drawings of the setup for measurements at
the BL-11B and BL-7C\cite{Yoshida17}. In the BL-11B experiments, we attached our vacuum
chamber to the beamline through a Kapton polyimide film and installed the Be window and
a photo diode inside the chamber to prevent absorption by the air as shown in
Figure~\ref{fig:BL11B}. The incident X-rays were measured using a Ti foil as a
photocathode beam monitor provided by the facility. In the BL-7C experiments, we
measured incident and transmitted X-rays by two ion chambers filled with N$_2$ gas
provided by the facility in front and back of the Be window placed in the air.

The contamination by higher-order X-rays are removed by a focusing double mirror in the
BL-7C measurements. The fraction of the remaining third-order X-ray contamination at
4~keV is negligible of about $3 \times 10^{-5}$ \cite{7c}. The contamination has not
been reported in the BL-11B in past experiments at least above 2.3~keV, hence we
consider that they are negligible at the higher energy range in our measurements.

\begin{figure}[H]
 \begin{center}
  \includegraphics[width=0.9\columnwidth]{sche.pdf}
  \caption[Setup for measurements at the BL-11B and BL-7C]{Schematic drawings of the
  setup for measurements at the BL-11B and BL-7C\cite{Yoshida17}.}
  \label{fig:PF_sche}
 \end{center}
\end{figure}

\begin{figure}[H]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{BL11B.pdf}
  \caption[The photo of the BL-11B]{right: The photo of the BL-11B. left:  Interior of
  the vacuum chamber. Be window and a photo diode as the back-side detector are
  installed.}
  \label{fig:BL11B}
 \end{center}
\end{figure}

In both experiments, we measured current of the detectors (i.e. the ion chambers, the
photo-diode, and the Ti foil), as the intensity of the incident and transmitted X-rays.
Here, we assume that $I_0(E)$ and $I_1(E)$ are the
incident and transmitted X-ray intensities with the Be window, and $I'_0(E)$ and
$I'_1(E)$ are those without the Be window at a given energy \textit{E}. The dark current
was measured and subtracted from them. We derived the transmission by
\begin{equation}
T(E)=\frac{I_{1}(E) / I_{0}(E)}{I_{1}^{\prime}(E) / I_{0}^{\prime}(E)}.
\label{eq:trans}
\end{equation}
With this double normalization, we can correct for any instability of the beam as well
as the energy dependence caused by the measurement setup. We repeated this measurement
with a step of 2~eV in eight energy ranges overlapping with each other at 2.35--2.6,
2.55--3.0, 2.95--3.5, 3.45--4.0~keV bands in BL-11B, and 4.0--5.0, 4.9--8.0, 7.9--10.0,
9.9--12.0~keV bands in BL7C. The overlapped range is 50 or 100~eV.

The upper panels of Figure~\ref{fig:11B_7C} show derived transmission curves of the top
position of the Be window and model curves in 2.1--4.5~keV at the BL-11B (left) and in
4.0--12.0~keV at the BL-7C (right). The ratio between two adjacent energy ranges was
corrected by multiplying a constant to the data set in the lower energy of the two. This
was repeated to construct a continuous transmission curve over the entire energy. The
seven constant values provide an estimate of the systematic uncertainty in the
measurements. The mean and the standard deviation of them are 1.00000 and 0.00081,
respectively. We thus evaluate the systematic uncertainty of the measurements to be
0.08\%.

\begin{figure}[H]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{curve_11B_7C.pdf}
  \caption[Transmission curve in 1.8--12.0~keV]{Upper panel: A model transmission curve
  assuming a photoelectric absorption by Be of a 250~\textmu m\xspace thickness and the transmission
  measured at the top position of the Be window obtained at the BL-11B (left) and BL-7C
  (right). Lower panel: The ratio of the measurement curve to the model curve. The anomalous
  increase of the transmission below 2.4~keV is attributable to higher-order X-ray
  reflections in the DCM.}
  \label{fig:11B_7C}
 \end{center}
\end{figure}

The lower panels of Figure~\ref{fig:11B_7C} show the ratio of the measurement against
the model assuming photoelectric absorption by Be of a 250~\textmu m\xspace thickness. Here, we used
the thickness only for displaying purpose to show the residual structures in a limited
range. We found an apparent increase in the transmission below 2.4~keV. After the
delivery of FM Be window, we made an additional measurement using another model (CM) to
investigate this issue at the beamline BL-11 of the Hiroshima Synchrotron Radiation
Center (HiSOR)\cite{hisor} in October 2019. In HiSOR, we examined the contribution of
higher-order X-rays by detuning the two crystal angles in the DCM. The result suggested
some influence of higher-order X-rays in the measurements at KEK PF. We conclude that
the anomalous increase of the transmission below 2.4~keV (Figure.~\ref{fig:11B_7C}) is
likely caused by the effect of higher-order X-ray reflections \cite{Midooka20} and
decided to ignore the data below 2.6~keV in the following analysis.

\subsection{Spatial Dependence of the X-ray Transmission of the Be Window} \label{sec:spatial}
We evaluated the spatial dependence of the transmission using the FM Be window. We used
the beamline BL-11B using almost the same setup in Section~\ref{sec:E_Be}. The beam size
is adjusted to 0.8~mm (horizontal) $\times$ 1.2~mm (vertical).  We measured the X-ray
transmission at some fixed energy over the entire Be window by a scan. We chose four
energies: 2700, 3000, 3456, and 4000~eV. The scan was continuous in the horizontal
direction at a speed of 500~\textmu m\xspace~s$^{-1}$ and step-wise in the vertical direction with a
step size of 1~mm. The time series in the horizontal scan was converted to the positions
by the speed and rebinned into a grid map. At the beginning and at the end of each
horizontal scan, we obtained reference transmission without the Be window for 10~s.

Figure~\ref{fig:BL11B_raster} shows the 2D maps of the transmission of the FM Be window.
In all energies, transmission at the lower left of the Be window is higher than the
average value. We speculate that this is caused in the manufacturing process. In the
4000~eV map, the mean and 3\,$\sigma$ scatter of the transmission is
$66.8^{+1.0}_{-1.1}$ \%, which translates to a thickness of 251.5$\pm$10.5~\textmu m\xspace. The
non-uniformity is $\pm$4.2\% at 3\,$\sigma$.

\begin{figure}[H]
 \begin{center}
  \includegraphics[width=0.78\columnwidth]{BL11B_raster.pdf}

  \caption[Spatial non-uniformity maps of the Be window transmission]{Spatial
  non-uniformity maps of the transmission of the FM Be window at four selected
  energies. The figures are aligned to the spacecraft coordinate: i.e, $+$X, $+$Y
  directions in the figure are aligned to SAT$+$X  and SAT$+$Y. These maps were faced in
  the direction viewed from the incident X-rays both in the experiment and also to the
  in-orbit setup. The fuzzy boundary is due to an artifact as a result of circular masking.}
  \label{fig:BL11B_raster}
 \end{center}
\end{figure}

\section{DISCUSSION} \label{sec:discussion}

\subsection{Modeling Transmission of the Stainless Steel Mesh} \label{sec:model_mesh}
In Figure~\ref{fig:mesh_curve}, the measured transmission (Table~\ref{tab:meshtable}) is
shown with the red points with errors, whereas the model curve adopted in the SXS CALDB
is with the blue curve. The mesh is made of SUS~304, which we assumed to have 70\% Fe,
20\% Cr, and 10\% Ni. The transmission is modeled by a product of the photoelectric
absorption and scattering attenuation of these elements.

The model curve adopted in the SXS was found inconsistent with the measurements, in
particular, for the low energy end in Ti K\rm $\alpha$\ (4.51~keV) and the high energy end
in Ag K\rm $\alpha$\ (22.17~keV). One reason is that the cross section of the mesh wire was
assumed a square with 0.2~mm $\times$ 0.2~mm in SXS, while it is round in reality, hence
the actual average depth is smaller than 0.2~mm.  The transmission is non-linear to the
depth, so we performed numerical integration assuming that the wire depth is
$2\sqrt{(t^2/4 - x^2 )}$ at $x (<t/2)$ away from the center of the wire with a diameter
of $t=0.2$~mm. Using the mass attenuation coefficient of a constituent metal
$\mu_\mathrm{tot,metal} (E)$, we derived the total attenuation coefficient $\alpha (E)$
as
\begin{equation}
\alpha(E)=\sum_{\mathrm{metal}} A_{\mathrm{metal}} * \rho_{\mathrm{metal}} * \mu_{\mathrm{tot,metal}}(E),
\end{equation}
where $A_\mathrm{metal}$ is the mass fraction and $\rho_\mathrm{metal}$ is the mass
density of the constituents (Fe, Cr, Ni). The transmission of the stainless steel mesh
$T_\mathrm{mesh}$ ($E$) is then calculated as:
\begin{equation}
T_{\mathrm{mesh}}(E)=f+(1-f) \int_{0}^{\frac{t}{2}} \exp \left(-2 \sqrt{t^{2} / 4-x^{2}} * \alpha(E)\right) d x,
\end{equation}
where \textit{f} is the opening fraction of the mesh. We note that the contribution of
the overlapped parts of the mesh wires can be ignored compared to the measurement error
in Table~\ref{tab:meshtable}.

First, we fixed \textit{f} to 0.71 estimated based on the design in
Section~\ref{sec:gv}, which is shown with the orange curve in
Figure~\ref{fig:mesh_curve}. This is still inconsistent with the
measurements. Considering some uncertainty in the estimate of \textit{f}, we treated it
as a free parameter. The best-fit value was 0.723$\pm$0.003. The resultant model is
shown with the green curve in Figure~\ref{fig:mesh_curve}, which is now consistent with
the measurements. We adopt this model for \textit{Resolve}.

\begin{figure}[H]
 \begin{center}
  \includegraphics[width=0.8\columnwidth]{mesh_curve.pdf}
  \caption[Comparison of three transmission models for the stainless steel
  mesh]{Comparison of three transmission models for the stainless steel mesh. Red points
  describe the measurement results with errors. Blue, orange, green curves show the SXS
  CALDB model, a round cross-section model with a fixed 71\% opening fraction, and the
  \textit{Resolve} model.}
  \label{fig:mesh_curve}
 \end{center}
\end{figure}

\subsection{Modeling Transmission of the Be Window} \label{sec:model_Be}
First, we investigated contamination by impurities. In the lower panel of Figure~\ref{fig:11B_7C}, we identified features associated with the photoelectric absorption edges of Cr-K, Mn-K, Fe-K, Ni-K, and Cu-K, but not others. We derived their areal densities by local fitting and summarize the result in Table~\ref{tab:minor_thick}.
The result is compared with the report by Materion Corp.

\begin{table}[H]
\caption[Areal density of impurities of the FM Be window]{Best-fit areal density of the impurities based on the transmission measurements of the FM Be window.}
\label{tab:minor_thick}
\begin{center}
\begin{threeparttable}
\begin{tabular}[t]{cccccc}\hline
Impurity  & Edge energy & Fitting energy & Best-fit areal\tnote{a} & Reference areal\tnote{b}\\
element & [eV] & range [eV] &  density [\textmu g\xspace/cm$^2$] &  density [\textmu g\xspace/cm$^2$]\\\hline
Cr & 5989 & 5500--6050 & 0.8 $\pm$ 0.1 & 1\\
Mn & 6539 & 6510--6900 & $<$ 0.2 & 1\\
Fe & 7112 & 7000--7500 & 11.4 $\pm$ 0.2 & 9\\
Ni & 8333 & 8200--8450 & 7.7 $\pm$ 0.2 & 5\\
Cu & 8979 & 8800--9180 & 0.7 $\pm$ 0.1 & 1\\\hline
\end{tabular}
\begin{tablenotes}
\item[a] \small{The uncertainties are at the 1$\sigma$ confidence level.}
\item[b] \small{The values are derived from the fraction of materials measurements by Materion Corp.}
\end{tablenotes}
\end{threeparttable}
\end{center}
\end{table}

We also detected other edge-like features at 3460, 6057~eV, and so on. They are physically interpreted as the Bragg diffraction in the Be window \cite{Yoshida17,Hoshino17}. Here, the Be window is a poly-crystal, in which a part of the material satisfies the Bragg condition at any given energy. An ideal case is called the powder diffraction and the diffracted X-rays form a ring called the Debye-Scherrer ring. Hereafter, we call these edge-like features as Bragg diffraction features (BDFs).
When the incident X-ray energy satisfies the Bragg condition expressed in
Equation~(\ref{eq:bragg}) for the lattice spacing $d_{hkl}$ of the Miller index ($hkl$),
X-rays are coherently scattered into a Debye-Scherrer ring. The ring disappears when the
incident energy becomes lower than the BDF energy $E_{hkl}(\theta = \pi/2)$, as the
scattered angle cannot exceed $2\theta = 180^{\circ}$ . This makes a discontinuity in
the measured transmission curve at $E_{hkl}(\theta = \pi/2)$ of all Miller indices
except for those with systematic absence (Appendix~\ref{sec:appendix}).
\begin{equation}
E_{hkl}(\theta)=\frac{h c}{2 d_{hkl} \sin \theta},
\label{eq:bragg}
\end{equation}
where $\theta$ is the diffraction angle, $h$ is the Planck constant, and $c$ is the speed of light.

The interpretation was confirmed in an additional measurement using CM. We performed an
X-ray diffraction measurement using a multi-purpose X-ray diffractometer (Rigaku Ultima
IV) in the Tokyo Metropolitan University. Figure~\ref{fig:xrd_fit} shows the intensity
of the scattered Cu K$\alpha$ line X-rays as a function of the diffraction angle. All
peaks are identified with a Miller index and correspond to a Debye-Scherrer ring. The
peak heights in the diffraction is in an agreement with the edge depths in the
transmission for each Miller index. The result endorses our interpretation that the
features in the transmission curve are indeed the BDFs.

\begin{figure}[H]
\begin{center}
\includegraphics[keepaspectratio, scale=0.6]{xrd_fit_index.pdf}
\caption[Results of the X-ray diffraction measurements and fitting]{Results of the X-ray diffraction measurements and fitting. All peaks are identified with a Miller index at an expected diffraction angle. All peaks have double-peak structures due to mixing of the Cu K\rm $\alpha$$^1$ and K\rm $\alpha$$^2$ lines in the incident X-rays.}
\label{fig:xrd_fit}
\end{center}
\end{figure}

We took the following approach in modeling the BDFs. We used the thickness of the photoelectric absorption, coherent scattering, and incoherent scattering as phenomenological fitting parameters.
The energy dependence of the mass attenuation coefficients for the three interactions was fixed to the values by the NIST database\cite{xcom}. This approach was motivated by the fact that the scattering cross section effectively changes in lattice from the values in atoms in dynamical diffraction theories \cite{Authier12}. Considering that the BDF is caused by the energy dependence of the anisotropy of coherent scattering, we changed the coherent scattering thickness in a step-wise manner at all BDFs. In between two BDFs, we assumed that the thickness is proportional to the inverse square of the incident energy (Figure~\ref{fig:thick_coh}). This is derived from the discussion\cite{oxford} on the cross section of the powder diffraction as

\begin{equation}
\label{e05}
\frac{d\sigma}{d\Omega} = N_{c} \left(\frac{d_{\tau}^3}{v_0}\right) M_{\tau}|F(\tau)|^{2}\frac{\tan{\theta}}{2\pi} \delta(\gamma),
\end{equation}
in which $N_c$ is the number of unit cells in the powder crystal, $v_0$ is the volume of
a unit cell, $\tau$ is the reciprocal vector corresponding to a Miller index, $F(\tau)$
and $M_{\tau}$ are respectively the structure factor and the reflection multiplexity
factor for the reflection $\tau$, and $d_{\tau}$ is the lattice spacing. The Bragg condition is met at $hc/E = 2d_{\tau}\sin{\theta}$. The angle $\gamma=2\theta$ and $d\Omega = 2\pi \sin{\gamma}d\gamma$.
The Equation~(\ref{e05}) is derived for the constant-wave case, but it applies to the
constant-energy case as well as wavelength-dispersive case with variable
conversions. The scattered photon counts per second integrated over the entire solid
angle (or entire energy range through the Bragg condition) is then given by
\begin{equation}
P_{\tau}(E) = \int dE' \frac{d\sigma}{d\Omega}\frac{d\Omega}{dE'} \Phi(E') = N_{c} \left(\frac{d_{\tau}}{v_0}\right) M_{\tau}|F(\tau)|^{2}\Phi(E)\frac{h^2c^2}{4E^2},
\end{equation}
in which $\Phi(E)$ is the incident photon counts per second per unit area per energy
bin. From these, $P_{\tau}(E)/\Phi(E) \propto E^{-2}$, which translates to the thickness
dependence on $E$ when the scattering depth is optically thin. The thickness was left
constant against energy beyond $\sim$9.2~keV, where there are numerous BDFs and we
observed no discontinuous structures in the data.

Table~\ref{tab:coh_depth} shows best-fit coherent scattering thickness of the
\textit{Resolve} FM Be window, which accounts for BDFs.  Considering these values and
the dependency of $E^{-2}$, we made the plot of the blue curve in Figure~\ref{fig:thick_coh}.

\begin{table}[H]
\begin{center}
\caption[Fitting results of the coherent scattering thickness of the \textit{Resolve} FM
 Be window]{Fitting results of the coherent scattering thickness at each BDF for the
 \textit{Resolve} FM Be window.}
\label{tab:coh_depth}
\begin{threeparttable}
\scalebox{0.8}{
\begin{tabular}[t]{ccccc}\hline
\begin{tabular}{c}Miller index\\\relax $hkl$\end{tabular} & \begin{tabular}{c}BDF energy\\\relax [keV]\end{tabular} & \begin{tabular}{c}Fitting range\\\relax [keV]\end{tabular} & \begin{tabular}{c}Coherent scattering \\\relax thickness $t_{\rm Be, coh}$ [\textmu m]\end{tabular} & \begin{tabular}{c}$t_{\rm Be, coh}$ difference \tnote{b}\\\relax above/below BDF energy [\textmu m]\end{tabular}\\\hline
100 & 3.1326  & N/A & 0\tnote{a} & 0\\
002 & 3.4596  & 3.47--3.57 & 115.2$\pm$1.6 & 115.2\\
101 & 3.5785  & 3.59--4.66 & 150.5$\pm$0.8 & 42.8\\
102 & 4.6671  & 4.67--5.42 & 123.5$\pm$0.8 & 35.0\\
110 & 5.4259  & 5.43--6.05 & 88.5$\pm$0.7 & 0.02\\
103 & 6.0616  & 6.07--6.25 & 158.2$\pm$1.1 & 85.0\\
200 & 6.2653  & 6.27--6.43 & 156.0$\pm$0.9 & 7.9\\
112 & 6.4350  & 6.44--6.49 & 159.4$\pm$1.4 & 11.3\\
201 & 6.4999  & 6.51--6.91 & 167.9$\pm$0.6 & 11.9\\
004 & 6.9191  & 6.93--7.1 & 171.2$\pm$0.6 & 23.0\\
202 & 7.1569  & 7.16--7.59 & 178.0$\pm$0.5 & 18.1\\
104 & 7.5953  & 7.6--8.13 & 172.7$\pm$0.4 & 14.6\\
203 & 8.1353  & 8.14--8.78 & 178.6$\pm$0.9 & 28.1\\
114 & 8.7929  & 8.8--9.19 & 192.9$\pm$0.5 & 40.1\\
105 & 9.1988  & 9.2--12.0 & 208.0$\pm$0.2 & 31.7\\\hline
\end{tabular}
}
 \begin{tablenotes}
  \item[a] We set $t_{\rm Be, coh}=0$ below the lowest BDF energy at the (002) plane.
  \item[b] These values are used to plot blue curve in Figure~\ref{fig:thick_coh}.
 \end{tablenotes}
\end{threeparttable}
\end{center}
\end{table}

\begin{figure}[H]
\begin{center}
    \includegraphics[width=0.9\columnwidth]{t_coh.pdf}
    \caption[Best-fit coherent scattering thickness of the Be with step-wise
 changes]{Best-fit thickness of coherent scattering of Be with step-wise changes. Purple
 lines represent the BDF energies apparent in the data, and the numbers in parentheses
 represent the crystal planes $(hkl)$ that cause the BDF. We summarize the actual
 $t_{\mathrm{Be, coh}}(E)$ values in Table~\ref{tab:coh_depth}.}
	\label{fig:thick_coh}
\end{center}
\end{figure}

Our model consists of the following components:
\begin{itemize}
    \item Photoelectric absorption and incoherent scattering by the bulk Be to account for the overall shape of the transmission curve.
    \item Photoelectric absorption by impurities (Cr, Mn, Fe, Ni, and Cu) to account for the absorption edges.
    \item Coherent scattering with step-wise changes in thickness to account for the BDFs of Be poly-crystals.
\end{itemize}

The total transmission is given by
\begin{equation}
T_{\mathrm{Be}}(E)=T_{\mathrm{Be, photo+inc}}(E) * T_{\mathrm{imp, photo}}(E) * T_{\mathrm{BDF}}(E),
\label{eq:three_models}
\end{equation}
where $T_\mathrm{Be,photo+inc} (E)$ is the transmission due to the photoelectric absorption and incoherent scattering by the bulk Be expressed as:
\begin{equation}
T_{\mathrm{Be, photo+inc}}(E)=\exp \{-t_{\mathrm{Be, photo+inc}} * \rho_{\mathrm{Be}} * (\mu_{\mathrm{Be, photo}}(E) + \mu_{\mathrm{Be, inc}}(E)) \}.
\label{eq:photoinc}
\end{equation}
Here, $t_\mathrm{Be,photo+inc}$ is the thickness of the bulk Be common to the photoelectric absorption and the incoherent scattering. $\rho_\mathrm{Be}$ is the Be density. $\mu_\mathrm{Be,photo} (E)$ and $\mu_\mathrm{Be,inc} (E)$ are, respectively, the mass attenuation coefficient of photoelectric absorption and incoherent scattering as a function of incident X-ray energy $E$.
$T_\mathrm{imp,photo} (E)$ is the transmission due to the photoelectric absorption by the impurities, which is calculated by
\begin{equation}
T_{\mathrm{imp,photo}}(E)=\prod_{\mathrm{imp}} \exp \{-\sigma_{\mathrm{imp}} * \mu_{\mathrm{imp,photo}}(E)\},
\label{eq:imp}
\end{equation}
where $\sigma_{\mathrm{imp}}$ is the areal density of the impurities present in the bulk Be. We modeled the transmission of the BDFs by
\begin{equation}
T_{\mathrm{BDF}}(E)=\exp \{-t_{\mathrm{Be, coh}}(E) * \rho_{\mathrm{Be}} * \mu_{\mathrm{Be, coh}}(E)\}.
\label{eq:bdf}
\end{equation}

Figure~\ref{fig:Be_best} shows the best-fit model for the \textit{Resolve} FM Be
window.
The best-fit parameters are: $T_{\mathrm{Be, photo+inc}}(E)=252.5$ \textmu m\xspace for Equation~\ref{eq:photoinc}, values in Table~\ref{tab:minor_thick} for Equation~\ref{eq:imp}, and those in Table~\ref{tab:coh_depth} for Equation~\ref{eq:bdf}.
The difference between the data and the model is smaller than $\sim$0.3\% over a
wide energy range, which satisfies the calibration requirements described in
Section~\ref{sec:intro}.  The mass attenuation coefficients by NIST have some
uncertainties. However, they are accounted to some extent by being multiplied by a free
parameter. In total, we estimate the transmission uncertainty smaller than $\sim$0.3\%
in 2.6--12.0~keV.

As shown in Figure~\ref{fig:BL11B_raster}, the Be thickness has a spatial non-uniformity
of $\pm 4.2\%$. We plan to include this in CALDB as a spatial map. Other spatial
differences, such as the absorption edges by the impurities and the BDFs, were found
small enough not to affect the observation results based on the comparison of the
measurements at three positions (center, top, and bottom) of the Be window.

\begin{figure}[H]
\begin{center}
	\includegraphics[width=0.8\columnwidth]{2600_12000_best.pdf}
	    \caption[Data and best-fit transmission model of the FM Be window in 2.6--12.0~keV]{Data and best-fit transmission model of the FM Be window in 2.6--12.0~keV (top). The difference between the data and the model (bottom) is smaller than $\sim$0.2\% in the entire measured energy range. The edges by the contaminants are shown in green, while the BDFs of the labeled Miller index is in purple.}
	    \label{fig:Be_best}
\end{center}
\end{figure}

Figure~\ref{fig:eachmodel} shows the three transmission components expressed in
Equation~(\ref{eq:three_models}). Photoelectric absorption and incoherent scattering of
Be (blue) dominates the overall transmission curve. The BDFs (green) and photoelectric
absorption of impurities (orange) account for local edge-like features expanded in the
inset.

\begin{figure}[H]
\begin{center}
 \includegraphics[width=0.8\columnwidth]{eachmodel.pdf}
 \caption[Decomposition of the three transmission model components]{Comparison of the three transmission model components in Equation~(\ref{eq:three_models}). Inset is an enlarged view of the area where the y-axis is greater than 0.875.}
 \label{fig:eachmodel}
\end{center}
\end{figure}

The BDF model in this study is quite similar with that developed for SXS by Yoshida et
al. (2017)\cite{Yoshida17}. Both have step-wise discontinuities followed by a $E^{-2}$
dependency on photon energy $E$. We gave some physical explanation for this behavior and
could impose some constraints global to the energy range. Together with improved data
quality, we could better constrain the model parameters now.

\subsection{Modeling Total Transmission of the GV}
We derived the best-fit model for the stainless steel mesh $T_\mathrm{mesh} (E)$ in
Section~\ref{sec:model_mesh} and the Be window $T_\mathrm{Be} (E)$ in
Section~\ref{sec:model_Be}. We finally multiply them to obtain the total transmission
curve of the GV $T (E)$ as
\begin{equation}
T(E) = T_{\mathrm{Be}}(E)\left\{f+(1-f) T_{\mathrm{mesh}}(E)\right\},
\end{equation}
where \textit{f} is fixed to the best-fit value of 0.723. The model was extrapolated to cover the full range from 10 to 40000~eV and interpolated to fill in the 0.25~eV pitch to deliver as a \textit{Resolve} CALDB product.

The ground calibration result presented here has some limitations, which need to be
verified with in-orbit observations. We will (1) evaluate any systematics caused by
differences in the measurement setup and  actual observation such as beam angles, the
detector solid angle, etc, and  (2) check the validity of the linear extrapolation of
the transmission model out of the measured energy range. We have several free parameters
in our model. Fine-tuning these parameters will be achieved by observing continuum
sources such as pulsar wind nebulae and blazers in reasonable telescope times before and
after the GV open \cite{Midooka20, Miller20}.

\subsection{Revisit of Crab Spectrum Obtained with Hitomi SXS}
We have so far discussed the transmission modeling in the \textit{Resolve} GV. Our
approach in the GV modeling can be applied to the ground measurement data taken for the
SXS and be compared to the observed data with the SXS in orbit. For the stainless
steel mesh, no measurements were done for the SXS, so we fixed the opening fraction
$f=0.71$ based on the drawings. This is because, as seen in Figure~\ref{fig:mesh_curve},
the mesh is optically thick in the spectral fitting range 2--12~keV and the influence on
the residuals is small. For the Be window, we used the ground measurement results of the
SXS FS model \cite{Yoshida17} obtained in the same setup with Section~\ref{sec:E_Be} and
applied our modeling in Section~\ref{sec:model_Be}. We derived the best-fit Be
photoelectric absorption thickness to be 272.5~\textmu m\xspace.

We generated a Crab spectrum based on the IACHEC model of the Crab nebula established
for cross-calibration purposes\cite{Weisskopf10}. The model is a simple power-law of a
photon index 2.10 and a normalization of 8.70 attenuated by an interstellar absorption
of an equivalent hydrogen column density of $0.42\times10^{22}$~cm$^{-2}$. We
convolved the model with the redistribution matrix function (RMF) of the SXS CALDB and
an ancillary response file (ARF) and compared it to the observed spectrum. We modified
the ARFs by allowing some deviations in the Be window
thickness. This is because that the FM Be window used in orbit is different from the
FS Be window used in the ground measurement and that the difference can be as large as
the level of spatial non-uniformity found in the \textit{Resolve} FM Be window of
$\pm$4.2\% (Section~\ref{sec:spatial}). We made a comparison with the data for the
modified ARFs of a Be thicknesses of 261, 266, 271, 276, and 281~\textmu m\xspace and computed the
$\chi^{2}$ values.

Figure~\ref{fig:Crab} shows the comparison between the data and the model in the SXS
CALDB (blue) and the best one (261~\textmu m\xspace) among our deviated thicknesses (red). We
confirmed that the best model better describes the Crab observation results. Further
improvement is hampered by other sources of calibration uncertainties such as Au M edge
features by the X-ray mirrors.

\begin{figure}[H]
 \begin{center}
  \includegraphics[width=0.8\columnwidth]{Crabfit_comparison2.pdf}
  \caption[Comparison of simulated Crab spectra with two different ARFs and the observed
  data with the SXS]{Comparison of simulated Crab spectra with two different ARFs (blue
  and red) and the observed
  data with the SXS (black). The blue and red curves respectively show the models using
  the latest SXS CALDB and the one developed in this study, where the best Be thickness
  is 261~\textmu m\xspace.}
  \label{fig:Crab}
 \end{center}
\end{figure}

\section{SUMMARY AND CONCLUSION} \label{sec:summary}
We performed the X-ray transmission measurements of the two main components of the
\textit{Resolve} GV: the Be window and the stainless steel mesh. The Be window of the
\textit{Resolve} FM was measured in the 2.0--12.0~keV using the synchrotron facility at
KEK PF. The FS stainless steel mesh of the SXS was measured at six discrete energies
using the X-ray beamline at ISAS. The stainless steel mesh was measured for the first
time in this study.

We developed the X-ray transmission model of the GV. For the stainless steel mesh, the
wire cross section was changed from a square in the SXS model to a circle. The mesh
opening fraction was derived by fitting the data rather than being fixed to the design
value. For the Be window, the following factors were modeled: (1) the photoelectric
absorption and incoherent scattering of Be, (2) the photoelectric absorption of
impurities (Cr, Mn, Fe, Ni, and Cu), and (3) the coherent scattering of Be. The physical
interpretation of the transmission discontinuity caused by the Bragg diffraction was
presented and incorporated into the model. We also measured a spatial non-uniformity in
the Be window transmission. All these measurements and models satisfy the calibration
requirements.

In addition to constructing the transmission model of the \textit{Resolve} GV, we
revisited the models of the SXS GV based on the approach developed in this work. We
performed the reanalysis of the Crab spectrum observed with the SXS when the GV was
closed using the ground measurement data and revising the model parameters. We found
that the new model better explains the Crab observation results.

In this study, we performed an elaborated X-ray transmission modeling of the Be window
considering the beryllium's poly-crystal structure. We also showed utility of measuring
X-ray diffraction pattern in addition to the X-ray transmission. We expect that the
presented result will serve as a useful reference for all future X-ray spectrometers
with a Be window.

\clearpage


\section{Introduction}

It is commonly known that stable bright and dark solitons exist in
one-dimensional uniform self-focusing and self-defocusing nonlinear media,
respectively \cite{KivsharAgrawal,DS}. However, the situation may be totally
different in higher dimensions, where bright solitons are destabilized by
the collapse \cit
{LB1998,STA,Malomed2005,Chen2012,review2016,review2019_1,review2019_2,SREV}.
One of solutions to this problem is offered by the use of spatially
inhomogeneous settings \cite{LL1,LL2,LL3,LL4,LL5,WGA}. First, linear
periodic potentials \cite{Brazhnyi,Morsch,Vysloukh,Vysloukh2,Yang,LPP},
which are induced by optical lattices in Bose-Einstein condensates (BECs),
or photonic lattices in optical waveguides, help to stabilize various types
of self-trapped modes, including fundamental \cite{LL6}, dipole \cite{LL7},
and vortex \cite{Baizakov,LL2,Baizakov2,LL8} solitons. Another option is
the use of periodic spatial modulation of the local nonlinearity strength,
alias nonlinear lattices \cite{NLREV}. Many types of solitons are maintained
in nonlinear lattices \cite{NL1,NL2,NL3}, \textit{viz}., fundamental ones,
as well as dipole, tripole, and vortex solitons. In addition, combinations
of linear and nonlinear lattices can be used to stabilize various species of
solitons \cite{CLN1,Sakaguchi,CLN2}. Inhomogeneous self-defocusing
nonlinearity \cite{SDN1,SDN2,Dror} is another setting in which many species
of stable self-trapped modes have been predicted, such as multipole \cit
{SDN1,SDN2} and vortex \cite{SDN3,Tian,Wu,SDN4} solitons, vortex clusters
\cite{SDN5,SDN6}, skyrmions \cite{SDN7}, soliton gyroscopes \cite{SDN8}, and
flat-top solitons \cite{SDN9,SDN10}. Inhomogeneous nonlinearities can be
implemented by means of creation of a spatially nonuniform distribution of
dopants in photorefractive optical media \cite{FRM}. In BEC, similar setups
can be created too, using the Feshbach resonance controlled by spatially
nonuniform external fields \cite{FR1,FR2,FR3}.

In this work, we aim to address the possibility of using nonlinear lattices
as an ingredient of the fractional Schr\"{o}dinger equation (FSE), which was
introduced by Laskin two decades ago in the context of quantum mechanics
\cite{Lask1,Lask2,Lask3,EXP3}. The main parameter that defines the
fractionality of the FSE is the L\'{e}vy index (LI), $\alpha $, see Eq. (\re
{NLFSE}) below. The experimental implementation of the FSE in
condensed-matter \cite{EXP1,EXP2} and optical \cite{EXP4} setups, where
nonlinearity is a natural feature, has drawn interest to the possibility of
existence of solitons in fractional dimensions \cite{Frac1,LiRRP,Frac2,Frac3
. In particular, \textquotedblleft accessible solitons" \cite{Frac4,Frac5}
and self-trapped states of vectorial \cite{Frac6}, gap \cite{Frac7},
nonlocal \cite{Frac8}, vortical \cite{Frac9}, and multi-peak types \cit
{Frac10} have been predicted in FSE models, as well as soliton clusters \cit
{Frac11,Frac11b}, symmetry breaking of solitons \cite{Frac12,Frac13}, coupled
solitons \cite{Frac13b} and dissipative solitons in a fractional
complex-Ginzburg-Landau model \cite{Frac14}. In the case of the ubiquitous
cubic (Kerr) self-focusing, the solitons are unstable at $\alpha \leq 1$,
as the combination of such values with the Kerr nonlinearity gives rise to
the collapse.

The objective of this work is to identify domains of existence and stability
of fundamental and multipole (dipole and tripole) solitons in FSE with a
nonlinear lattice combining the cubic self-focusing and quintic defocusing
(vortices produced by the FSE including the uniform cubic-quintic (CQ)
nonlinearity were recently addressed in Ref. \cite{Frac9}). While
fundamental solitons produced by this model are not strongly different from
their counterparts in the one with the uniform nonlinearity, it opens a way
to create stable multipoles, i.e., complexes of solitons with opposite signs
placed at adjacent sites of the lattice.

The rest of the paper is organized as follows. In Sec. \ref{sec2} we present
the model and the method used for the linear-stability analysis. Then, in
Sec. \ref{sec3}, we summarize extensive numerical results that make it
possible to outline the stability domains for the fundamental solitons and
multipoles, \textit{viz}., dipole and tripole bound states. These results,
based on the computation of eigenvalues for small perturbations, are
corroborated by direct simulations of perturbed propagation dynamics of the
solitons. The paper is concluded in Sec. \ref{sec4}.

\section{The model}

\label{sec2} The FSE modeling the propagation of a light beam with field
amplitude $E$ under the action of the fractional diffraction in the medium
with the CQ\ spatially modulated nonlinearity is written in the normalized
form:
\begin{equation}
i\frac{\partial E}{\partial z}=\frac{1}{2}\left( -\frac{\partial ^{2}}
\partial x^{2}}\right) ^{\alpha /2}E-\mathrm{g}(x)\left\vert E\right\vert
^{2}E+\xi (x)\left\vert E\right\vert ^{4}E,  \label{NLFSE}
\end{equation
where $z$ is the propagation distance, while the fractional diffraction
operator with LI $\alpha $ is defined by the integral expression \cit
{Lask1,Zhong}
\begin{equation}
\begin{split}
\left( -\frac{\partial ^{2}}{\partial x^{2}}\right) ^{\alpha /2}E=& \frac{1}
2\pi }\int_{-\infty }^{+\infty }|p|^{\alpha }dp \\
\times & \int_{-\infty }^{+\infty }d\xi \exp \left[ ip\left( x-\xi \right)
\right] E(\xi ),
\end{split
\end{equation
the limit of $\alpha =2$ corresponding to the usual second derivative. In
the case of the cubic-only self-focusing nonlinearity, LI is limited by
values $\alpha >1$, as the collapse occurs at $\alpha \leq 1$, as mentioned
above. However, the inclusion of the quintic term with the defocusing sign
suppresses the collapse and makes it possible to consider values $\alpha <1$
as well \cite{Frac9}, although the solutions tend to be unstable in the
latter case, see Fig. \ref{fig2}(b) below. Note that, in the case of
constant nonlinearity coefficients $\mathrm{g}$ and $\xi $, Eq. (\ref{NLFSE
) admits continuous-wave (CW) solutions with an arbitrary real wavenumber $k$
and arbitrary amplitude $E_{0}$
\begin{eqnarray}
E_{\mathrm{CW}} &=&E_{0}\exp \left( ib_{\mathrm{CW}}z+ikx\right) ,  \notag \\
b_{\mathrm{CW}} &=&-\frac{1}{2}|k|^{\alpha }+\mathrm{g}E_{0}^{2}-\xi
E_{0}^{4}.  \label{CW}
\end{eqnarray}

The nonlinear lattice is introduced, replacing the constant coefficients by
the spatially-modulated ones:
\begin{equation}
\mathrm{g}(x)=\mathrm{cos}^{2}x,\,\xi (x)=\mathrm{sin}^{2}x.  \label{CQNL}
\end{equation
Here, the lattice period and amplitudes are fixed, severally, to be $\pi $
and $1$, by means of rescaling. Stationary solutions with real propagation
constant $b$ are sought for as $E(x,z)=U(x)\mathrm{exp}(ibz)$ (cf. the CW
solution (\ref{CW})), where the stationary field, $U(x)$, satisfies the
following equation:
\begin{equation}
-bU=\frac{1}{2}\left( -\frac{\partial ^{2}}{\partial x^{2}}\right) ^{\alpha
/2}U-\left( \mathrm{cos}^{2}x\right) \left\vert U\right\vert ^{2}U+\left(
\mathrm{sin}^{2}x\right) \left\vert U\right\vert ^{4}U.  \label{NLFSES}
\end{equation
Soliton solutions to Eq. (\pageref*{NLFSES}) are characterized by their
integral power,
\begin{equation}
P=\int_{-\infty }^{+\infty }\left\vert U(x)\right\vert ^{2}dx.  \label{POWER}
\end{equation
Equation (\ref{NLFSES}) is solved below by dint of the modified
squared-operator method \cite{Lakoba}.

The linear-stability analysis will be performed for a perturbed solution
taken as
\begin{equation}
E=[U(x)+p(x)\mathrm{exp}(\lambda z)+q^{\ast }(x)\mathrm{exp}(\lambda ^{\ast
}z)]\mathrm{exp}(ibz),  \label{PERB}
\end{equation
where $p(x)$ and $q^{\ast }(x)$ represent the eigenmode of small
perturbations corresponding to eigenvalue $\lambda $ (the instability growth
rate). The eigenmodes are solutions of the linearized equations produced by
the substitution of ansatz (\ref{PERB}) in Eq. (\ref{NLFSE}):
\begin{gather}
i\lambda p=+\frac{1}{2}\left( -\frac{\partial ^{2}}{\partial x^{2}}\right)
^{\alpha /2}p+bp-\left( \mathrm{cos}^{2}x\right) U^{2}(2p+q)  \notag
\label{LAS} \\
+\left( \mathrm{sin}^{2}x\right) U^{4}(3p+2q),  \notag \\
\\
i\lambda q=-\frac{1}{2}\left( -\frac{\partial ^{2}}{\partial x^{2}}\right)
^{\alpha /2}q-bq+\left( \mathrm{cos}^{2}x\right) U^{2}(2q+p)  \notag \\
-\left( \mathrm{sin}^{2}x\right) U^{4}(3q+2p).  \notag
\end{gather
The solitons are stable if $\mathrm{Re}(\lambda )=0$ for all eigenvalues.
The so predicted (in)stability will be verified through direct numerical
simulations of the full evolution equation (\ref{NLFSE}), performed by means
of the split-step Fourier method.

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig1.eps}
\end{center}
\caption{Profiles of fundamental solitons: (a) with different values of LI (
\'{e}vy index) $\protect\alpha $, for propagation constant $b=1$; (b) with
different values of $b$, for $\protect\alpha =1.5$. Note that $\protec
\alpha =2$ corresponds to the usual local (non-fractional) Schr\"{o}dinger
equation.}
\label{fig1}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig2.eps}
\end{center}
\caption{Soliton's power $P$ (a) and the maximum real part of eigenvalues
(instability growth rate) $\protect\lambda _{\mathrm{Rm}}$ (b) versus
\protect\alpha $ for fundamental solitons at $b=2$. In panel (b), the
vertical arrow indicates the boundary of the stabillity region, with very
weak instability to the left of it. Power $P$ (c) and $\protect\lambda _
\mathrm{Rm}}$ (d) versus $b$ for fundamental solitons at $\protect\alpha =1.5
$. The arrow in panel (d) indicates the respective boundary of the stability
region. Panels (a) and (c) include the best-fit approximation (red dashed
lines) provided by Eq. (\protect\ref{scaling}), with $\mathrm{const}=9$ and
12$, respectively. Panels (b) and (d) include the imaginary part, $\protec
\lambda _{\mathrm{Im}}$ (gray dashed lines), corresponding to $\protec
\lambda _{\mathrm{Rm}}$. Simulations of perturbed evolution of the
fundamental solitons marked by labels A1--A4 are displayed in Figs. \protect
\ref{fig3}(a-d), respectively.}
\label{fig2}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig3.eps}
\end{center}
\caption{The propagation of perturbed fundamental solitons: (a) a stable
soliton with $\protect\alpha =1.6$, $b=2$; (b) an unstable onen with
\protect\alpha =1.1$, $b=2$; (c) a stable soliton with $\protect\alpha =1.5
, $b=1$; (d) an unstable one with $\protect\alpha =1.5$, $b=4.5$.}
\label{fig3}
\end{figure}

\section{Numerical results}

\label{sec3} In this Section, we report numerical findings for solitons in
the FSE with the CQ lattice, including their profiles, stability domains,
and perturbed propagation. The Section is divided in two parts, which
address, separately, fundamental and multipole (dipole and tripole) solitons.

\subsection{Fundamental solitons}

Typical profiles of fundamental stationary solitons are presented in Fig.
\ref{fig1}. Naturally, the maximum of the local intensity is located at $x=0
, where the self-focusing and defocusing coefficients in Eq. (\ref{NLFSES})
attain, respectively, a maximum and minimum (zero). Accordingly, the soliton
is stabilized by being placed at the bottom of an effective nonlinear
potential well induced by the spatially-modulated CQ nonlinearity. In fact,
Fig. \ref{fig1} demonstrates that the fundamental solitons are completely
confined to a spatial interval $|x|<3\pi /4$, hence Eq. (\ref{NLFSES})
suggests that the shape of the solitons is chiefly determined by the cubic
term. Then, taking, for an approximate estimate, a constant value of the
respective coefficient, $\mathrm{g}$ (see Eq. (\ref{NLFSE})), Eq. (\re
{NLFSES}) predicts a scaling relation between the propagation constant,
amplitude, and width ($W$): $b\sim W^{-\alpha }\sim A^{2}$, which determines
the approximate scaling of the integral power,
\begin{equation}
P\sim A^{2}W\sim \left( \mathrm{const}\cdot b\right) ^{1-1/\alpha }.
\label{scaling}
\end{equation}

As LI $\alpha $ increases, the amplitude $A$ of the fundamental soliton
decreases, while its width increases, as seen in Fig. \ref{fig1}(a). In this
figure, the amplitudes are $A(\alpha =0.8)=2.34$, $A(\alpha =1.1)=1.86$, and
$A(\alpha =2)=1.73$. On the other hand, it is natural that the amplitude
increases and width decreases with the growth of propagation constant $b$,
as shown in Fig. \ref{fig1}(b), in agreement with the above-mentioned
scaling relations. In the latter figure, the amplitudes of the fundamental
soliton are $A(b=1)=1.73$ and $A(b=2)=2.36$.

Next, we present the characteristics of families of fundamental solitons. First,
in Fig. \ref{fig2}(a) we display the relation between the soliton's power $P
, calculated as per Eq. (\ref{POWER}), and LI $\alpha $, for a fixed value
of the propagation constant. It is seen that the dependence $P(\alpha )$ is
a nearly linear one. In fact, this dependence may be approximately explained
by Eq. (\ref{scaling}), as shown by the fitting line in Fig. \ref{fig2}(a).

The stability is determined by the dependence of the maximum real part of
eigenvalues, $\lambda _{\mathrm{Rm}}$, obtained from the numerical solution
of Eq. (\ref{LAS}), on $\alpha $, as shown in Fig. \ref{fig2}(b), which
identifies the respective stability region
\begin{equation}
\alpha > \alpha _{\mathrm{crit}}^{\mathrm{(fund)}}\approx 1.45,
\label{alpha-fund}
\end{equation
as defined by condition $\lambda _{\mathrm{Rm}}=0$. As mentioned above, the
plots in Figs. \ref{fig2}(a) and (b) can be extended towards $\alpha <1$,
but in that area the fundamental solitons are definitely unstable, as
suggested by Fig. \ref{fig2}(b).

It is well known that a necessary (but, generally, not sufficient) stability
condition for bright solitons is provided by the Vakhitov-Kolokolov (VK)
criterion, which is formulated in terms of the dependence of the integral
power on the propagation constant: $dP/db>0$ \cite{Vakh,LB1998}. The relation
$P(b)$ is displayed in Fig. \ref{fig2}(c), which clearly complies with the
VK criterion at all values of $b$. The $P(b)$ dependences can be also
explained in a quasi-analytical form, making use of scaling relation (\re
{scaling}), as shown by means of fitting line in Fig. \ref{fig2}(c). In
particular, for the value of $\alpha =1.5$ adopted in Fig. \ref{fig2}(c) Eq.
(\ref{scaling}) yields $P\sim b^{1/3}$.

\begin{figure}[!tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig4.eps}
\end{center}
\caption{Typical profiles of dipole bound states: (a) with different values
of $\protect\alpha $ at $b=1$; (b) with different values of $b$ at $\protec
\alpha =1.8$. Profiles of tripoles: (c) with different values of $\protec
\alpha $ at $b=1$; (d) with different values of $b$ at $\protect\alpha =1.8
. }
\label{fig4}
\end{figure}

\begin{figure}[!tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig5.eps}
\end{center}
\caption{Dipole bound states. Integral power $P$ (a) and the maximum real
part of perturbation eigenvalues $\protect\lambda _{\mathrm{Rm}}$ (b) versus
$\protect\alpha $ at a fixed propagation constant, $b=2$. Power $P$ (c) and
\protect\lambda _{\mathrm{Rm}}$ (d) versus $b$ for at fixed LI, $\protec
\alpha =1.8$. Similar to Fig. \protect\ref{fig2}, arrows in panels (b) and
(d) indicate boundaries of the respective stability regions, the instability
being very weak to the left of the boundary in (b). A cardinal difference
from the fundamental solitons (cf. Fig. \protect\ref{fig2}(d)) is that the
stability region for the dipoles in panel (d) is bounded on the side of
small $b$. The red dashed line in (a) displays the data
for the fundamental solitons from Fig. \ref{fig2}(a), with the same
propagation constant, $b=2$, multiplied by $2$. Panels (b) and
(d) also include the imaginary part $\protect\lambda _{\mathrm{Im}}$
corresponding to $\protect\lambda _{\mathrm{Rm}}$ (gray dashed lines).
Perturbed propagation of the dipoles marked by B1--B4 are displayed in
Figs. \protect\ref{fig7}(a-d), respectively.}
\label{fig5}
\end{figure}

The stability analysis is continued by data displayed in Fig. \ref{fig2}(d),
which displays numerically exact results for the maximum instability growth
rate, $\lambda _{\mathrm{Rm}}(b)$, in the same interval of $b$ that is
shown in Fig. \ref{fig2}(c). It is seen that the stability actually takes
place in a part of the interval,
\begin{equation}
0\leq b<b_{\max }^{\mathrm{(fund)}}\approx 2,  \label{b-fund}
\end{equation
and unstable at $b>b_{\max }^{\mathrm{(fund)}}$. The instability of the
narrow solitons corresponding to large $b$ may be explained by the fact that
the effectively nonlocal fractional-diffraction operator cannot provide the
stability of very narrow modes, unlike the usual local operator
corresponding to $\alpha =2$.

Predictions for the stability produced by the computation of $\lambda _
\mathrm{Rm}}$ have been verified by direct simulations of the evolution of
the fundamental solitons performed in the framework of the full nonlinear
equation (\ref{NLFSE}), as shown in Fig. \ref{fig3}. In the simulations,
random initial perturbations at the $1\%$ amplitude level were added to the
input. In Fig. \ref{fig3} the left and right column display the stable and
unstable propagation, respectively. In the former case, the solitons with
\lambda _{\mathrm{Rm}}=0$ demonstrate completely stable evolution up to
z=1000$, which corresponds, with $A\simeq 2$ (see Fig. \ref{fig1}) to $\sim
5000$ characteristic nonlinearity lengths, $z_{\mathrm{nonlin}}\sim A^{-2}$.
On the other hand, unstable solitons with $\lambda _{\mathrm{Rm}}>0$
oscillate irregularly in Fig. \ref{fig3}(b) or very slowly decay through the
emission of small-amplitude \textquotedblleft radiation" in Fig. \ref{fig3
(d). In fact, even the unstable solitons may be categorized as objects that
will seem sufficiently robust ones in the experiment, as the instability
development does not destroy them. Additional simulations demonstrate that
the instability, if any, keeps the same mild character at other values of
the parameters.

To conclude this subsection, it is relevant to mention that the simplified
model, without the quintic term and a constant coefficient in front of the
cubic one, provide similar results for the shape and stability of the
fundamental solitons. Howvere, such a model is not able to produce dipole
and tripole modes presented below.

\subsection{Dipoles and tripoles}

As mentioned above, the nonlinear lattice offers the possibility to create
multipole bound states. Typical profiles of dipole and tripole complexes,
with opposite signs of adjacent fundamental solitons that build the
dipoles, are presented, severally, in Figs. \ref{fig4}(a,b) and (c,d). The
dipoles are composed, roughly speaking, of two fundamental solitons placed
at adjacent local minima of the potential energy, $x=0$ and $\pi $ (see Eq.
\ref{CQNL})). Dipoles of this type are usually called intersite-centered
(alias densely-packed) ones, without an empty site in the middle, which
would be seen in \textit{onsite-centered} modes. They may be stable due to
the balance between repulsion between the constituents with opposite signs
and pinning to the underlying lattice \cite{Kapitula}. On the other hand,
lattice-based complexes built of adjacent fundamental solitons with
identical signs are unstable in accordance with the general analysis \cit
{Kapitula}, therefore they are not reported here. The tripoles are also
built as tightly-packed complexes.

According to Figs. \ref{fig4}(a) and (c), amplitudes and widths of both the
dipole and tripole families gradually decrease and increase, respectively,
with the growth of LI $\alpha $, similar to what is observed above for the
fundamental solitons, cf. Fig. \ref{fig1}(a). For example, the dipole's
amplitudes are $A(\alpha =0.8)=2.36$, $A(\alpha =1.1)=1.9$ and $A(\alpha
=2)=1.76$. Further, it is seen in Figs. \ref{fig4}(b) and (d) that the
amplitude and width, respectively, gradually increase and decrease with the
growth of propagation constant $b$, also similar to properties of the
fundamental solitons shown in Fig. \ref{fig1}(b). Accordingly, typical
values of the dipole's amplitude are $A(b=1)=1.75$, and $A(b=2)=2.36$.

\begin{figure}[!tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig6.eps}
\end{center}
\caption{Characteristics of the tripole complexes. Integral power $P$ (a)
and the maximum real part of perturbation eigenvalues $\protect\lambda _
\mathrm{Rm}}$ (b) versus $\protect\alpha $ at a fixed propagation constant,
b=2$. Power $P$ (c) and $\protect\lambda _{\mathrm{Rm}}$ (d) versus $b$ at
fixed LI, $\protect\alpha =1.8$. Similar to Fig. \protect\ref{fig5}, arrows
in (b) and (d) indicate boundaries of the respective stability regions.
The red dashed line in (a) displays the data
for the fundamental solitons from Fig. \ref{fig2}(a), with the same
propagation constant, $b=2$, multiplied by $3$.
Panels (b) and (d) also include the  imaginary
part $\protect\lambda _{\mathrm{Im}}$ (gray dashed lines) corresponding to
\protect\lambda _{\mathrm{Rm}}$. Simulations of the perturbed propagations
of the tripole states marked by C1--C4 are displayed in Figs. \protect\re
{fig8}(a-d), respectively.}
\label{fig6}
\end{figure}

\begin{figure}[!tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig7.eps}
\end{center}
\caption{The propagation of perturbed dipole solitons: (a) a stable soliton
at $\protect\alpha =1.7$, $b=2$; (b) an unstable one at $\protect\alpha =1.3
, $b=2$; (c) a stable soliton at $\protect\alpha =1.8$, $b=2$; (d) an
unstable one at $\protect\alpha =1.8$, $b=0.5$.}
\label{fig7}
\end{figure}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig8.eps}
\end{center}
\caption{The propagation of perturbed tripole solitons: (a) a stable one at
\protect\alpha =1.9$, $b=2$; (b) an unstable soliton at $\protect\alpha =1.4
, $b=2$; (c) a stable one at $\protect\alpha =1.8$, $b=1.5$; (d) an unstable
soliton at $\protect\alpha =1.8$, $b=0.5$.}
\label{fig8}
\end{figure}

Proceeding to characteristics of the family of dipole states, the respective
relation between the integral power $P$, as defined by Eq. (\ref{POWER}),
and LI $\alpha $ is displayed in Fig. \ref{fig5}(a), which shows a linear
growth of $P$ with $\alpha $. In fact, this $P(\alpha )$ curve is
nearly identical to its counterpart for the fundamental solitons, displayed
for the same fixed value of the propagation constant, $b=2$, in Fig. \ref{fig2}(a),
multiplied by $2$. Similar to the case of the fundamental
solitons (cf. Fig. \ref{fig2}(a)), this dependence can be explained,
although with poorer accuracy, by scaling relation (\ref{scaling}), as shown
by the fitting curve in \ref{fig5}(a). Next, the maximum real part of the
perturbation eigenvalues for the dipoles, $\lambda _{\mathrm{Rm}}$, produced
by the numerical solution of Eq. (\ref{LAS}), is displayed, as a function of
$\alpha $ in Fig. \ref{fig5}(b), demonstrating the stability region
\begin{equation}
\alpha >\alpha _{\mathrm{crit}}^{\mathrm{(dip)}}\approx 1.6,
\label{alpha-dip}
\end{equation
which is slightly narrower than its counterpart for the fundamental solitons
given by Eq. (\ref{alpha-fund}).

Another characteristic of the dipole family is $P(b)$ dependence, displayed
in versus $b$ in Fig. \ref{fig5}(c), which satisfies the VK criterion. Note
that the $P(b)$ dependence is still very accurately approximated by scaling
relation (\ref{scaling}), as shown by the fitting curve in the figure. The
corresponding dependence of the instability growth rate, $\lambda _{\mathrm
Rm}}$, on $b$ is displayed in Fig. \ref{fig5}(d). This plot demonstrates a
drastic difference from its counterpart for fundamental solitons in Fig. \re
{fig2}(d), as the stability interval for the dipoles is a \emph{window},
bounded both from above and from below,
\begin{equation}
b_{\min }^{\mathrm{(dip)}}\approx 0.9<b<b_{\max }^{\mathrm{(dip)}}\approx
2.9,  \label{b-dip}
\end{equation
cf. Eq. (\ref{b-fund}). The instability of broad dipoles corresponding to
small values of the propagation constant, $b<b_{\min }^{\mathrm{(dip)}}$, is
explained by the fact that the strong overlap between broad constituent
solitons that build the dipole leads to strong interaction between them,
which produces a destabilizing effect.

For the tripole complexes the relation between the integral power, $P$, and
LI $\alpha $ is presented in Fig. \ref{fig6}(a). This nearly linear
$P(\alpha )$ dependence is very close to its counterpart for the fundamental
solitons, displayed for the same fixed $b=2$, in Fig. \ref{fig2}(a),
multiplied by $3$ (again satisfying the VK criterion), similar to the above
notice that the $P(\alpha )$ curve in Fig. \ref{fig5}(a) is closely
approximated by its counterpart for the fundamental solitons from Fig. \re
{fig2}(a) multiplied by $3$. Next, in Fig. \re
{fig6}(b) the maximum real part of eigenvalues $\lambda _{\mathrm{Rm}}$,
produced by the numerical solution of Eq. (\ref{LAS}), is displayed as a
function of $\alpha $, the respective stability region bein
\begin{equation}
\alpha >\alpha _{\mathrm{crit}}^{\mathrm{(trip)}}\approx 1.8,
\label{alpha-trip}
\end{equation
which is narrower than its counterpart (\ref{alpha-dip}) for the dipoles.
The same eigenvalue is plotted versus the propagation constant
$b$ in Fig. \ref{fig6} (d). In terms of $b$, the stability window is
\begin{equation}
b_{\min }^{\mathrm{(trip)}}\approx 1<b<b_{\max }^{\mathrm{(trip)}}\approx
1.8,  \label{b-trip}
\end{equation
which is considerably smaller than a similar one (\ref{b-dip}) for the
dipoles. It is not surprising that the more complex bound states (tripoles)
is more fragile than simpler ones, arranged as dipoles.

Direct simulations of the perturbed propagation of dipole states are
presented in Fig. \ref{fig7}, in which the left and right columns refer,
respectively, to the stable and unstable propagation. According to the left
column, stable dipoles keep their shapes in the course of the very long
propagation, up to $z=1000$. On the other hand, those dipoles that are
predicted to be unstable by means of the linear-stability analysis feature
irregular oscillations and gradual decay through the emission of
small-amplitude waves, in the right column of Fig. \ref{fig7}. Similarly, in
Fig. \ref{fig8}, stable tripoles keep their shapes in the left column, while
their unstable counterparts also demonstrate irregular oscillations and
gradual delay.

\section{Conclusion}

\label{sec4} In this work, a one-dimensional model combining the fractional
diffraction and nonlinear CQ (cubic-quintic) lattice is introduced, and
families of fundamental and multipole (dipole and tripole) solitons in it
are reported. The shapes and stability of these states are established by
means of the systematic numerical investigation. Some results are explained
with the help of the analytical approximation, which is based on the
consideration of scaling properties of the underlying equation. The solitons
are stable above a critical value of the LI (L\'{e}vy index), $\alpha $. In
terms of the propagation constant, $b$, the fundamental solitons are stable
below a critical value, while the dipole and tripole complexes are stable in
intervals of $b$ limited from above and below. The stability areas of all
the states, predicted by the calculation of the perturbation eigenvalues,
are corroborated by direct simulations of the perturbed evolution. In the
simulations, the instability of some states seems weak, which suggests that
they may also represent relevant objects for experimental observation.
Naturally, the stability regions, in terms of both $\alpha $ and $b$, are
widest for the fundamental solitons, and narrowest for the most complex
states, \textit{viz., the tripoles.}

As an extension of the present analysis, it may be interesting to consider
mobility of the solitons on top of the underlying nonlinear lattice, and
collisions between moving ones. A straightforward extension may also address
bound states of more than three solitons, and dynamics of multi-soliton
arrays pinned to the nonlinear CQ lattice.

\section*{Declaration of Competing Interest}
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

\section*{Funding}
National Major Instruments and Equipment Development Project of National Natural Science Foundation of China (No. 61827815); National Natural Science Foundation of China (No. 62075138); Science and Technology Project of Shenzhen (JCYJ20190808121817100, JCYJ20190808164007485); Israel Science Foundation (grant No. 1286/17).


\section{Introduction}
\label{sec:intro}

Fundamental frequency (F0) is an essential acoustic feature in human speech communication and human-human interactions. As a key feature of speech prosody, F0 plays an important role in every aspect of speech communication: it conveys linguistic information (F0 helps to clarify the syntactic structure of an utterance or is used for semantic emphasis), para-linguistic information such as emotion or social attitude, and is even a part of the speaker identity through his speaking style.
Accordingly, generative F0 modeling can be extremely useful in the fields of text-to-speech \cite{wang2018style}, voice identity conversion (VC) \cite{Kaneko_2018_CGANVC}, and  expressive voice conversion \cite{Luo_2017}, by allowing a direct and parametric control of F0 to manipulate the expressivity of a voice (such as speaking style or emotions).
By nature, F0 variations occur over different temporal scales each associated with specific functions, ranging from micro-variations to macro-contours such as accentuation, emotions, and modalities. To cover these specificities of F0 modelling, stylization methods \cite{Teu08, Obi18} and multi-level modelling \cite{Obi11c, Yin16, Wang2017-xu} have been proposed.

Notably in VC, generative models such as Gaussian mixture models \cite{VeauxGMM} and LSTM-based Sequence-to-Sequence models \cite{robinson2019} were used to learn F0 transformations from neutral to expressive speech. Lately, various works focused on the use of Continuous Wavelet Transform (CWT), as an intermediary representation of the F0, on which Generative Adversarial Networks (GAN) models such as Dual-GAN \cite{Luo_2019}, Cycle-GAN \cite{zhou2020_CG}, VAW-GAN \cite{zhou2020_VAWG} or VA-GAN \cite{luo_chen_takiguchi_ariki_2019} are trained to learn transformations. A majority of those models is learnt on parallel data and by emotion pairs, which allows to learn a direct mapping between two different emotional versions of an utterance while preserving a fixed and controlled linguistic content.

A promising approach called CWT Adaptive Scales (CWT-AS) was proposed by Luo et al. \cite{Luo_2019}. The CWT computes a decomposition of the F0 signal over wavelet kernels which allows a representation of F0 over different temporal scales \cite{Ming15}, with various application in expressive voice conversion \cite{Ming15, Ming2016-tp, Luo_2017, Luo_2019, zhou2020_VAWG}. F0 modelling with CWT has been specified more recently upgraded with the possibility to compute the decomposition on arbitrary linguistic scales (e.g., phoneme, syllable, word, and utterance as described in \cite{Luo_2017}). An Adaptive-Scale (AS) algorithm \cite{Luo_2019} is described to select an optimal CWT representation for each pair of emotions, by selecting the scales that maximize in average the distance between the emotions in the CWT space.

From these selected scales, the CWT decomposition of the F0 contours is computed. Finally, the transformation function between each pair of emotion is learned from those representation using a Dual-GAN.
Though this approach appears promising, it suffers from two main limitations: 1) The scale selection is only based on the maximization of the distance between the emotions, but ignores their reconstruction ability of the F0 signal. This may lead to poor F0 reconstruction which in turn would degrade the quality and the naturalness of the transformation; 2) The CWT-AS decomposition of the F0 signal and the dual-GAN are optimized independently which constitutes a bottleneck for training. Consequently, the CWT decomposition may not be optimal in the sense of the dual-GAN objective.

To overcome those limitations, we propose a end-to-end architecture to learn efficiently F0 transformation between emotions. The proposed neural architecture brings together the F0 decomposition and the dual-GAN into a single network, so that the CWT decomposition is optimized in the  sense of the dual-GAN objective, and combining separation and reconstruction losses of the resulting decomposition. An application to the voice conversion of social attitudes shows that the proposed approach significantly improves the quality of the transformation by comparison with the CWT-AS approach.

\section{PROPOSED METHOD}
\label{sec:meth}

In this section, we introduce our proposal based on CWT-AS \cite{Luo_2019} and show how it differs by integrating the F0 stylization part on top of the transformation learning process which we refer to as end-to-end method for voice f0 conversion in \ref{Framework}. Our contributions' concepts and technical details are given in \ref{WKCE} and \ref{Model}.

\subsection{Framework overview}
\label{Framework}

As our proposed VC system requires parallel data, sets of utterances $X^{a}$ and $X^{b}$ respectively relative to expressivity $a$ and $b$ are considered. A pair of utterances is then sampled and F0 sequences, source $\mathbf{x}^{a}$ and target $\mathbf{x}^{b}$, are extracted. Aside from the expressivity, each utterance in a pair has the same content (linguistic content, speaker identity). The source and target F0 are given to what we called a Wavelet Kernel Convolutional Encoder (WKCE) denoted $W_{e}$. A classifier, denoted $C$, whose objective is to predict the expressivity is fed with WKCE outputs. As shown in Figure \ref{dual_gan}, these two modules must be seen as a pre-network ($pN$) for Dual-GAN ($DG$) that can be pre-trained as well as trained along with Dual-GAN forming an end-to-end system for f0 conversion.

\begin{figure*}[h]
    \centering
    \includegraphics[width=\textwidth]{Figs/framework_mod.pdf}
    \caption{End-to-end neural architecture for F0 voice conversion. On the right, F0 decomposition over 4 of the learned scales obtained for the source (red) and target (green) expressivities.}
    \label{dual_gan}
\end{figure*}

\subsection{Wavelet Kernel Convolutional Autoencoder}
\label{WKCE}

As a multiscale modelling method, CWT is entirely fitting when trying to represent both long and short-term dependencies, prosody is influenced by. As CWT can only be applied to continuous functions, a simple linear interpolation between voiced F0 segments is needed to obtain a continuous phrase-related F0 function which can then be sampled in a vector $\mathbf{x} \in [0, 1]^{T}$.

Our WKCE performs convolutions between the F0 signal $\mathbf{x}$ and a wavelet kernel based on a mother wavelet $\boldsymbol{\psi_{s}} \in \mathbb{R}^{T}$ defined for a time vector $\mathbf{t} \in \mathbb{R}^{T}$ as

\begin{equation}
\small
\boldsymbol{\psi}_{s} = \frac{2\pi^{-\frac{1}{4}}}{\sqrt{3}}(1-\left(\frac{\mathbf{t}}{s}\right)^{2})e^{-\frac{1}{2} \left(\frac{\mathbf{t}}{s}\right)^{2}}
\end{equation}

Considering a wavelet kernel depending on $N$ learnable parameters $s$ that control the width of each wavelet composing the kernel, the contribution  $\mathbf{h}_{\mathbf{x}}^{s}$ of the temporal level $s$ to the F0 signal $\mathbf{x}$ is the convolution between $\mathbf{x}$ and $\boldsymbol{\psi}_{s}$. Therefore, feed with $\mathbf{x}$, our WKCE module will output $W_{e}(\mathbf{x}) = [\mathbf{h}_{\mathbf{x}}^{s_{0}}, ..., \mathbf{h}_{\mathbf{x}}^{s_{N}}] \in \mathbb{R}^{N \times T}$. If we denote $W_{r}$ the reconstruction operation, then the reconstructed signal $\hat{\mathbf{x}}$ is given by

\begin{equation}
\small
	\hat{\mathbf{x}} = W_{r}(W_{e}(\mathbf{x})) = \frac{d_{j}\sqrt{d_{t}}}{C_{d} Y_{0}}\sum_{i=0}^{N-1} \mathbf{h}_{\mathbf{x}}^{s_{i}} + \mathbf{\bar{x}}
	\label{cwt_rec}
\end{equation}

with $\mathbf{\bar{x}}$ the average of $\mathbf{x}$, $d_{t}=1.2$, $d_{j}=0.125$, $C_{d}=3.541$ and $Y_{0} = 0.867$ (for details, see \cite{Torrence1998APG}).

If we denote $\mathbb{E}$, the mathematical expectation and consider $\mathbf{x}^{a}$ and $\mathbf{x}^{b}$ sampled from source and target distributions $P(\mathbf{x}^{a})$ and $P(\mathbf{x}^{b})$ respectively, this module can be trained for reconstruction objective with respect to $L_{1}$ loss formulated as follows

\vspace*{-\baselineskip}
\vspace{0.8mm}

\begin{equation}
    \small
    \begin{split}
    L_{rec} = &\mathbb{E}_{\mathbf{x}^{a} \sim P(\mathbf{x}^{a})}(||W_{r}(W_{e}(\mathbf{x}^{a}))) - \mathbf{x}^{a}||_{1}) +\\
            &\mathbb{E}_{\mathbf{x}^{b} \sim P(\mathbf{x}^{b})}(||W_{r}(W_{e}(\mathbf{x}^{b})) - \mathbf{x}^{b}||_{1})
    \end{split}
    \label{rec}
\end{equation}

A constraint of classification on the CWTs latent space can be added, $W_{e}$ and $C$ are trained with respect to $L_{cl}$, the cross-entropy (CE) loss between the predicted source expressivity $\hat{a} = C(W_{e}(\mathbf{x}^{a}))$ and the true value $a$ summed with the CE between $\hat{b}$ and $b$.

\begin{equation}
\small
\begin{split}
    L_{cl} = &\mathbb{E}_{\mathbf{x}^{a} \sim P(\mathbf{x}^{a})}[a*C(W_{e}(\mathbf{x}^{a}))] +                       \mathbb{E}_{\mathbf{x}^{b} \sim P(\mathbf{x}^{b})}[b*C(W_{e}(\mathbf{x}^{b}))] \\
              + &\mathbb{E}_{\mathbf{x}^{b} \sim P(\mathbf{x}^{b})}(1-a)[1 - log(C(W_{e}(\mathbf{x}^{a})))] \\
              + &\mathbb{E}_{\mathbf{x}^{a} \sim P(\mathbf{x}^{a})}(1-b)[1 - log(C(W_{e}(\mathbf{x}^{b}))))]
\end{split}
\label{adv}
\end{equation}

\subsection{Model}
\label{Model}

In this paper we focus on a specific GAN network called Dual-GAN which is capable of learning a mapping between parallel pairs of data. This network is based on two concepts: 1) Adversarial learning \cite{NIPS2014_5423}, which is to train a generative model to find a solution in a min-max game between two neural networks, called as generator $G$ and discriminator $D$. 2) Dual supervised learning \cite{Xia_2017} which is to train the models of two dual tasks simultaneously exploiting the probabilistic correlation between them to regularize the training process. Combining those breakthroughs allows to take advantage of the GAN's ability to produce realistic transformations as well as the significant improvements due to dual supervised learning.

This second point implies that both forward and inverse transformations, respectively $G_{a \rightarrow b} : (W_{e}(\mathbf{x}^{a}), z^{a}) \rightarrow \mathbf{x}^{b}$ and $G_{b \rightarrow a} : (W_{e}(\mathbf{x}^{b}, z^{b}) \rightarrow \mathbf{x}^{a}$, are learned jointly, where $z^{a}$ and $z^{b}$ are random independant noises provided in the form of dropout at each layer of $G_{a}$ and $G_{b}$. A first loss $L_{a \leftrightarrow b}$ is required to train $G_{a \rightarrow b}$, $G_{b \rightarrow a}$ and $W_{e}$.

\begin{equation}
\small
    \begin{split}
    L_{a \leftrightarrow b} = &\mathbb{E}_{(\mathbf{x}^{a}, \mathbf{x}^{b}) \sim P(\mathbf{x}^{a}, \mathbf{x}^{b})}(||W_{r}(G_{a \rightarrow b}(W_{e}(\mathbf{x}^{a}))) - \mathbf{x}^{b}||_{1}) \\
            + &\mathbb{E}_{(\mathbf{x}^{a}, \mathbf{x}^{b}) \sim P(\mathbf{x}^{a}, \mathbf{x}^{b})}(||W_{r}(G_{b \rightarrow a}(W_{e}(\mathbf{x}^{b}))) - \mathbf{x}^{a}||_{1})
    \end{split}
    \label{rec}
\end{equation}

In the same time, $D_{a}$ discriminates between converted outputs $\hat{X^{b}}$ of $G_{a \rightarrow b}$ and real samples of domain $X^{b}$, $D_{b}$ does analogously the same to complete the adversarial mechanism. The adversarial loss $L_{ADV}$ is required to train $G_{a \rightarrow b}$, $G_{b \rightarrow a}$, $D_{a}$, $D_{b}$ and $W_{e}$

\begin{equation}
\small
\begin{split}
    L_{adv} = &\mathbb{E}_{\mathbf{x}^{a} \sim P(\mathbf{x}^{a})}[D_{a}(W_{e}(\mathbf{x}^{a}))] +                      \mathbb{E}_{\mathbf{x}^{b} \sim P(\mathbf{x}^{b})}[D_{b}(W_{e}(\mathbf{x}^{b}))] \\
              + &\mathbb{E}_{\mathbf{x}^{b} \sim P(\mathbf{x}^{b})}[1 - log(D_{a}(G_{b \rightarrow a}(W_{e}(\mathbf{x}^{b}))))] \\
              + &\mathbb{E}_{\mathbf{x}^{a} \sim P(\mathbf{x}^{a})}[1 - log(D_{b}(G_{a \rightarrow b}(W_{e}(\mathbf{x}^{a}))))]
\end{split}
\label{adv}
\end{equation}

A third constraint called Dual loss is added so as to strengthen the intrinsic connection between $G_{a \rightarrow b}$ and $G_{b \rightarrow a}$, it can be understood as a regularization of the process.

\begin{equation}
\small
\begin{split}
    L_{dual} = \mathbb{E}_{(\mathbf{x}^{a}, \mathbf{x}^{b}) \sim P(\mathbf{x}^{a}, \mathbf{x}^{b})}(||&W_{e}(\mathbf{x}^{a})*G_{a \rightarrow b}(W_{e}(\mathbf{x}^{a})) \\
    - &W_{e}(\mathbf{x}^{b})*G_{b \rightarrow a}(W_{e}(\mathbf{x}^{b}))||_{1})
\end{split}
\label{adv}
\end{equation}

Therefore two final losses can be formulated for pre-Net pretraining and proper Dual-GAN training, respectively $L_{pN}$ and $L_{DG}$ with $\alpha$, $\beta$, $\lambda$ and $\gamma$ respectively weighting reconstruction, classification, transformation and dual objectives.

\vspace*{-\baselineskip}
\vspace{0.8mm}

{\footnotesize
\begin{align}
\small
L_{pN} &= \alpha L_{rec} + \beta L_{cl} \\
 L_{DG} &= \lambda L_{a \leftrightarrow b} +  L_{adv} + \gamma L_{dual}
\end{align}
}%

\section{Experiments}
\label{sec:exp}

\subsection{Dataset}

For our experiments we used the freely available speech database Att-HACK \cite{LeMoine_sp}. The database comprises 25 speakers interpreting 100 utterances in 4 social attitudes : friendly, distant, dominant and seductive, later denoted FR, DIST, DOM and SED respectively. With 3 to 5 repetitions each per attitude for a total of around 30 hours of speech, the database offers a wide variety of prosodic strategies in the expression of attitudes.
Prosodic features conveying expressivity has been shown to be speaker dependent \cite{Sisman_18}, for this reason, two speakers were selected and used independently for both training and validation : a female (F08) and a male (M07), representing almost 400 utterances each. The \textit{train/valid} split was 80/20 \% and has been done linguistically.

\subsection{Implementation details}

\subsubsection{Input pipeline}

We extracted fundamental frequency from the speech signal by using SWIPEP algorithm \cite{Cam07}. All F0 sequences are sampled to 1ms (as recommended in \cite{Luo_2019}), passed to log(F0) and a linear interpolation has been processed between voiced segments. For each pair, a mapping between syllables starting and ending times has been done to align source and target. Once pairs are aligned syllable-wise, F0 sequences are padded with zeros up to a value $T=4000$.

\subsubsection{Architecture design}

Our WKCE denoted $W_{e}$, with $N=32$ learnable scales, has been implemented as a custom layer. A constraint of growth has been added on the range of scales to ensure the continuity of the learned CWTs. The output of $W_{e}$ of shape [32, 4000] is unpadded and temporally sliced to form batches of shape [batch\_size, 32, 512].

The classifier is built using convolutional blocks composed of three convolutional layers starting from 32 up to 128 filters. Each block uses 3 × 3 convolutions with \textit{ReLU} activation, a dropout of 0.2, padding mode same and pooling operations using both strides and max pooling 2D with values 2 and 4, respectively, reducing features and time. Those blocks are followed with a flatten layer and two dense layers with respectively, 1000 units and a \textit{ReLU} activation, and 2 units and a \textit{softmax} activation.

Two configurations of our pre-network $pN$ can be distinguished : \begin{itemize}
    \item  \textit{config\_A} : $pN = \{W_{e}\}$ learns the CWT scales regarding the CWT reconstruction objective ($\alpha=1$, $\beta=0$)
    \item \textit{config\_B} : $pN = \{W_{e} + C\}$ learns the CWT scales regarding both the CWT reconstruction and the CWT related attitude classification objectives ($\alpha=10$, $\beta=1$)
\end{itemize}

The architecture of the Dual-GAN itself as well as the contribution of each module, $G_{a \rightarrow b}$, $G_{b \rightarrow a}$, $D_{a}$ and $D_{b}$, in training process are taken from \cite{Luo_2019}.

\vspace{0.8mm}

\subsubsection{Training procedure}

$W_{e}$ is being fed with batches of size 1 so as to allow the Dual-GAN to process on phrase related batches. $pN$ is pre-trained minimizing $L_{pN}$ depending on the considered configuration. CWTs are learned with respect to voiced segments by using a voicing binary mask of 1 (voiced) and 0 (unvoiced). Then $\{DG + W_{e}\}$ is trained minimizing $L_{DG}$ with $\lambda=5$ and $\gamma=15$. For \textit{config\_B}, the classification scores after pre-training were passed as sample weights for Dual-GAN training. ADAM optimizer with $0.0001$ as learning rate has been used. All codes are written in Python-Tensorflow 2.1, the baseline has been re-implemented.

\section{RESULTS AND DISCUSSIONS}
\label{sec:res}

\subsection{Objective Evaluation}
\label{sssec:obj}

Considering 4 attitudes, we evaluated 12 transformations (6 forward, 6 invert), the overall results are shown in Table \ref{RMSE}, using the RMSE between 1) the original F0 and the reconstruction obtained from its representation, 2) the converted F0 and the corresponding target F0. An example of conversion is depicted in Figure \ref{conv}.

First of all, we observe that the proposed end-to-end system for f0 voice conversion outperforms the traditional \textit{baseline}, in average over all categories, for both configurations. Secondly our most sophisticated configuration \textit{config\_B} achieves slightly better results than \textit{config\_A}.

\begin{table}[h]
\centering
\begin{tabular}{c|cc}

    & \multicolumn{2}{c}{\textbf{RMSE (Hz)}} \\

    \hline
    \textbf{Models} & \textbf{Reconstruction} & \textbf{Transformation} \\

    \hline
    \textit{baseline} & 17.32 & 21.71 \\
    \textit{config\_A} & 9.16 & 19.15 \\
    \textit{config\_B} & 13.4 & \textbf{18.83} \\

\end{tabular}
\caption{A comparison of the RMSE results of the \textit{baseline}, \textit{config\_A} and \textit{config\_B} for reconstruction and transformation}
\label{RMSE}
\end{table}

\vspace{-5.5mm}

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.49\textwidth]{Figs/F0_conv_.pdf}
    \caption{Example of F0 transformation from distant to dominant for speaker M07 for \textit{baseline} and ours \textit{config\_A} and \textit{config\_B}.}
    \label{conv}
  \end{figure}

\vspace*{-\baselineskip}

\subsection{Subjective Evaluation}
\label{sssec:sub}

We further conducted a listening experiment to compare the two proposed configurations with the baseline CWT-AS in terms of attitude similarity. We perform XAB test to assess the emotion similarity by asking listeners to choose between two converted utterances (\textit{baseline} and one of our configurations) the one which sounds more similar to the original target. The overall results are reported in Figure \ref{XAB}. Our proposed system in \textit{config\_B} outperforms the baseline for each of the considered transformations. The results for \textit{config\_A}, slightly lower, still outperforms the \textit{baseline}.

\vspace*{-\baselineskip}

  \begin{figure}[h!]
    \centering
    \includegraphics[width=0.35\textwidth]{Figs/overall_trans.pdf}
    \caption{The XAB preference results with 95\% confidence interval between the \textit{baseline} and ours \textit{config\_A} and \textit{config\_B} regarding attitude similarity.}
    \label{XAB}
  \end{figure}

\vspace*{-\baselineskip}

\subsection{Scales distributions}

This part provides an a posteriori comparison of the F0 scales distribution as obtained by the description presented in parts \ref{sssec:obj} and \ref{sssec:sub}. Each transformation (forward and backward) between a pair of attitudes is associated with a set of temporal scales that are used to compute the CWT representations used for the conversion. Consequently, each transformation can be described by a distribution of the temporal scales that are used to convert the F0 optimally.  Figure \ref{dist} presents the distribution selected by the baseline CWT-AS algorithm and learned by our proposed contribution, as obtained for speaker F08 for the six pairs of attitudes.

The best performance being obtained with the \textit{config\_B} of the proposed contribution, we further investigated by comparison of the F0 scales distribution.
First, one can clearly observe that the temporal distribution of the \textit{config\_B} is wider than the others, the transformation covering a wide range of temporal scales from the micro variations over the phonemes to the global contours of the sentence. Additionally, the distributions associated with the \textit{baseline} and the \textit{config\_A} appear mostly independent with respect to the transformation pair, while the distribution associated with the \textit{config\_B} tend to be more varied depending on the transformation pair. This suggests that the\textit{config\_B} may adapt more efficiently to the specificities of each pair

\vspace*{-\baselineskip}
\vspace{1.5mm}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{Figs/scales_distrib.png}
    \caption{CWT scale distributions for the three considered algorithms: \textit{baseline}, \textit{config\_A} and \textit{config\_B}. For interpretation, P, SY, W and SE markers respectivelly denote average duration of the phone, syllable, word and sentence linguistic units.}
    \label{dist}
\end{figure}

\vspace{-4.5mm}
\section{Conclusion}
\label{sec:con}

\vspace{-1.5mm}

In this paper, we propose a end-to-end framework for F0 transformation in the context of expressive voice conversion, bringing together the F0 decomposition in different temporal levels and its transformation in a single network. Both objective and subjective evaluations showed our method can achieve better performance than the baseline. We aim at generalizing for multi-speaker F0 conversion and to avoid pair-learning by building an expressive embedding. An online page featuring conversion examples is available at \href{http://recherche.ircam.fr/anasyn/VC_demo/index.html}{http://recherche.ircam.fr/anasyn/VC\_demo/index.html}.

\section*{Acknowledgements}
\label{sec:akn}
This research is supported by the MoVe project: ‘‘MOdelling of speech attitudes and application to an expressive conversationnal agent'', and funded by the Paris Region Ph2D grant.

\bibliographystyle{IEEEbib}