11institutetext: 1School of Information Science and Technology,
ShanghaiTech University, Shanghai, China
2School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai,
China
3Institute of Brain-Intelligence Technology, Zhangjiang Laboratory
4Shanghai Engineering Research Center of Intelligent Vision and Imaging,
ShanghaiTech University, Shanghai, China
11email: {yujingyi<EMAIL_ADDRESS>
# IREM: High-Resolution Magnetic Resonance Image Reconstruction via Implicit
Neural Representation

Qing Wu1 Yuwei Li1 Lan Xu1 Ruiming Feng2 Hongjiang Wei2 Qing Yang3 Boliang Yu1
Xiaozhao Liu1 Jingyi Yu1,(🖂) and Yuyao Zhang1,4,(🖂)  Yuyao Zhang and Jingyi Yu
are co-corresponding author.

###### Abstract

For collecting high-quality high-resolution (HR) MR image, we propose a novel
image reconstruction network named IREM, which is trained on multiple low-
resolution (LR) MR images and achieve an arbitrary up-sampling rate for HR
image reconstruction. In this work, we suppose the desired HR image as an
implicit continuous function of the 3D image spatial coordinate, and the
thick-slice LR images as several sparse discrete samplings of this function.
Then the super-resolution (SR) task is to learn the continuous volumetric
function from a limited observations using a fully-connected neural network
combined with Fourier feature positional encoding. By simply minimizing the
error between the network prediction and the acquired LR image intensity
across each imaging plane, IREM is trained to represent a continuous model of
the observed tissue anatomy. Experimental results indicate that IREM succeeds
in representing high frequency image feature, and in real scene data
collection, IREM reduces scan time and achieves high-quality high-resolution
MR imaging in terms of SNR and local image detail.

###### Keywords:

MRI High Resolution Super Resolution Implicit Neural Representation

##  1 Introduction

High Resolution (HR) medical images provide rich structural details to
facilitate early and accurate diagnosis. In MRI, unpredictable patient
movements result in difficulties to collect artifact-free high-resolution 3D
images at sub-millimeter within a single scan. Thus, multiple anisotropic low-
resolution (LR) MR scans with thick slice are always conducted at different
scan orientations with a relatively short scan time to eliminate motion
artifacts. MRI Super-Resolution (SR) techniques are then used to enhance the
image spatial resolution from one or more LR images. According to the number
of processed images, related works can be divided into two categories:
multiple image super-resolution (MISR) [19, 2] or single image super-
resolution (SISR) [17].

MRI SR techniques were originally proposed as MISR algorithms [14]. Multiple
thick-sliced LR images of the same subject were acquired to reconstruct an HR
image. Gholipour et al. [5] developed a model-based MISR technique (always
known as SRR) that enabled isotropic volumetric image reconstruction from
arbitrarily oriented slice acquisitions. SRR achieved great success and was
broadly applied for motion-prone MR acquisitions (i.e., fetus brain
construction [4]). Jia et al. [9] proposed a Sparse Representation based SR
approach using only the multiple LR scans to build the over-complete
dictionary. Previous studies [5, 21, 22, 24] showed that for MISR methods, in
a given acquisition time, when the alignment among multiple frames is known or
correctly estimated, the SR reconstructed image presented higher SNR (thick-
slice improves SNR in single voxel) when comparing with the isotropic image
acquired with similar spatial resolution. But it is challenging for MISR
methods to build submillimeter-level image detail, such as in brain
cerebellum.

On the other hand, SISR is attracting attentions due to its advantage of
shorter scan time to further reduce the motion artifacts. Recently, the
performance of SISR methods were significantly improved by Deep convolution
neural network (CNN) because of its non-linearity to imitate transformation
between LR and HR images. Various network architectures for MRI SISR have been
proposed [11, 12, 3, 15, 1]. However, there have to exist a large HR image
dataset as ground truth to train the SR network and the learned model lacks
generalization for the input images with different qualities, e.g., different
spatial resolutions, SNRs. As a result, this kind of methods cannot handle SR
tasks when no HR ground truth is available (i.e., fetus brain reconstruction).

In this paper, we propose IREM, which is short for Implicit REpresentation for
MRI. Our work intends to further improve MISR reconstruction performance via
deep neural network. Instead of CNN architecture, we train a fully-connected
deep network (MLP) to learn and represent the desired HR image. Specifically,
we suppose the HR image as an implicit continuous function of the 3D image
spatial coordinate, and the 2D thick-slice LR images as several sparse
discrete samplings of this function. Then the SR task is to learn the
continuous volumetric function from a limited 2D observations using an
extended MLP. The MLP takes the positional encoded 3D coordinate
$\left(x,y,z\right)$ as input and outputs the image intensity
$I\left(x,y,z\right)$ of the voxel at that location. The positional encoding
is condected via Fourier feature mapping [20] to project the input 3D
coordinates into a higher dimensional space, for enabling the network to learn
sufficient high frequency image features from the LR image stacks. By
minimizing the error between the network prediction and the acquired LR image
intensity across each LR imaging plane, IREM is trained to reconstruct a
continuous model of the observed tissue anatomy. We demonstrate that our
resulting HR MRI reconstruction method quantitatively and qualitatively
outperforms state-of-the-art MISR methods. The major advantage of IREM is
summarized as below:

  * •

No need for large amount of HR data set for training the network;

  * •

No constraints for input LR image resolution and up-sampling scale;

  * •

Comparing with the realistic HR image acquisition with similar total scan
time, the SR image provides equivalent level of image quality and higher SNR.

##  2 Method

Inspired by a very recent implicit scene representation work [13], we build a
deep learning model that learn to represent an image function $I(x,y,z)$ from
input voxel location $(x,y,z)$. In model training, the spatial coordinate of
the LR images are input into a fully-connected network with positional
encoding, while the model is trained to output the image intensity $I(x,y,z)$
of the voxel at that location with minimizing error comparing with the LR
image intensity. Then, the fine-trained network reconstruct the desired HR
image with arbitrary up-sampling SR rate (with arbitrary input $(x,y,z)$). The
spatial locations that are not observed in the LR images are thus the “test
data” generated by the proposed SR model.

Figure 1: An overview of the proposed IREM approach.

An overview of IREM is depicted in Figure 1, which is presented in three
stages: a) Image Spatial Normalization; b) Model Optimization; c) HR Image
Reconstruction.

###  2.1 Image Spatial Normalization

As demonstrated in Figure 1 (a), we first build a set of anisotropic LR
(thick-sliced MR scanning) image stacks $I_{L}=\\{{I}_{C},I_{A},I_{S}\\}$. For
simplifying, we use three orthogonal scanning orientations (coronal, axial,
sagittal) to represent the LR stacks. Then we randomly select one from the
three LR scans as fixed image and rigidly register all the LR scans to the
fixed image to generate transformation matrices
$T_{L2N}=\\{T_{C2N},T_{A2N},T_{S2N}\\}$. Finally, we utilize transformation
matrices $T_{L2N}$ to transfer each LR image from its original space to the
normalized image space $N$.

In the normalized space, the image intensities from different LR image stacks
at same coordinate $(x,y,z)$ represents coherent observations of the same
tissue anatomy in different imaging orientations. Note that the pair of voxel
coordinate and intensity in the normalized space $N$:
$\left(x,y,z\right){\rightarrow}I_{N}\left(x,y,z\right)$ are transformed from
LR scans at position $\left(x_{L},y_{L},z_{L}\right)\rightarrow
I_{L}\left(x_{L},y_{L},z_{L}\right)$ by transformation matrix (i.e.,
$\left(x,y,z\right)=\left(x_{L},y_{L},z_{L}\right)\ast T_{L2N}$ and
$I_{N}\left(x,y,z\right)=I_{L}\left(x_{L},y_{L},z_{L}\right))$.

###  2.2 Model Optimization

The process of model optimization is illustrated in Figure 1 (b). After the
image spatial normalization, the training data set for building IREM model is
prepared. We feed the voxel coordinate $\left(x,y,z\right)$ into IREM to
compute the predicted voxel intensity $I_{N}^{\ast}\left(x,y,z\right)$. Then,
we optimize IREM by minimizing the mean square error loss function between the
predicted voxel intensity and the real observed voxel intensity across the LR
imaging planes. The loss function $L$ is denoted as:

$L(\theta)=\frac{1}{\mathcal{K}}\sum_{(x,y,z)\in
N}\left\|I_{N}(x,y,z)-I_{N}^{*}(x,y,z)\right\|^{2}$ (1)

where $\mathcal{K}$ is mini-batch size, $N$ is normalized image space, and
$\theta$ denotes IREM’s parameters.

####  2.2.1 Architecture of IREM.

As illustrated in Figure 2, IREM consists of a position encoding section (via
Fourier feature mapping [20]), and a fully-connected network (MLP). It takes
3D spatial coordinate as input and outputs the intensity of the voxel at that
location.

Figure 2: Architecture of IREM.

##### Fourier feature mapping.

Theoretically, a simple multi-layer perceptron can approximate any complicated
functions [7]. Recently, Rahaman et al. [16] found that deep learning networks
are biased towards learning lower frequency functions in practical training.
To this end, Mildenhall et al. proposed Fourier feature mapping [20] that maps
the low-dimensional input to a higher dimensional space and thus enable
network to learn higher frequency image feature. In IREM, we perfrom Fourier
feature mapping [20] to map the 3D voxel coordinates to a higher dimensional
space $\mathbb{R}^{2L}(2L\ >3)$ before passing them to the fully-connected
network. Let $\gamma\left(\cdot\right)$ denotes Fourier feature mapping from
the space $\mathbb{R}^{3}$ to $\mathbb{R}^{2L}$ and it is calculated by

$\gamma\left(\mathcal{P}\right)=\left[\cos{\left(2\pi\mathcal{BP}\right)},\sin{\left(2\pi\mathcal{BP}\right)}\right]^{T}$
(2)

where $\mathcal{P}=\left(x,y,z\right)\in\mathbb{R}^{3}$ and each element in
$\mathcal{B}\in\mathbb{R}^{L\times 3}$ is sampled from gaussian distribution
$\mathcal{N}\left(0,1\right)$.

##### Fully-connected network.

The network has eighteen fully-connected layers. Each fully-connected layer is
followed by a batch normalization layer [8] and a ReLU activation [6]. In
order to ease the difficulty of optimizing IREM, we include two skip
connections that concatenate the input of the fully-connected network to the
6th layer’s activation and the 12th layer’s activation, respectively. And
except that there are $2L$, $2L+256$, $2L+256$ neurons in the 1st, 7th, and
13th fully-connected layers. Other fully-connected layers have all 256
neurons.

###  2.3 HR Image Reconstruction

The sketch map of an HR image reconstruction by a fine-trained IREM is shown
in Figure 1 (c). We build a dense grid (i.e., a coordinate set with more
voxels than in LR stacks) in the normalized image space $N$, with isotropic
image resolution. Then, we pass each voxel coordinate $\left(x,y,z\right)$ in
the grid into IREM to produce predicted the voxel intensity
$I_{N}^{\ast}\left(x,y,z\right)$. An HR image thus is reconstructed.

##  3 Experiments

###  3.1 Data

We conduct experiments on three datasets. The details are shown in Table 1.
Dataset #A consists of T1-weighted (T1w) HR brain MR scans from five healthy
adults on a 7T MR scanner and #B consists of 3T T2 flair images from two
patients with lesions in brain white matter. Dataset #C consists of four T1w
brain MR scans (An HR scan used for reference and three orthogonal LR scans
used for input images) from the same volunteer.

Evaluation for the accuracy of implicit HR image representation on simulation
data. For dataset #A and #B, we employ the original MR image as ground truth
(GT) HR image and downsample the image by the factor of $\\{4,\ 8\\}$ in the
three image dimensions, to simulate three orthogonal LR scans. We then use the
built LR stacks to train our model, and compare the reconstructed HR image
with the GT to evaluate the ability of IREM to learn the high frequency image
content in the implicit HR image.

Evaluation for the performance in real data collection protocol. We conduct a
real data collection scene using IREM on dataset #C. In this case, no actual
GT can be built. We scan a HR T1w HR brain image with $0.7$×$0.7$×$0.7$
$mm^{3}$ isotropic in 30.9 mins as HR reference. Then three anisotropic LR
images are scanned in coronal, axial, and sagittal orientations in about 10
mins each, so the total scan time of 3 LR stacks is about 2 mins shorter than
the HR scan.

Table 1: Details of the MR image datasets used in the experiments.

Dataset Name | # Image | Number of Slices | Matrix Size | Voxel Size ($mm^{3}$) | Scan Time
---|---|---|---|---|---
#A | sub1-sub5 | 244 | 244×244 | 0.8×0.8×0.8 | -
#B | sub1-sub2 | 64 | 256×256 | 1×1×2 | -
#C | HR reference | 254 | 368×345 | 0.7×0.7×0.7 | 30.9 mins
axial scan | 80 | 320×320 | 0.7×0.7×2.8 | 10.3 mins
coronal scan | 80 | 320×320 | 0.7×2.8×0.7 | 10.3 mins
sagittal scan | 64 | 320×320 | 2.8×0.7×0.7 | 8.2 mins

###  3.2 Implementation details.

We adopt Adam optimizer [10] to train IREM through back-propagation with a
mini-batch size of 2500, and the hyperparameters of the Adam are set as
follows: $\beta_{1}=0.9,\ \beta_{2}=0.999,\ \varepsilon={10}^{-8}$. The
learning rate starts from ${10}^{-4}$ and decays by factor 0.5 every 500
epochs. For fair comparison, we implement two MISR methods, including Super-
Resolution Reconstruction [4] (SRR) and B-spline interpolation, to compare
with proposed work. We assess quantitatively the performance of the three
methods in terms of peak signal-to-noise ratio (PSNR), structural similarity
(SSIM) [23].

###  3.3 Results

Figure 3: Evolution of the performance of IREM on the different dimensions
$2L$ of Fourier feature mapping [20] for dataset #A. Here the factor $k$ of
down-sampling is 8.

####  3.3.1 Effect of the Fourier feature mapping.

Based on dataset #A, we investigate the effect of Fourier feature mapping [20]
on the performance of IREM. As shown in Figure 3, the performance of IREM
improves with the dimension $2L$ of Fourier feature mapping [20], but its
growth rate gradually decreases, which is consistent with the previous work
[20]. To achieve the balance between efficiency and accuracy, the dimension
$2L$ is set as 256 in all the experiments below.

(a) On dataset #A.

(b) On dataset #B.

Figure 4: Quantitative results (PSNR(dB) / SSMI) of B-spline interpolation,
SRR [4], and IREM on dataset #A and #B. Here $k$ denotes the factor of down-
sampling. Figure 5: Qualitative results of B-spline interpolation, SRR [4],
and IREM on dataset #A (sub1) and #B (sub1). Here the factor $k$ of down-
sampling is 4.

####  3.3.2 Performance on simulation validation.

On dataset #A, the quantitative evaluation results of B-spline interpolation,
SRR [4], and IREM are shown in Figure 4 (a). We indicate that IREM
consistently outperforms the two baselines at all factors $k$ of down-sampling
in terms of all metrics. Figure 5 (a) demonstrates that the representative
slices of an HR image (#A sub1) reconstructed from the three methods, and the
results from IREM visually are closest to GT HR image.

On dataset #B, the quantitative and qualitative evaluation results of the
three methods are shown in Figure 4 (b) and Figure 5 (b), respectively. We can
see that IREM achieves the best performance, which is consistent with the
experiments on dataset #A.

Figure 6: Results of B-spline interpolation, SRR [4], and IREM on dataset #C.

####  3.3.3 Performance in real data collection protocol.

Figure 6 shows the results of the three methods on dataset #C. All MISR
algorithms are conducted after a unique image spatial normalization step.
Comparing with IREM, the image built from B-Spline is more blurry and SRR [4]
yields an image with artifacts. As indicated in the enlarged part, the image
built from IREM achieves equivalent qualitative image details comparing with
the HR reference in the cerebellum, which is one of the most complicated
anatomy in human brain. Besides, benefitting from the multiple anisotropic
thick slice scanning strategy, the SNR in each voxel of the LR image stack is
about 16 times (definition and computational detail of SNR can be found in
[18]) higher than that in the reference HR image. Thus the image contrast
between white matter and gray matter in the reconstructed image is better than
that in the reference HR image. The result suggests that IREM is a more
effective and reasonable pipeline to achieve high-quality HR image comparing
with scanning directly an isotropic HR image.

##  4 Conclusion

In this paper, we proposed IREM, a novel implicit representation based deep
learning framework to improve the performance of multiple image super-
resolution task. In IREM, we combined spatial encoding and fully connected
neural network, then trained a powerful model that precisely predict MR image
intensity from input spatial coordinates. The HR reconstruction result on two
set of simulated image data indicated the ability of IREM to accurately
approach the implicit HR space. While the real scene data collection
demonstrated that based on multiple low-resolution MR scans with thicker
slice, IREM reconstructed HR image with both improved anatomy detail and image
contrast.

##  5 Acknowledgements

This study is supported by the National Natural Science Foundation of China
(No. 62071299, 61901256).

## References

  * [1] Chen, Y., Christodoulou, A.G., Zhou, Z., Shi, F., Xie, Y., Li, D.: Mri super-resolution with gan and 3d multi-level densenet: smaller, faster, and better. arXiv preprint arXiv:2003.01217 (2020)
  * [2] Chen, Y., Shi, F., Christodoulou, A., Zhou, Z., Xie, Y., Li, D.: Efficient and Accurate MRI Super-Resolution Using a Generative Adversarial Network and 3D Multi-level Densely Connected Network, pp. 91–99 (09 2018). https://doi.org/10.1007/978-3-030-00928-1_11
  * [3] Delannoy, Q., Pham, C.H., Cazorla, C., Tor-Díez, C., Dollé, G., Meunier, H., Bednarek, N., Fablet, R., Passat, N., Rousseau, F.: Segsrgan: Super-resolution and segmentation using generative adversarial networks—application to neonatal brain mri. Computers in Biology and Medicine 120, 103755 (2020)
  * [4] Ebner, M., Wang, G., Li, W., Aertsen, M., Patel, P., Aughwane, R., Melbourne, A., Doel, T., Dymarkowski, S., David, A., Deprest, J., Ourselin, S., Vercauteren, T.: An automated framework for localization, segmentation and super-resolution reconstruction of fetal brain mri. NeuroImage 206, 116324 (11 2019). https://doi.org/10.1016/j.neuroimage.2019.116324
  * [5] Gholipour, A., Estroff, J., Warfield, S.: Robust super-resolution volume reconstruction from slice acquisitions: Application to fetal brain mri. IEEE transactions on medical imaging 29, 1739–58 (10 2010). https://doi.org/10.1109/TMI.2010.2051680
  * [6] He, K., Zhang, X., Ren, S., Sun, J.: Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. IEEE International Conference on Computer Vision (ICCV 2015) 1502 (02 2015). https://doi.org/10.1109/ICCV.2015.123
  * [7] Hornik, K., Stinchcomb, M., White, H.: Multilayer feedforward networks are universal approximator. IEEE Transactions on Neural Networks 2 (01 1989)
  * [8] Ioffe, S., Szegedy, C.: Batch normalization: Accelerating deep network training by reducing internal covariate shift (02 2015)
  * [9] Jia, Y., He, Z., Gholipour, A., Warfield, S.: Single anisotropic 3-d mr image upsampling via overcomplete dictionary trained from in-plane high resolution slices. IEEE journal of biomedical and health informatics 20 (08 2015). https://doi.org/10.1109/JBHI.2015.2470682
  * [10] Kingma, D., Ba, J.: Adam: A method for stochastic optimization. International Conference on Learning Representations (12 2014)
  * [11] Lyu, Q., Shan, H., Steber, C., Helis, C., Whitlow, C., Chan, M., Wang, G.: Multi-contrast super-resolution mri through a progressive network. IEEE transactions on medical imaging 39(9), 2738–2749 (2020)
  * [12] Lyu, Q., Shan, H., Wang, G.: Mri super-resolution with ensemble learning and complementary priors. IEEE Transactions on Computational Imaging 6, 615–624 (2020)
  * [13] Mildenhall, B., Srinivasan, P., Tancik, M., Barron, J., Ramamoorthi, R., Ng, R.: NeRF: Representing Scenes as Neural Radiance Fields for View Synthesis, pp. 405–421 (11 2020). https://doi.org/10.1007/978-3-030-58452-8_24
  * [14] Peled, S., Yeshurun, Y.: Superresolution in mri: application to human white matter fiber tract visualization by diffusion tensor imaging. Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine 45, 29–35 (02 2001). https://doi.org/10.1002/1522-2594(200101)45:18216;29::aid-mrm10058217;3.0.co;2-z
  * [15] Pham, C.H., Ducournau, A., Fablet, R., Rousseau, F.: Brain mri super-resolution using deep 3d convolutional networks. In: 2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017). pp. 197–200. IEEE (2017)
  * [16] Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F., Bengio, Y., Courville, A.: On the spectral bias of neural networks. In: International Conference on Machine Learning. pp. 5301–5310. PMLR (2019)
  * [17] Rueda, A., Malpica, N., Romero, E.: Single-image super-resolution of brain mr images using overcomplete dictionaries. Medical Image Analysis 17, 113–132 (10 2012). https://doi.org/10.1016/j.media.2012.09.003
  * [18] Scherrer, B., Gholipour, A., Warfield, S.K.: Super-resolution reconstruction to increase the spatial resolution of diffusion weighted images from orthogonal anisotropic acquisitions. Medical image analysis 16(7), 1465–1476 (2012)
  * [19] Shi, F., Cheng, J., Wang, L., Yap, P.T., Shen, D.: Lrtv: Mr image super-resolution with low-rank and total variation regularizations. IEEE Transactions on Medical Imaging 34 (05 2015). https://doi.org/10.1109/TMI.2015.2437894
  * [20] Tancik, M., Srinivasan, P.P., Mildenhall, B., Fridovich-Keil, S., Raghavan, N., Singhal, U., Ramamoorthi, R., Barron, J.T., Ng, R.: Fourier features let networks learn high frequency functions in low dimensional domains. arXiv preprint arXiv:2006.10739 (2020)
  * [21] Ur Rahman, S., Wesarg, S.: Combining short-axis and long-axis cardiac mr images by applying a super-resolution reconstruction algorithm. Fraunhofer IGD 7623 (03 2010). https://doi.org/10.1117/12.844356
  * [22] Van Reeth, E., Tan, C.H., Tham, I., Poh, C.L.: Isotropic reconstruction of a 4-d mri thoracic sequence using super-resolution. Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine 73 (02 2015). https://doi.org/10.1002/mrm.25157
  * [23] Wang, Z., Bovik, A., Sheikh, H., Member, S., Simoncelli, E.: Image quality assessment: From error measurement to structural similarity. IEEE Trans. Imgage Process. 13 (11 2003)
  * [24] Yang, J., Wright, J., Yu, L.: Image super-resolution via sparse representation. Image Processing, IEEE Transactions on 19, 2861 – 2873 (12 2010). https://doi.org/10.1109/TIP.2010.2050625


11institutetext:  Dr. Karl Remeis-Observatory & ECAP, Friedrich-Alexander
University Erlangen-Nürnberg, Sternwartstr. 7, 96049 Bamberg, Germany
11email<EMAIL_ADDRESS>22institutetext: Institut für Physik und
Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476
Potsdam-Golm, Germany 33institutetext: Department of Physics, University of
Warwick, Gibbet Hill, Road, Coventry, CV4 7AL, UK 44institutetext: Department
of Astrophysics/IMAPP, Radboud University, P. O. Box 9010, NL-6500 GL
Nijmegen, The Netherlands 55institutetext: University, and National
Observatory, of the Chinese Academy of Sciences, Beijing 100012, China
66institutetext: Anton Pannekoek Institute for Astronomy & GRAPPA, University
of Amsterdam, Postbus 94249, 1090 GE Amsterdam, The Netherlands
77institutetext: Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-
Straße 1, 85741 Garching, Germany

# Discovery of a highly magnetic He-sdO star from a double-degenerate binary
merger

M. Dorsch 11 N. Reindl 22 I. Pelisoli 33 U. Heber 11 S. Geier 22 A. G. Istrate
44 S. Justham 556677

(Received ; accepted )

Helium-rich hot subdwarf stars of spectral type O (He-sdO) are considered
prime candidates for stellar merger remnants. Such events should lead to the
generation of strong magnetic fields. However, no magnetic He-sdO has yet been
unambiguously discovered despite the high magnetic rate (20 %) among white
dwarf stars, the progeny of hot subdwarfs. Here we present the discovery of a
strong magnetic field ($B=353\pm 10$ kG) from Zeeman-split hydrogen, helium,
and metal lines in the optical X-SHOOTER spectrum of an He-sdO and present the
first spectroscopic analysis of any magnetic hot subdwarf. For this we used
line-blanketed Tlusty non-local thermodynamic equilibrium models and assumed a
simple homogeneous magnetic field. The derived atmospheric parameters
$T_{\mathrm{eff}}$ = $44900\pm 1000$ K and $\log g$ = $5.93\pm 0.15$ are
typical for He-sdO stars, while the star is less hydrogen-poor than most He-
sdOs at $\log n(\mathrm{He})/n(\mathrm{H})$ = $+0.28\pm 0.10$. The star is a
slow rotator ($v_{\mathrm{rot}}\sin i$ $<$ 40 km s-1). Its chemical
composition is N-rich and C- and O-poor, and the Si and S abundances are close
to solar. Combining the atmospheric parameters with Gaia parallax and
photometry, the stellar radius and luminosity are found to be typical for He-
sdOs and place the star on the helium main sequence in the Hertzsprung-Russell
diagram. Its mass of $0.93^{+0.44}_{-0.30}$ $M_{\odot}$, although uncertain,
appears to be remarkably high. The strong magnetic field along with the
atmospheric parameters and metal abundances provide overwhelming evidence for
the double-degenerate merger scenario.

###### Key Words.:

stars: individual (Gaia DR2 5694207034772278400) — subdwarfs — stars: magnetic
field

##  1 Introduction

Most hot subdwarf-O (sdO) and subdwarf-B (sdB) stars are helium-burning stars
with hydrogen envelopes that are too thin to sustain hydrogen shell burning
(Heber, 2009, 2016). The bulk of sdB stars form the hot extreme of the
horizontal branch at masses close to half solar and radii of 0.1 – 0.2
$R_{\odot}$. Various types of binary interaction can explain the majority of
hot subdwarf stars (Pelisoli et al., 2020). This includes stable Roche
overflow to a main sequence companion star, common-envelope ejection by white
dwarf (WD) companions to red giant stars, and the merging of two helium-core
white dwarfs (He-WDs; Han et al., 2002, 2003; Webbink, 1984). Because hot
subdwarf stars have very thin envelopes, they evolve directly to the WD
sequence without an excursion to the asymptotic giant branch. A smoking gun
for the merger scenario would be the detection of strong magnetic fields,
which are predicted by detailed merger models for massive stars (Schneider et
al., 2019) and carbon-oxygen (CO) WD mergers (Ji et al., 2013; Zhu et al.,
2015) and are likely to occur for helium WD mergers as well.

About 20 % of all WD stars within 20 pc of Earth are known to host surface
magnetic fields, with strengths ranging from a few kilogauss to several
hundred megagauss (e.g. Bagnulo & Landstreet, 2020, 2021). In contrast,
magnetic fields have not yet been directly observed for any hot subdwarf star.
Spectropolarimetry has failed to provide evidence for the magnetic fields of
more than 40 sdB and sdO stars, with upper limits of about 1 to 2 kG
(Landstreet et al., 2012; Mathys et al., 2012; Randall et al., 2015; Bagnulo
et al., 2015). Indirect signs of magnetism have been claimed from light and
radial velocity variations tentatively attributed to magnetic spots (Jeffery
et al., 2013; Geier et al., 2015; Balona et al., 2019; Momany et al., 2020;
Vos et al., 2021), but lack independent confirmation.

A first helium-rich hot O-type subdwarf (He-sdO) star with a significant
magnetic field ($\approx$300 – 700 kG) was identified spectroscopically by
Heber et al. (2013), although this star has not been analysed to date. Here we
present the first spectroscopic analysis of a hot subdwarf star with a strong
magnetic field, the He-sdO Gaia DR2 5694207034772278400 (henceforth
J0809-2627).

##  2 Observations

J0809-2627 was identified as a candidate hot subdwarf by Geier et al. (2019).
We obtained follow-up spectroscopy with the Intermediate Dispersion
Spectrograph (IDS) spectrograph at the Isaac Newton Telescope (INT) to confirm
this classification and identified the star as a helium-rich sdO star. At a
resolution of $\Delta\lambda$ $\approx$ 4 Å and signal-to-noise ratio (S/N) of
50, the IDS spectrum shows strongly broadened hydrogen and helium lines, which
we tentatively interpreted as rapid rotation.

Because hot subdwarf stars are known to be slow rotators (Geier & Heber,
2012), the strong line-broadening of J0809-2627 is unusual. Therefore, the
star was observed again in April 2021 with the medium-resolution spectrograph
X-SHOOTER at the ESO-VLT. We used the 3070 – 7400 Å range of the UVB and VIS
channels, which covers all relevant spectral lines at a mean S/N of about 80
and at a resolving power of $R$ = $\lambda/\Delta\lambda$ $\approx$ $10000$.
At this resolution, which is eight times better than that of IDS, strong
Zeeman-split multiplets are clearly visible, demonstrating that the line
broadening is not due to rapid rotation but is caused by the presence of a
magnetic field (see Fig. 1). The 3070 – 7400 Å spectral range is shown in Fig.
8. All exposure times can be found in Appendix A.1.

Figure 1: Zeeman-split He ii 4686 Å line in the radial-velocity-corrected
X-SHOOTER spectrum of J0809-2627 (grey). The best-fit model is shown in red.
The dashed lines illustrate the positions of the three components that depend
on the magnetic field strength due to the linear Zeeman effect.

Figure 2:  Strongest He ii and H i lines, as well as the best-fitting He i
lines in the X-SHOOTER spectrum of J0809-2627 (grey). The best-fit model is
shown in red. Labels indicate He ii and H i line positions at $B=0$, as well
as lower and upper LS terms for He i.

The light curve of J0809-2627 has been observed by the TESS satellite but does
not show significant variability at a $\approx$ 0.42 parts per thousand (ppt)
detection limit (see Appendix A.2). However, J0809-2627 is located in a
relatively crowded field, so contamination by nearby stars may hide any
variability.

##  3 Models

We used a grid of atmospheric structures computed with the Tlusty code as the
basis of our spectroscopic analysis. A description of the code is given in
Hubeny & Lanz (2017a, b, c). These models are plane-parallel, homogeneous,
hydrostatic, and include H, He, C, N, O, Si, P, S, Fe, and Ni in non-local
thermodynamic equilibrium.

The observed line splittings can be explained by the linear Zeeman effect
caused by the presence of a strong surface magnetic field. The magnetic field
was assumed to be homogeneous across the visible surface, that is, of uniform
strength and direction. Its axis was allowed to be inclined at an angle $\psi$
with respect to the line of sight. Smaller inclinations lead to stronger
displaced Zeeman components relative to the central component. The magnetic
field in J0809-2627 is not strong enough to have a large effect on the
atmospheric structure beyond the additional metal line opacity in the UV
region (Wickramasinghe & Martin, 1986; Tremblay et al., 2015)111 Although iron
and nickel lines are not observed in the optical spectrum of J0809-2627, they
were still included in our models at abundances of three times their solar
value. This is necessary because the opacity of Zeeman-spilt iron group lines
in the UV region likely has a significant effect on the atmospheric structure.
For the same reason, we used a microturbulent velocity of 6 km s-1 for the
atmospheric structure calculations only.. In addition, the linear Zeeman
effect dominates over the quadratic Zeeman effect at strengths below about 1
MG (Garstang & Kemic, 1974). Therefore, our approach was to consider only the
linear Zeeman effect and only in the spectrum synthesis. We modified the
spectrum synthesis code Synspec (Hubeny & Lanz, 2017a) to include linear
Zeeman multiplets for hydrogen, helium, and detectable metal lines, similar to
the method Kawka & Vennes (2011) used to model the cool magnetic WD NLTT
10480, which has a magnetic field strength ($\approx$500 kG) similar to that
of J0809-2627. Polarised radiative transfer in the lines was not included in
our simple models. More details of our model for the magnetic field are given
in Appendix B.

##  4 Spectral fits

Table 1: Atmospheric parameters derived from spectroscopy. Elemental abundance
results are stated by number fraction. Upper limits are given as best-fit
values, while their uncertainties represent values that can clearly be
excluded.

| Spectral fit
---|---
$T_{\mathrm{eff}}$ (K) | $44900\pm 1000$
$\log g$ | $\phantom{+}5.93\pm 0.15$
$\log n(\mathrm{He})/n(\mathrm{H})$ | $+0.28\pm 0.10$
$B$ (kG) | $353\pm 10$
$\psi$ ($\degr$) | $\phantom{0}64\pm 25$
$v_{\mathrm{rad}}$ (km s-1) | $\phantom{0}33\pm 2$
$v_{\mathrm{rot}}\sin i\leavevmode\nobreak\ $ (km s-1) | $<40$

Element | $\log n_{\mathrm{X}}/\sum_{i}n_{i}$
---|---
H | $-0.46\pm 0.07$
He | $-0.18\pm 0.03$
C | ¡$-3.74^{+0.40}$
N | $-2.98\pm 0.20$
O | ¡$-3.92^{+0.30}$
Si | $-4.26\pm 0.30$
S | $-4.16\pm 0.30$

Figure 3: Zeeman-split N iii-iv, Si iv, and S iv lines in the X-SHOOTER
spectrum (grey) and the best-fit model (red) of J0809-2627. The projected
rotational velocity is set to $v_{\mathrm{rot}}\sin i$ = 20 km s-1. Line
positions at zero magnetic field are labelled.

We performed a global $\chi^{2}$ fit to the X-SHOOTER spectrum of J0809-2627
to determine its atmospheric parameters, including $B$ and $\psi$. Some of the
strongest hydrogen and helium lines that could be used for the fit are shown
in Fig. 2. The $\chi^{2}$ fit prefers a projected rotational velocity of
$v_{\mathrm{rot}}\sin i=34$ km s-1. However, $v_{\mathrm{rot}}\sin i$ is not
well constrained because broadening of the displaced Zeeman components may
also be caused by a more complex magnetic field geometry. Possibly for the
same reason, the central components of several helium lines are sharper in the
observation than predicted by our model. The strength of the central
components of Zeeman triplets may also be somewhat increased by magneto-
optical effects (Martin & Wickramasinghe, 1981), which are not included in our
models. Most observed metal lines appear sharp given the limited resolution of
the X-SHOOTER spectra. The final atmospheric parameters as derived from the
X-SHOOTER spectra are listed in the left part of Table 1. All uncertainties
are estimated because they are dominated by systematic effects rather than by
noise.

As demonstrated in Fig. 3, the X-SHOOTER spectrum of J0809-2627 shows strong N
iii-iv, Si iv, and S iv lines. The determination of abundances was complicated
not only by the magnetic field but also because most transitions originate
from high-lying energy levels that are hard to model, even in non-magnetic sdO
stars. This meant that some predicted lines could not be used for the
abundance determination, for example the C iii triplet at 4070 Å and N iii
4379 Å. Metal abundances were therefore estimated by eye, by comparing the
observed spectra with models with different abundances and keeping the
atmospheric parameters fixed to the best-fit values. The derived abundances
are listed on the right side of Table 1. J0809-2627 is enriched in nitrogen at
close to 20 times the solar abundance (by number). Carbon and oxygen seem to
be at least somewhat subsolar because no lines from these elements are clearly
detected – in particular, the C iii 4159 Å, C iv 5805 Å, and O iii 3760 Å
multiplets. The photosphere therefore likely consists of material processed by
hydrogen fusion in the CNO cycle. Silicon and sulphur are enriched at about
two and six times their solar value, respectively.

The X-SHOOTER spectra give us a precise radial velocity for J0809-2627. To
test whether J0809-2627 is radial velocity variable or not, follow-up
spectroscopic observations are necessary because the IDS spectra are of
insufficient quality.

By combining the spectroscopic surface gravity, the Gaia parallax, and the
spectral energy distribution (SED), the mass, radius, and luminosity of
J0809-2627 can be derived. This analysis is described in Appendix C. While the
stellar radius, $0.184\pm 0.011$ $R_{\odot}$, and luminosity,
$123^{+19}_{-16}$ $L_{\odot}$, are typical for He-sdO stars, the derived mass
for J0809-2627, $0.93^{+0.44}_{-0.30}$ $M_{\odot}$, is rather high. It should
be noted that the uncertainty of the mass is considerable due the limited
accuracy of the spectroscopic surface gravity. These parameters place
J0809-2627 close to the zero-age helium main sequence at $\approx$0.8
$M_{\odot}$ (Paczyński, 1971), which means it consists of a helium core with
next to no hydrogen envelope left.

##  5 Evolutionary status

In general, the presence of a strong magnetic field can be considered the
smoking gun for stellar mergers (Schneider et al., 2019). Given that
J0809-2627 is an apparently single magnetic helium-rich sdO star and likely
more massive than usual for that class of star (Schindewolf et al., 2018;
Dorsch et al., 2019, 2020, 2021), the merging of two He-WD stars can be
considered a natural scenario for its formation. In this scenario an object
with a mass as high as 0.93 $M_{\odot}$ would be unlikely, but the scenario is
predicted to produce significant numbers of He-sdO stars with masses of up to
about 0.8 $M_{\odot}$ (Han et al., 2003), which would be well within the
uncertainty of our mass for J0809-2627. A related formation channel for He-sdO
stars is the merger of an He-WD with a post-sdB WD, which would produce a
hybrid WD with a CO core and a thick helium envelope, a so-called HeCO-WD
(Justham et al. 2011, following pioneering models from Iben 1990). This
channel is a natural consequence of one of the common-envelope formation
channels for sdB stars described by Han et al. (2002) and seems likely to
contribute substantially to the population of single He-sdO stars. Moreover,
the inferred effective temperature and surface gravity of J0809-2627 are close
to those of the densest region of the population predictions in Justham et al.
(2011).

Based on the three-dimensional simulations of CO WD mergers of Yoon et al.
(2007) and Lorén-Aguilar et al. (2009), García-Berro et al. (2012) showed that
strong magnetic fields can be generated in the hot, convective, differentially
rotating corona that surrounds the primary during and shortly after the
merging. These magnetic fields are predicted to be frozen to the outer layers
of the final merger product. Magneto-dipole radiation rapidly spins down the
newly formed magnetic WD if magnetic and rotation axes are not aligned. These
predictions appear to be consistent with the recent results of Bagnulo &
Landstreet (2021), who found that young WDs with fields exceeding 1 MG are
more massive than canonical WDs.

The CNO-processed chemical composition and the relatively slow projected
rotation of J0809-2627 provide additional important constraints on merger
models. Evolutionary calculations for mergers of He-WDs with He-WDs were
performed by Zhang & Jeffery (2012), Hall & Jeffery (2016), and Schwab (2018).
Zhang & Jeffery (2012) predict nitrogen-rich surfaces for slow double He-WD
mergers, similar to what we find for J0809-2627. Such systems retain the
initial composition of the secondary He-WD, which includes a CNO-cycle
signature, as observed for J0809-2627. Because our upper limit for the carbon
abundance of J0809-2627 is only slightly subsolar, the composite merger model
of Zhang & Jeffery (2012) cannot be excluded for remnant masses below about
0.7 $M_{\odot}$. This model assumes that more than half of the secondary mass
is rapidly accreted and forms a hot corona, after which the remaining mass is
accreted slowly. As shown in Fig. 4, the 0.8 $M_{\odot}$ slow merger model of
Zhang & Jeffery (2012) matches the effective temperature and luminosity of
J0809-2627 well at the onset of core helium burning.

For hybrid WD mergers, it seems plausible that surface abundances consistent
with J0809-2627 are also possible as long as the He-WD is disrupted rather
than the hybrid HeCO-WD (in other words, as long as the He-WD is the less
massive component of the merging binary). Of the hybrid merger models in
Justham et al. (2011), the luminosity and surface temperature of J0809-2627
are matched by mergers in which a hybrid WD of only 0.35 $M_{\odot}$ accretes
0.15–0.3 $M_{\odot}$ from a merger with a He-WD. Those models are also
consistent with the surface gravity of J0809-2627. However, only the more
massive of those merger models are within the mass range we infer for
J0809-2627.

The best-fit hydrogen abundance for J0809-2627 is higher than that observed
for most (extreme) He-sdO stars (Stroeer et al., 2007; Schindewolf et al.,
2018). Hydrogen was not included in the merger models of Justham et al. (2011)
or Zhang & Jeffery (2012), while both hydrogen and helium were included in the
models of Schwab (2018). The Schwab (2018) models predict that the surface
hydrogen abundance depends on the rotational velocity after the merging –
models that include rotation result in very hydrogen-poor surfaces, which is
not what we find for J0809-2627. Specifically, Schwab (2018) predict
rotational velocities of about 30 km s-1 once a 0.5 $M_{\odot}$ merger remnant
reaches the core-helium-burning phase, which would be consistent with the
elemental abundance pattern of J0809-2627. However, a more massive 0.7
$M_{\odot}$ merger remnant is predicted to be a fast rotator, with
$v_{\mathrm{rot}}$ $\approx$ 100 km s-1. This would require a low inclination
for J0809-2627 given its low projected rotational velocity ($<$ 40 km s-1).
Alternatively, a misalignment of rotation and magnetic axes might provide
further loss of angular momentum through magneto-dipole radiation (García-
Berro et al., 2012), which was not included in the calculations of Schwab
(2018). In the 0.7 $M_{\odot}$ model of Schwab (2018), the surface is
predicted to be enriched in carbon, which, however, is not observed for
J0809-2627.

Figure 4: Position of J0809-2627 (red) in the Hertzsprung-Russell diagram
compared to several evolutionary tracks. Pink lines show two double He-WD
slow-merger tracks from Zhang & Jeffery (2012) for metallicity $Z=0.02$. The
final masses for both tracks are labelled. Black lines show the zero-age
helium main sequence from Paczyński (1971), where black filled circles
indicate masses of 0.5, 0.7, 1.0, and 1.5 $M_{\odot}$. The dashed-dotted grey
lines indicate radii of 0.1 and 1 $R_{\odot}$.

J0809-2627 will evolve through a helium-shell-burning phase and eventually
become a WD. If the magnetic flux is conserved until J0809-2627 reaches the WD
stage ($B\sim R^{-2}$), its surface field strength will reach about 120 MG at
a radius of 0.01 $R_{\odot}$.

##  6 Summary and conclusions

We performed a detailed spectral analysis of the magnetic He-sdO J0809-2627
and used a simple homogeneous model for the magnetic field to derive a mean
magnetic field strength of $353\pm 10$ kG at an inclination of $\psi=64\pm
25\,\degr$. Apart from its strong magnetic field, J0809-2627 has a super-solar
helium abundance at $\log n(\mathrm{He})/n(\mathrm{H})$= $+0.28\pm 0.10$,
which is low compared to typical non-magnetic He-sdO stars. Its effective
temperature, $T_{\mathrm{eff}}$ = $44900\pm 1000$ K, and surface gravity,
$\log g$ = $5.93\pm 0.15$, are not unusual. The He-sdO is nitrogen-rich and
seems to be carbon- and oxygen-poor, although only upper limits could be
derived for C and O. We combined astrometry, photometry, and spectroscopy to
derive stellar parameters of $R$ = $0.184\pm 0.011$ $R_{\odot}$, $L$ =
$123^{+19}_{-16}$ $L_{\odot}$, and $M$ = $0.93^{+0.44}_{-0.30}$ $M_{\odot}$.
Placing it in the Hertzsprung-Russell diagram, J0809-2627 is located on the
helium main sequence, which implies that it has completely lost any hydrogen
envelope.

The TESS light curve of J0809-2627 does not show significant variability (see
Appendix A.2). Assuming a radius of $0.184$ $R_{\odot}$ and a rotational
velocity of less than 40 km s-1, a magnetic spot would be expected to result
in a modulation at a period of more than 5.6 h. The observed lack of
photometric variability may be caused by a low inclination of the rotational
axis, although a low inclination of the magnetic field axis seems to be
excluded by spectroscopy. However, many magnetic main sequence stars are
oblique rotators, which means that their magnetic axis is not aligned with
their rotational axis (Landstreet & Borra, 1978). Photometric observations
that are more sensitive than the $\approx$0.42 ppt detection limit achieved by
the current TESS data would help to further constrain the presence of magnetic
spots.

In summary, the detection of the strong magnetic field of J0809-2627, its
unusually high mass, its chemical composition, its slow projected rotation,
and the lack of evidence for a close stellar companion provide overwhelming
evidence for J0809-2627 being the result of a double-degenerate merger that
produced a stably helium-burning star. Nonetheless, radial velocity monitoring
still has to be done to exclude a close binary nature222He-sdOs with masses as
high as that of J0809-2627 or more could also form from intermediate mass
stars that had their envelope stripped by a compact companion (Götberg et al.,
2018). However, the stripping process is unlikely to generate a strong
magnetic field, so the magnetic field would need to be primordial or generated
during the red giant evolution of the progenitor. In addition, the stripped
stars of Götberg et al. (2018) are predicted to be much more luminous than
J0809-2627.. As noted by Schwab (2018), magneto-hydrodynamic simulations would
be required for suitable double He-WD mergers or He-WD+HeCO-WD mergers,
similar to those performed for CO-WD mergers by Ji et al. (2013) and Zhu et
al. (2015).

The double He-WD and He-WD+HeCO-WD merger scenarios are also both viable for
the more common non-magnetic He-sdO stars as well. Therefore, the question
remains as to why the vast majority of He-sdO stars are apparently non-
magnetic while J0809-2627 has a strong magnetic field. One possibility would
be that the surface magnetic field observed for J0809-2627 is not stable and
weakens quickly after a merging event.

While our simple homogeneous model for the magnetic field generally results in
good fits for most lines in the X-SHOOTER spectrum of J0809-2627, the magnetic
field geometry might be more complicated in reality. Spectropolarimetric
follow-up observations would provide more information about the magnetic field
geometry. The current X-SHOOTER spectrum excludes projected rotation
velocities of more than about 40 km s-1, which is still consistent with
rotation velocities predicted for some double He-WD merger remnants (Schwab,
2018). Higher-resolution monitoring would further constrain this value.
Ultraviolet observations would allow us to improve on the determination of
abundances of carbon and oxygen, for which only upper limits are available at
present, and would also provide abundance measurements for many more metals,
including the iron group. A complete abundance pattern is crucial for
establishing the composition of the accreted material.

###### Acknowledgements.

We thank the referee, John Landstreet, for his very useful comments. We thank
Ylva Götberg for discussions on stripped massive stars and Max Pritzkuleit for
his support with the HOTFUSS project. UH and MD acknowledge funding by the
Deutsche Forschungsgemeinschaft (DFG) through grants IR190/1-1, HE1356/70-1
and HE1356/71-1. IP acknowledges support from the UK’s Science and Technology
Facilities Council (STFC), grant ST/T000406/1. AGI acknowledges support from
the Netherlands Organisation for Scientific Research (NWO). SJ acknowledges
funding from the Netherlands Organisation for Scientific Research (NWO), as
part of the Vidi research program BinWaves (project number 639.042.728, PI: de
Mink). This research has made use of NASA’s Astrophysics Data System. Based on
observations made with ESO Telescopes at the La Silla Paranal Observatory
under programme ID 105.206H.001 and observations with the Isaac Newton
Telescope (under programme ID ING.NL.19B.005) operated on the island of La
Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del
Roque de los Muchachos of the Instituto de Astrofísica de Canarias. This paper
includes data collected by the TESS mission. Funding for the TESS mission is
provided by the NASA’s Science Mission Directorate. This work has made use of
data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Aller et al. (2020) Aller, A., Lillo-Box, J., Jones, D., Miranda, L. F., & Barceló Forteza, S. 2020, A&A, 635, A128
  * Bagnulo et al. (2015) Bagnulo, S., Fossati, L., Landstreet, J. D., & Izzo, C. 2015, A&A, 583, A115
  * Bagnulo & Landstreet (2020) Bagnulo, S. & Landstreet, J. D. 2020, A&A, 643, A134
  * Bagnulo & Landstreet (2021) Bagnulo, S. & Landstreet, J. D. 2021, MNRAS, 507, 5902
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58
  * Balona et al. (2019) Balona, L. A., Handler, G., Chowdhury, S., et al. 2019, MNRAS, 485, 3457
  * Condon & Shortley (1935) Condon, E. U. & Shortley, G. H. 1935, The Theory of Atomic Spectra
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, 2246
  * DENIS Consortium (2005) DENIS Consortium. 2005, VizieR Online Data Catalog: The DENIS database
  * Dorsch et al. (2021) Dorsch, M., Jeffery, C. S., Irrgang, A., Woolf, V., & Heber, U. 2021, A&A, 653, A120
  * Dorsch et al. (2019) Dorsch, M., Latour, M., & Heber, U. 2019, A&A, 630, A130
  * Dorsch et al. (2020) Dorsch, M., Latour, M., Heber, U., et al. 2020, A&A, 643, A22
  * Drlica-Wagner et al. (2021) Drlica-Wagner, A., Carlin, J. L., Nidever, D. L., et al. 2021, ApJS, 256, 2
  * El-Badry et al. (2021) El-Badry, K., Rix, H.-W., & Heintz, T. M. 2021, MNRAS, 506, 2269
  * Feinstein et al. (2019) Feinstein, A. D., Montet, B. T., Foreman-Mackey, D., et al. 2019, PASP, 131, 094502
  * Fitzpatrick et al. (2019) Fitzpatrick, E. L., Massa, D., Gordon, K. D., Bohlin, R., & Clayton, G. C. 2019, ApJ, 886, 108
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * García-Berro et al. (2012) García-Berro, E., Lorén-Aguilar, P., Aznar-Siguán, G., et al. 2012, ApJ, 749, 25
  * Garstang & Kemic (1974) Garstang, R. H. & Kemic, S. B. 1974, Ap&SS, 31, 103
  * Geier & Heber (2012) Geier, S. & Heber, U. 2012, A&A, 543, A149
  * Geier et al. (2015) Geier, S., Kupfer, T., Heber, U., et al. 2015, A&A, 577, A26
  * Geier et al. (2019) Geier, S., Raddi, R., Gentile Fusillo, N. P., & Marsh, T. R. 2019, A&A, 621, A38
  * Götberg et al. (2018) Götberg, Y., de Mink, S. E., Groh, J. H., et al. 2018, A&A, 615, A78
  * Hall & Jeffery (2016) Hall, P. D. & Jeffery, C. S. 2016, MNRAS, 463, 2756
  * Han et al. (2003) Han, Z., Podsiadlowski, P., Maxted, P. F. L., & Marsh, T. R. 2003, MNRAS, 341, 669
  * Han et al. (2002) Han, Z., Podsiadlowski, P., Maxted, P. F. L., Marsh, T. R., & Ivanova, N. 2002, Monthly Notices of the Royal Astronomical Society, 336, 449
  * Heber (2009) Heber, U. 2009, Annual Review of Astronomy and Astrophysics, 47, 211
  * Heber (2016) Heber, U. 2016, Publications of the Astronomical Society of the Pacific, 128, 082001
  * Heber et al. (2013) Heber, U., Geier, S., & Gaensicke, B. 2013, in European Physical Journal Web of Conferences, Vol. 43, European Physical Journal Web of Conferences, 04002
  * Heber et al. (2018) Heber, U., Irrgang, A., & Schaffenroth, J. 2018, Open Astronomy, 27, 35
  * Hönl (1925) Hönl, H. 1925, Zeitschrift für Physik, 31, 340
  * Hubeny & Lanz (2017a) Hubeny, I. & Lanz, T. 2017a, ArXiv e-prints [arXiv:1706.01859]
  * Hubeny & Lanz (2017b) Hubeny, I. & Lanz, T. 2017b, ArXiv e-prints [arXiv:1706.01935]
  * Hubeny & Lanz (2017c) Hubeny, I. & Lanz, T. 2017c, ArXiv e-prints [arXiv:1706.01937]
  * Iben (1990) Iben, Icko, J. 1990, ApJ, 353, 215
  * Irrgang et al. (2021) Irrgang, A., Geier, S., Heber, U., et al. 2021, A&A, 650, A102
  * Jeffery et al. (2013) Jeffery, C. S., Ramsay, G., Naslim, N., et al. 2013, MNRAS, 429, 3207
  * Ji et al. (2013) Ji, S., Fisher, R. T., García-Berro, E., et al. 2013, ApJ, 773, 136
  * Justham et al. (2011) Justham, S., Podsiadlowski, P., & Han, Z. 2011, MNRAS, 410, 984
  * Kawka & Vennes (2011) Kawka, A. & Vennes, S. 2011, A&A, 532, A7
  * Landstreet et al. (2012) Landstreet, J. D., Bagnulo, S., Fossati, L., Jordan, S., & O’Toole, S. J. 2012, A&A, 541, A100
  * Landstreet & Borra (1978) Landstreet, J. D. & Borra, E. F. 1978, ApJ, 224, L5
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4
  * Lorén-Aguilar et al. (2009) Lorén-Aguilar, P., Isern, J., & García-Berro, E. 2009, A&A, 500, 1193
  * Magnier et al. (2020) Magnier, E. A., Schlafly, E. F., Finkbeiner, D. P., et al. 2020, ApJS, 251, 6
  * Martin & Wickramasinghe (1981) Martin, B. & Wickramasinghe, D. T. 1981, MNRAS, 196, 23
  * Mathys et al. (2012) Mathys, G., Hubrig, S., Mason, E., et al. 2012, Astronomische Nachrichten, 333, 30
  * Moehler et al. (2014) Moehler, S., Modigliani, A., Freudling, W., et al. 2014, A&A, 568, A9
  * Momany et al. (2020) Momany, Y., Zaggia, S., Montalto, M., et al. 2020, Nature Astronomy, 4, 1092
  * Onken et al. (2019) Onken, C. A., Wolf, C., Bessell, M. S., et al. 2019, PASA, 36, e033
  * Paczyński (1971) Paczyński, B. 1971, Acta Astron., 21, 1
  * Pelisoli et al. (2020) Pelisoli, I., Vos, J., Geier, S., Schaffenroth, V., & Baran, A. S. 2020, A&A, 642, A180
  * Randall et al. (2015) Randall, S. K., Bagnulo, S., Ziegerer, E., Geier, S., & Fontaine, G. 2015, A&A, 576, A65
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3
  * Schindewolf et al. (2018) Schindewolf, M., Németh, P., Heber, U., et al. 2018, A&A, 620, A36
  * Schlafly et al. (2018) Schlafly, E. F., Green, G. M., Lang, D., et al. 2018, ApJS, 234, 39
  * Schneider et al. (2019) Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211
  * Schwab (2018) Schwab, J. 2018, MNRAS, 476, 5303
  * Stroeer et al. (2007) Stroeer, A., Heber, U., Lisker, T., et al. 2007, Astronomy and Astrophysics, 462, 269
  * Tremblay et al. (2015) Tremblay, P. E., Fontaine, G., Freytag, B., et al. 2015, ApJ, 812, 19
  * Vos et al. (2021) Vos, J., Pelisoli, I., Budaj, J., et al. 2021, A&A, 655, A43
  * Webbink (1984) Webbink, R. F. 1984, The Astrophysical Journal, 277, 355
  * Wickramasinghe & Martin (1986) Wickramasinghe, D. T. & Martin, B. 1986, MNRAS, 223, 323
  * Yoon et al. (2007) Yoon, S. C., Podsiadlowski, P., & Rosswog, S. 2007, MNRAS, 380, 933
  * Zhang & Jeffery (2012) Zhang, X. & Jeffery, C. S. 2012, Monthly Notices of the Royal Astronomical Society, 419, 452
  * Zhu et al. (2015) Zhu, C., Pakmor, R., van Kerkwijk, M. H., & Chang, P. 2015, ApJ, 806, L1

##  Appendix A Observations

###  A.1 Spectroscopy

Details of the spectroscopic observations are listed in Table 2.

Table 2: Dates and exposure times for the spectroscopic observations of
J0809-2627.

Spectrograph | Mid exposure BJD | Exposure time (s)
---|---|---
INT/IDS | 2458837.6695 | 750
INT/IDS | 2458837.6798 | 750
VLT/X-SHOOTER | 2459314.5691 | 1800

###  A.2 TESS

J0809-2627 was in the TESS field during sectors 7, 8, and 34. The star was not
targeted for short-cadence observations, but the full-frame images are
available. They have an integration time of 30 min for sectors 7 and 8, and 10
min for sector 34. The total time span of the light curve is 760 days. We
performed the photometry of all available images using eleanor (Feinstein et
al. 2019). As the field is relatively crowded (see Fig. 5), we used the small
aperture option, which limits aperture size to 8 pixels333The TESS pixel size
is 21”.. We then carried out the Fourier transform of all the data up to the
Nyquist frequency, as shown in Fig. 6. In order to determine a detection
limit, we randomly assigned one of the measured fluxes to each measured time
and recalculated the Fourier transform a thousand times, each time recording
the maximum amplitude. We selected the maximum overall value as the detection
limit; as it happens once every 1000 runs, this corresponds to a false-alarm
probability of $1/1000$. No peaks are detected above this threshold;
therefore, we find no periodic variability for periods between $\approx 20$
min and $\approx 380$ days down to a detection limit of $\approx 0.42$ ppt. It
is important to emphasise, however, that the large pixel size can result in
contamination from nearby constant stars, which can hide any variability.
Follow-up with better image resolution is encouraged.

Figure 5: Image from TESS sector 7 created with tpfplotter (Aller et al. 2020)
showing the location of J0809-2627 (white cross) and other stars with a
magnitude difference of six or less (red circles). Figure 6: Fourier transform
of the TESS data. The amplitude is shown in ppt. The adopted detection limit
of 0.42 ppt is shown as a dashed red line, and no peaks are detected above
this limit.

##  Appendix B Linear Zeeman effect

We assume LS coupling for all metal line transitions given that the heaviest
observed element is sulphur. An LS coupling state is described by its orbital
angular momentum, $L$, its spin angular momentum, $S$, and its total angular
momentum, $J$. For each included line, the upper and lower LS states were
split into $2J+1$ components with magnetic quantum numbers from $m=-J$ to
$m=+J$. Only transitions with $\Delta m=0$ or $\pm$1 were considered.

For a transition between lower and upper states with magnetic quantum numbers
$m_{\mathrm{l}}$ and $m_{\mathrm{u}}$, the wavelength shift with respect to
the rest wavelength at zero magnetic field, $\lambda_{0}$, is given as

$\Delta\lambda=\frac{eB\lambda_{0}^{2}}{4\pi
m_{e}c^{2}}\cdot(m_{\mathrm{l}}g_{\mathrm{l}}-m_{\mathrm{u}}g_{\mathrm{u}}),$
(1)

where $e$ is the elementary charge, $B$ is the magnetic field strength,
$m_{e}$ is the electron mass, $c$ is the speed of light, and $g_{\mathrm{l}}$
and $g_{\mathrm{u}}$ are the Landé $g$-factors for the lower and upper level.

In the case of LS coupling, the Landé $g$-factor is given as

$g=1+\frac{J(J+1)+S(S+1)-L(L+1)}{2J(J+1)}.$ (2)

Transitions between levels of hydrogen and singly ionised helium represent
special cases in this model for the magnetic field because these levels are
sufficiently described by a main quantum number, $n$. For these levels, as
well as neutral helium, the Landé $g$-factors were set to unity.

Table 3: Relative intensities for the $\pi$ and $\sigma$ components following Hönl (1925). The quantum numbers $J$ and $m$ refer to the initial level, and $\Delta J=J_{\mathrm{final}}-J_{\mathrm{initial}}$. We note the difference from Eq. (5) of Kawka & Vennes (2011), which contains a sign error.  $\Delta J$ | $\pi$ | $\pm\sigma$
---|---|---
$0$ | $m^{2}$ | $\frac{1}{4}(J\mp m)(J\pm m+1)$
$+1$ | $(J+1)^{2}-m^{2}$ | $\frac{1}{4}(J\pm m+1)(J\pm m+2)$
$-1$ | $J^{2}-m^{2}$ | $\frac{1}{4}(J\mp m)(J\mp m-1)$

Analytic expressions for the relative intensities of the components of any
linear Zeeman multiplet were first presented by Hönl (1925). They are
summarised in Table 3. Unlike the Landé $g$-factors, these relative
intensities depend only on $J$ and $m$, and are therefore valid in any
coupling scheme.

The relative intensity ($I$) of the central ($\pi$) and displaced ($\sigma$)
components depends on the angle ($\psi$) between the magnetic field axis and
the line of sight. The angular dependence given by Condon & Shortley (1935) is
for the $\pi$ (central) component

$\frac{I(\psi)}{I(90\degr)}=\sin^{2}\psi$ (3)

and for the $\sigma$ components

$\frac{I(\psi)}{I(90\degr)}=1+\cos^{2}\psi.$ (4)

##  Appendix C Mass, radius, luminosity, and stellar population
characteristics

Figure 7: SED fit for J0809-2627. The grey line shows the final model of J0809-2627, and filter-averaged flux measurements are indicated by dashed horizontal lines. Photometric surveys are identified by colour codes: Pan-STARRS (dark red; Magnier et al. 2020), SkyMapper (dark gold; Onken et al. 2019), Gaia (cyan; Riello et al. 2021), DECam (gold; Schlafly et al. 2018; Drlica-Wagner et al. 2021), DENIS (orange; DENIS Consortium 2005), and 2MASS (red; Cutri et al. 2003).  Table 4: SED fit and stellar parameters. The stellar parameters result from the SED fit with a prescribed $T_{\mathrm{eff}}$ = $44900\pm 1000$ K. The mode and the highest density interval of each quantity are given for 1$\sigma$ probability (see Bailer-Jones et al. 2018). | SED fit
---|---
$T_{\mathrm{eff}}$ (K) | $47000^{+11000}_{-3000}$
$\log\Theta\,\mathrm{(rad)}$ | $-11.246^{+0.005}_{-0.004}$
$E_{44-55}$ (mag) | $0.073^{+0.004}_{-0.004}$
Radius $R=\Theta/(2\varpi)$ |  $0.184^{+0.011}_{-0.010}$ $R_{\odot}$
Mass $M=gR^{2}/G$ |  $0.93^{+0.44}_{-0.30}$ $M_{\odot}$
Luminosity $L/L_{\odot}=(R/R_{\odot})^{2}(T_{\mathrm{eff}}/T_{\mathrm{eff},\odot})^{4}$ | $123^{+19}_{-16}$

We constructed the SED of J0809-2627 by collecting photometric measurements
from several surveys. Two SED fits were then performed using the same grid as
for the spectroscopic analysis, with $T_{\mathrm{eff}}$ either as a free
parameter or fixed to the value derived by the spectroscopic fit (for details,
see Heber et al. 2018; Irrgang et al. 2021). The latter approach is useful
because $T_{\mathrm{eff}}$ $\gtrsim$ 35 000 K are not well constrained from
optical and infrared magnitudes alone. The surface gravity and helium
abundance were always fixed to spectroscopic values, and the angular diameter
$\Theta$ was always a free parameter. Interstellar extinction was accounted
for using the reddening law of Fitzpatrick et al. (2019) with the colour
excess $E_{44-55}$ as a free parameter. An extinction parameter of
$R(55)=3.02$ was assumed.

The angular diameter derived from the SED (with fixed $T_{\mathrm{eff}}$) was
combined with the parallax measurement provided by Gaia Early Data Release 3
(EDR3; Gaia Collaboration et al. 2016, 2021) and the spectroscopic surface
gravity to estimate the stellar parameters mass, radius, and luminosity for
J0809-2627. Here, the parallax measurement, $\varpi=0.68\pm 0.04$ mas, was
corrected for its zero point offset following Lindegren et al. (2021), and the
corresponding uncertainty was inflated using the function of El-Badry et al.
(2021). The SED fit and stellar parameters are summarised in Table 4. Figure 7
shows the SED fit with $T_{\mathrm{eff}}$ fixed to the spectroscopic value.

Because J0809-2627 is close to the Galactic disk, we expect it to belong to
the thin disk stellar population I. This can be confirmed by calculating its
current Galactic space velocity vector from the radial velocity (assumed to be
non-variable), parallax, and proper motions provided by Gaia EDR3 (Gaia
Collaboration et al. 2021). The resulting vector ($U$, $V$, $W$) = ($-17.3$
$\pm$ 2.0, 230.1 $\pm$ 2.7, 7.1 $\pm$ 0.8) km s-1, where $U$ is measured
towards the Galactic centre, $V$ in the direction of Galactic rotation, and
$W$ perpendicular to the disc, is consistent with a membership of the thin
disk population.

##  Appendix D Spectral comparisons

Here we show spectral comparisons for the X-SHOOTER spectra of J0809-2627,
specifically the UVB and VIS spectral arms up to 7400 Å. The strongest lines
are labelled at their position for zero magnetic field. Spectral regions that
are not well reproduced by our models were excluded from the fit. This
includes several metal lines, as well as some neutral helium lines. Especially
narrow He i lines from the triplet ($S$=1) system, such as He i 5876 Å, and
lines with strong forbidden components, such as He i 4472, 4922 Å, are not
well modelled. This is to be expected given the simple nature of our model for
the magnetic field, as well as the lack of helium broadening tables in the
presence of a magnetic field. The C iii lines at 4070 Å are poorly modelled
due to limitations of the Tlusty model ion. Telluric lines were removed using
the grid of transmission spectra provided by Moehler et al. (2014).

Figure 8: X-SHOOTER spectrum of J0809-2627 (grey) and the best-fit model
(red). Figure 9: X-SHOOTER spectrum of J0809-2627 (grey) and the best-fit
model (red). Figure 10: X-SHOOTER spectrum of J0809-2627 (grey) and the best-
fit model (red).


# dark-modes excitation in coupled nano Fabry-Perot structure induced by
symmetry breaking

Baptiste Fix<EMAIL_ADDRESS>Julien Jaeck  Patrick Bouchon  DOTA,
ONERA, Université Paris Saclay, F-91123 Palaiseau - France  Nathalie Bardou
Center for Nanoscience and Nanotechnology (C2N) - CNRS, Univ. Paris-Sud,
Université Paris-Saclay, 91460 Marcoussis, France  Benjamin Vest  Laboratoire
Charles Fabry, 2 avenue Augustin Fresnel, CNRS, Université Paris-Saclay, 91127
Palaiseau Cedex, France  Riad Haidar  Département de Physique, École
Polytechnique, 91128 Palaiseau, France  DOTA, ONERA, Université Paris Saclay,
F-91123 Palaiseau - France

###### Abstract

Absorbing metamaterials have attracted a widespread interest but suffer from
the metallic ohmic losses that limits their Q-factor to values around 10. In
this paper, we report the engineering and experimental characterization of a
high Q dark-mode in a metasurface resulting from the coupling of several
individual nano-Fabry-Perot structures supporting low-contrast resonances. The
resulting resonance has a Q-factor above 20 and a nearly total absorption.
Moreover, we also demonstrate that several dark-mode resonances can be excited
simultaneously, by exploiting the coupling of various pairs of nano Fabry-
Perot resonators within the same subwavelength period. Finally, the
resonance’s Q-factor depends monotonously on the thickness of the metasurface
and can reach values up to 85.

Absorbing metamaterials have the ability to turn an otherwise reflective
surface into a nearly perfect absorber for a given wavelength and
polarization. Among those metamaterials, the MIM (metal-insulator-metal)
resonator has attracted a widespread interestCui  _et al._ (2014) due to its
ease of fabrication and the possibility to use large scale fabrication
technique such as colloids Moreau  _et al._ (2012) or nanoimprint Cattoni  _et
al._ (2011). They have been used in various applications, e.g., infrared
detection Palaferri  _et al._ (2018), nonlinear optics, thermal emission
Makhsiyan  _et al._ (2015), biosensing Cattoni  _et al._ (2011). In most
cases, a noble metal is used, which limits the ohmic losses of the structure,
and leads to typical Q-factors ( corresponding to the ratio
$\frac{\lambda}{\Delta\lambda}$) of 10 in the optical range. For applications
taking advantages of light matter interactions or targeting either a narrower
spectral band, these losses can be detrimental Khurgin (2015).

Figure 1: (a) Spectral reflectivity under TM-polarised illumination of two
standard critically coupled MIM resonators ($h=45\,\mathrm{nm}$, $W_{1}=475$
and $W_{2}=525\,\mathrm{nm}$). Their combination inside a single periodic
pattern exhibits two resonances and leads to a widening of the absorbing
window. (b) Spectral reflectivity under TM-polarised illumination of two
under-coupled Fabry-Perot resonators ($h=280\,\mathrm{nm}$, $W_{A}=400$ and
$W_{B}=495\,\mathrm{nm}$). Their combination inside a single periodic pattern
exhibits an unexpected resonance at $\lambda_{r}=4.23\,\mathrm{\mu m}$ with a
minimum of reflectivity of $14\%$ and a Q-factor of 22.

An interesting solution to achieve high Q-factor resonances despite the
metallic losses is to engineer the structure in order to enable the excitation
of a dark-mode, which, due to its non-radiating nature, has inherently a long
lifetime and a high Q-factor. For plasmonic structures, a solution to ensure
the existence of a dark-mode is to couple two similar resonators. Thus, the
two individual modes hybridize into a symmetric radiant bonding mode and a
dark anti-bonding modeJain and El-Sayed (2010). Such coupled structures have
been evidenced for gold nano-particles Deng  _et al._ (2018); Sheikholeslami
_et al._ (2010), nanorods Gao  _et al._ (2018); Osberg  _et al._ (2014), grids
Christ  _et al._ (2008); Zhou and Odom (2011) and metamaterials Omaghali  _et
al._ (2014). Due to their anti-symmetry, dark anti-bonding modes cannot be
excited. Yet, it has been demonstrated that creating a grey-mode could enable
a free space coupling while maintaining, at least partially, the dark-mode
properties. Omaghali  _et al._ (2014); Dayal  _et al._ ; Al-Naib (2020)
Several excitation strategies have been used in the literature to achieve it,
such as spatially inhomogeneous Huang  _et al._ (2010), evanescent waves Yang
_et al._ (2010), high-angle-of-incidence excitation Zhou and Odom (2011).
Other strategies have been proposed such as coupling the dark-mode to a
bright-mode Gómez  _et al._ (2013) or to modify the symmetry of the structure
itself Jeyaram  _et al._ (2013); Omaghali  _et al._ (2014); Cao  _et al._
(2012).

Therefore, finding a way to use such a dark-mode in a MIM geometry is very
appealing to overcome Q-factor limitations experienced by usual Fabry-Perot-
like resonances of MIM structures. However, several studies have been done on
symmetry broken MIM metasurfaces, but none have reported the excitation of
dark-modes. It is even the other way around, as combination of MIM resonators
have been shown to behave nearly independently and used to tailor broadband
absorbers Cui  _et al._ (2011); Liu  _et al._ (2011); Bouchon  _et al._
(2012); Cui  _et al._ (2014). Here we introduce a general framework for dual
MIM to sustain an excitable anti-bonding dark-mode. It is based on the
coupling of two undercoupled MIM Fabry-Perot cavities similarly to what was
observed with coupled metallic groovesFix  _et al._ (2017). This dark-mode is
excitable even at normal incidence and we show that the Q-factor value can be
increased up to 85 choosing the appropriate configuration. Finally, we
demonstrate how several dark-modes can be tailored in the same metasurface.

The MIM resonator geometry, depicted in the insets of Fig 1, is made of a
stack of a continuous metallic layer, thick enough to be optically opaque, a
dielectric layer (thickness $h$) and a patterned metallic layer (thickness
$t$). Classically, the geometrical parameters are set so that near critical
coupling is achieved, i.e., $R\simeq 0$. For instance, Fig 1.a shows the
computed spectra of two critically coupled resonators centered at
$4.2\,\mathrm{\mu m}$ and $4.6\,\mathrm{\mu m}$ ($W_{1}=475\,\mathrm{nm}$,
$W_{2}=525\,\mathrm{nm}$, $h=45\,\mathrm{nm}$, $t=50\,\mathrm{nm}$, period
$P_{1}=885\,\mathrm{nm}$). At resonance, the full width at half maximum (FWHM)
of the first resonator is $\Delta\lambda=400\,\mathrm{nm}$, which corresponds
to a quality factor below 10. These two resonators can be combined in the same
subwavelength period ($P_{2}=2P_{1}=1770\,\mathrm{nm}$), leading to two
resonances with a broadening of the absorption band
($\Delta\lambda=710\,\mathrm{nm}$). The combined resonators spectrum is well
described by the product of the individual resonators spectra Koechlin  _et
al._ (2013).

In contrast, solid lines on Fig 1.b show the numerically simulated spectra of
two undercoupled nano Fabry-Perot resonators. A resonator is undercoupled when
the radiative losses are too high with regards to the internal losses Seok
_et al._ (2011). This is achieved by dramatically increasing the thickness of
the SiC layer up to $h=280\,\mathrm{nm}$, which leads to reflectivity higher
than 80% for both resonator A ($W_{A}=400\,\mathrm{nm}$ ,
$P_{1}=885\,\mathrm{nm}$) and resonator B ($W_{B}=495\,\mathrm{nm}$,
$P_{1}=885\,\mathrm{nm}$). Surprisingly, when both resonators are combined in
the same subwavelength period ($P_{2}=1770\,\mathrm{nm}$) to form a new
resonator geometry called AB bi-FP, the resulting spectrum is no longer a
simple product of the two individual resonators spectra. The reflectivity
spectrum shows a blue-shifted resonance near critical coupling ($R=14\%$),
with a narrower FWHM ($\Delta\lambda=190\,\mathrm{nm}$), thus a higher
Q-factor of 22. The optical characterizations of the resonators far-field
reflectivity (see dashed lines on Fig 1.b) confirm the numerical prediction.
Indeed, the agreement between the numerical prediction and the reflectivity
spectrum of the AB bi-FP structure for a TM polarised light source with an
incidence of $13^{\circ}$ is very good. The slight mismatch between
$4.5\,\mathrm{\mu m}$ to $5.5\,\mathrm{\mu m}$, is attributed to absorption of
the SiC layer in this spectral range.

Figure 2:  Computed reflectivity spectrum under TM-polarised illumination at
normal incidence and with regards to the incidence angle of (a-b) the A
single-FP, (d-e) the AA bi-FP and (h-i) the AB bi-FP stuctures. In dashed
lines, the computed value of the complex electromagnetic mode (reflectivity
poles) of the structures. (c) Illustration of the bonding and anti-bonding
mode for two coupled nano-FP cavities. (f-g), (j-k) Cartography of the $E_{z}$
field for the two modes of the bi-FP structures.

In the following of the paper, we discuss the origin, properties and
engineering of this narrow resonance in a variety of situations combining
individual FP resonators. single-FP, bi-FP, and tri-FP stands for periodic
structure with one, two and three ribbons in a given period. The widths of the
ribbons are $W_{A}$, $W_{B}$ and $W_{C}$ and for the sake of clarity, the
single-FP are named A, B and C, the bi-FP are named AB and BC and the tri-FP
ABC. In the specific case where two ribbons are identical, but not equally
spaced, the structure is named AA bi-FP structure. The single-FP structures
have a period $P_{1}=885\,\mathrm{nm}$ while the period of the bi-FP
structures is $P_{2}=2P_{1}$ with the ribbons evenly distributed in the
period.

We now discuss the origin of the narrow resonance of the AB structure on Fig
1.b. The complex electromagnetic modes of the A single-FP structure, and the
AA and AB bi-FP structures have been calculated as a function of the
illumination angles (see dashed lines in Fig 2.b,e,i).

These modes have been super-imposed with the structures reflectivity spectra
at these illumination angles. In the case of the single-FP structure, only the
usual bright gap-plasmon (green dashed line) mode exist and correspond to a
highly under-coupled resonance ($R_{min}\simeq$ 90%), as shown in Fig 2a).

For the bi-FP structures, two modes exist and correspond respectively to the
bonding and anti-bonding mode. They result from the hybridization of the
resonances of individual FP resonators combined in a given period (see Fig
2.c). A numerical search of the reflectivity poles of the metasurface allows
to identify these two modes. In the case of AA bi-FP structure, the bonding,
lower-energy mode (green dashed line on Fig 2.e), provides a low-Q resonance
centered around 4.5 um at normal incidence (see Fig 2.d), corresponding to the
single-FP resonance, blueshifting at higher incidences. The mode exhibits a
symmetric field with respect to the middle plane of a period (see Fig 2.f).

In contrast, the anti-bonding, higher-energy only exhibits a narrow
reflectivity drop at non-zero incidence (orange dashed line on Fig 2.d). This
is a direct consequence of the antisymmetric field distribution of this mode,
which make it dark under a symmetric, normal incidence excitation (see Fig
2.g).

Interestingly, symmetry breaking induced by oblique excitation allows to
couple far field radiation into the mode volume. The resulting resonance
exhibits a high Q-factor of about 50 to be compared with the bonding mode
resonance with a Q-factor around 10.

In the case of the AB bi-FP structure experimentally studied in the Fig 1.b,
the modification of one ribbon width break the middle plane of symmetry
(Reflectivity spectra, map and mode structures are depicted on Fig 2.h-k).
This enables the excitation of the "dark-mode" (now strictly a grey-mode) at
any incidence angle. The increased coupling between the anti-bonding mode and
the free space decreased the resonance Q-factor compared to the case of the
"purely dark-mode" of the AA bi-structure. Still, the Q-factor of the
resonance is about 25 which is still 2 to 3 times higher than the gap-plasmon
resonances.

Figure 3: (a) Measured (dashed) and computed (solid) reflection spectra of the
AB (red) and BC (orange) bi-FP structures at $13^{\circ}$ of incidence. The
resonance wavelength of the bi-FP structure is tunable by the modification of
one or two ribbon widths.
($\mathrm{W_{A}}$=400$\,\mathrm{nm}$,$\mathrm{W_{B}}$=495$\,\mathrm{nm}$,$\mathrm{W_{C}}$=600$\,\mathrm{nm}$)
(b) Measured reflectivity map of the AB and BC bi-FP structure for a TM
polarised illumination between $13^{\circ}$ and $70^{\circ}$ of incidence. The
dotted lines represent the evolution of the resonance wavelength obtained by
the numerical simulation.

The anti-bonding resonance has three appealing properties. First, the
resonance wavelength can be tuned by adjusting the dimensions of the pattern
while maintaining the critical coupling. Secondly, it is possible to introduce
more anti-bonding resonances by adding under-coupled FP cavities inside the
periodic pattern. Thirdly, the quality factor can be tuned independently from
the resonance wavelength.

The Fig 3.a demonstrates the tunability of the resonance: it compares the
numerically simulated and experimentally measured reflectivity spectrum at
$13^{\circ}$ of incidence of the AB and BC bi-FP structures, the two
structures only differing by the width of one of their two ribbons. Both
spectra have a similar shape, with a critically coupled resonance with $14\%$
(resp $21\%$) reflectivity and a Q-factor of 22 (resp 25) but at two different
resonant wavelengths $\lambda_{r}=4.23\,\mathrm{\mu m}$ and
$\lambda_{r}=4.61\,\mathrm{\mu m}$. The measured reflectivity diagram as a
function of the angle of incidence is depicted in Fig 3.b for TM-polarised
illuminations for the AB and BC bi-FP resonators. The dotted lines represent
the computed resonant wavelengths. Both resonators have a similar angular
behavior, mostly insensitive to the angle of incidence below 20 degrees, and
slightly red-shifting above.

Figure 4: (a) Scheme of the ABC coupled-FP structure with period
$P_{3}=3P_{1}$, $\mathrm{W_{A}}=405\,\mathrm{nm}$,
$\mathrm{W_{B}}=500\,\mathrm{nm}$, $\mathrm{W_{C}}=600\,\mathrm{nm}$. (b)
Measured and computed reflection spectra of the ABC coupled-FP structure at
$13^{\circ}$ incidence. Two resonances are located at
$\lambda_{1}=4.43\,\mathrm{\mu m}$ and $\lambda_{2}=5.12\,\mathrm{\mu m}$ with
minima of reflectivity of $15\%$ and Q-factors of 15 and 11 respectively. The
inset is an SEM image of the fabricated structure. Computed distribution of
the $E_{z}$ field at (c) $\lambda_{1}=4.43\,\mathrm{\mu m}$ and (d)
$\lambda_{2}=5.12\,\mathrm{\mu m}$.

It is possible to combine even more under-coupled nano Fabry-Perot inside the
same sub-wavelength period. As just seen, in Fig 3.a, the B single FP can be
coupled to two different FP (namely A and C), leading to two bi-FP resonators
with different resonance wavelengths. Therefore, three different ribbons can
be placed in the same subwavelength period giving birth to a tri-FP structure,
as illustrated in Fig 4a. This structure has been fabricated on the same
sample than the previously studied resonators, so the individual FP cavities
would still be under-coupled. A top view SEM image of the tri-FP is shown in
the inset of Fig 4b. The computed and measured reflectivity spectra of the ABC
tri-FP structure are plotted in Fig 4b. The experimental characterization of
the resonator reflectivity shows two resonances at
$\lambda_{1}=4.43\,\mathrm{\mu m}$ and $\lambda_{2}=5.12\,\mathrm{\mu m}$ with
a lower-than-$20\%$ reflectivity in both cases and quality factor $Q_{1}=15$
and $Q_{2}=10$. These quality factors are lower than expected theoretically
for a plane wave excitation ($Q_{1}=28$ and $Q_{2}=22$), which is due to the
focusing of the illumination. In order to reach a good agreement between the
experiment and the computation, a $\pm 6^{\circ}$ aperture of the focused beam
inside the FTIR has to be taken into account. The remaining 20% reflectivity
difference between 4.5 $\,\mathrm{\mu m}$ and 5 $\,\mathrm{\mu m}$ can be
attributed to the SiC absorption. The ABC tri-FP spectrum is not the mere
combination of the two bi-FP (AB and BC) responses of Fig 3. Indeed, the two
resonances are red-shifted with comparison to their bi-FP counterparts. It can
be explained by looking at the computed cartography of the field at the two
resonances (Fig 4.c,d). We can locate on the two maps the field distribution
of the anti-bonding modes resonances of AB or BC. They display a stronger
field amplitude unbalance than the corresponding resonances of Fig 2.k. It can
be inferred from this comparison that the presence of the third resonator
affects the anti-symmetry of the field distribution at resonance, thus further
"lightening" the initially dark mode.

In the following, we explore the possibility of adjusting the Q-factor to
achieve higher values. Fig 5 shows the calculated evolution of the Q-factor of
single-FP and bi-FP resonances as a function of the dielectric layer
thickness. The bi-FP resonance exists for a thickness above
180$\,\mathrm{nm}$, at which point the individual FP are undercoupled.From
this minimum thickness, the Q-factor of the bi-FP resonance increases linearly
up to values as high as 85 for a $500\,\mathrm{nm}$-thick SiC layer. The
experimental Q-factors measured on the different nanostructures are reported
on the Fig 5a and are in good agreement with the numerical simulation, the
small difference being here again due to the aperture of the focused incident
beam inside the FTIR. Fig 5b shows the example of one bi-FP and one tri-FP
structure for a $h=500\,\mathrm{nm}$ thick SiC layer. Compared to the previous
structures, the ribbons widths had to be re-evaluated to 450 nm and 610 nm for
the bi-FP and 260 nm, 470 nm and 630 nm for the tri-FP in order to maintain
the resonances around $4.5\,\mathrm{\mu m}$ and $5\,\mathrm{\mu m}$. Moreover,
since this property rely only on the insulator thickness while the resonance
wavelength is mostly driven by the ribbons width, a resonator is easily design
to exhibit a resonance at the desired wavelength and Q-factor.

Figure 5:  (a) Computed evolution of the Q-factor with regards to the SiC
layer thickness for the AB bi-FP (red) and B single-FP (blue) structures. The
measured Q-factors of the AB, BC bi-FP structures and the A,B and C single FP-
structures,reported on the graphic, are in good agreement to the simulation
and prove the high Q-factor accessible for the bi-FP structure. (b)
Reflectivity spectrum of a bi-FP ($W_{A}=450\,\mathrm{nm}$ and
$W_{B}=610\,\mathrm{nm}$) and a tri-FP ($W_{A}=260\,\mathrm{nm}$,
$W_{B}=470\,\mathrm{nm}$ and $W_{C}=630\,\mathrm{nm}$) with a higher SiC
thickness ($h=500\,\mathrm{nm}$). As expected, the quality factor of the
resonance is about 85 for all the resonances.

To conclude, we have demonstrated that the engineering of the coupling between
two nano FP-like nanotructures when using symmetry breaking enables the
excitation of a high Q-factor anti-bonding "dark-mode". We have confirmed
those predictions experimentally and observed anti-bonding resonances with
Q-factors as high as 20 for a set of resonators designed to provide various
resonant wavelengths near 4.5 microns. Moreover, we have demonstrated that
nano-FP like structures can be combined within a same subwavelength period to
provide several anti-bonding resonances. Finally, additional computation show
that the structure Q-factor is tunable through the thickness of the insulating
layer and can reach values up to 85 for both bi-FP and tri-FP structures.
Thanks to their simple fabrication process and their high Q-factor compared to
other MIM structures, we expect these nano-resonators to be useful for thermal
emission, chemical sensing or low dimension spectral filter.

Method The numerical simulation of these nanostructures are performed with the
B-spline modal method which allows fast and exact electromagnetic analysis on
a non-uniform mesh.Bouchon  _et al._ (2010). The gold is modeled using a Drude
formulaPalik and Ghosh (1985) :
$\epsilon\left(\lambda\right)=1-\left[(\lambda_{p}/\lambda+i\gamma)~{}~{}\lambda_{p}/\lambda\right]^{-1}$
with $\lambda_{p}=158.9\,\mathrm{nm}$ and $\gamma=0.0077$ while the SiC is
modeled with a real permitivity $\epsilon=7.29$. The complex electromagnetic
modes of the structures have been calculated using the RCWA method Hugonin and
Lalanne (2005). The nanostructures have been fabricated on a 2mm x 2 mm
footprint area on the same sample. The bottom continuous gold layer (200 nm
thick) was deposited by electron-beam evaporation on a silicon wafer. Then, a
280 nm thick layer of SiC was deposited by evaporation. The patterns were
realized by electron-beam lithography on a Ma-N 2403 resist follow by the
deposition of 2 nm of titanium and 50 nm of gold. Finally a liftoff was
performed using acetone. Using the Scanning Electron Microscope (SEM), we
measured the dimensions : $P_{2}=1770\,\mathrm{nm}$ and
$\mathrm{W_{A}}$=400$\,\mathrm{nm}$,$\mathrm{W_{B}}$=495$\,\mathrm{nm}$,
$\mathrm{W_{C}}$=600$\,\mathrm{nm}$. The ABC tri-FP structure geometrical
parameters have been measured using the SEM as :
$\mathrm{W_{A}}=405\,\mathrm{nm}$, $\mathrm{W_{B}}=500\,\mathrm{nm}$,
$\mathrm{W_{C}}=600\,\mathrm{nm}$. The optical characterizations were achieved
with a Bruker Vertex 70V Fourier transform infrared (FTIR) spectrometer. This
research was supported by a DGA-MRIS scholarship.

## References

  * Cui  _et al._ (2014) Y. Cui, Y. He, Y. Jin, F. Ding, L. Yang, Y. Ye, S. Zhong, Y. Lin, and S. He, “Plasmonic and metamaterial structures as electromagnetic absorbers,” Laser & Photonics Reviews 8, 495–520 (2014).
  * Moreau  _et al._ (2012) A. Moreau, C. Ciracì, J. J. Mock, R. T. Hill, Q. Wang, B. J. Wiley, A. Chilkoti, and D. R. Smith, “Controlled-reflectance surfaces with film-coupled colloidal nanoantennas,” Nature 492, 86–89 (2012).
  * Cattoni  _et al._ (2011) A. Cattoni, P. Ghenuche, A. Haghiri-Gosnet, D. Decanini, J. Chen, J. Pelouard, and S. Collin, “$\lambda$3/1000 plasmonic nanocavities for biosensing fabricated by soft uv nanoimprint lithography,” Nano Letters 11, 3557–3563 (2011).
  * Palaferri  _et al._ (2018) D. Palaferri, Y. Todorov, A. Bigioli, A. Mottaghizadeh, D. Gacemi, A. Calabrese, A. Vasanelli, L. Li, A. G. Davies, E. H. Linfield, _et al._ , “Room-temperature nine-$\mu$m-wavelength photodetectors and ghz-frequency heterodyne receivers,” Nature 556, 85–88 (2018).
  * Makhsiyan  _et al._ (2015) M. Makhsiyan, P. Bouchon, J. Jaeck, J.-L. Pelouard, and R. Haïdar, “Shaping the spatial and spectral emissivity at the diffraction limit,” Applied Physics Letters 107, 251103 (2015).
  * Khurgin (2015) J. B. Khurgin, “How to deal with the loss in plasmonics and metamaterials,” Nature nanotechnology 10, 2–6 (2015).
  * Jain and El-Sayed (2010) P. K. Jain and M. A. El-Sayed, “Plasmonic coupling in noble metal nanostructures,” Chemical Physics Letters 487, 153 – 164 (2010).
  * Deng  _et al._ (2018) T.-S. Deng, J. Parker, Y. Yifat, N. Shepherd, and N. F. Scherer, “Dark plasmon modes in symmetric gold nanoparticle dimers illuminated by focused cylindrical vector beams,” The Journal of Physical Chemistry C 122, 27662–27672 (2018), https://doi.org/10.1021/acs.jpcc.8b10415 .
  * Sheikholeslami  _et al._ (2010) S. Sheikholeslami, Y.-w. Jun, P. K. Jain, and A. P. Alivisatos, “Coupling of optical resonances in a compositionally asymmetric plasmonic nanoparticle dimer,” Nano Letters 10, 2655–2660 (2010), pMID: 20536212, https://doi.org/10.1021/nl101380f .
  * Gao  _et al._ (2018) Y. Gao, N. Zhou, Z. Shi, X. Guo, and L. Tong, “Dark dimer mode excitation and strong coupling with a nanorod dipole,” Photon. Res. 6, 887–892 (2018).
  * Osberg  _et al._ (2014) K. D. Osberg, N. Harris, T. Ozel, J. C. Ku, G. C. Schatz, and C. A. Mirkin, “Systematic study of antibonding modes in gold nanorod dimers and trimers,” Nano Letters 14, 6949–6954 (2014), pMID: 25411044, https://doi.org/10.1021/nl503207j .
  * Christ  _et al._ (2008) A. Christ, O. J. F. Martin, Y. Ekinci, N. A. Gippius, and S. G. Tikhodeev, “Symmetry breaking in a plasmonic metamaterial at optical wavelength,” Nano Letters 8, 2171–2175 (2008), pMID: 18578551, https://doi.org/10.1021/nl0805559 .
  * Zhou and Odom (2011) W. Zhou and T. W. Odom, “Tunable subradiant lattice plasmons by out-of-plane dipolar interactions,” Nature Nanotechnology 6, 423–427 (2011).
  * Omaghali  _et al._ (2014) N. E. J. Omaghali, V. Tkachenko, A. Andreone, and G. Abbate, “Optical sensing using dark mode excitation in an asymmetric dimer metamaterial,” Sensors 14, 272–282 (2014).
  * (15) G. Dayal, X. Y. Chin, C. Soci, and R. Singh, “High-q whispering-gallery-mode-based plasmonic fano resonances in coupled metallic metasurfaces at near infrared frequencies,”  4, 1295–1301.
  * Al-Naib (2020) I. Al-Naib, “Novel terahertz metasurfaces based on complementary coupled split ring resonators,” Optical Materials 99, 109596 (2020).
  * Huang  _et al._ (2010) J.-S. Huang, J. Kern, P. Geisler, P. Weinmann, M. Kamp, A. Forchel, P. Biagioni, and B. Hecht, “Mode imaging and selection in strongly coupled nanoantennas,” Nano Letters 10, 2105–2110 (2010), pMID: 20411912, https://doi.org/10.1021/nl100614p .
  * Yang  _et al._ (2010) S.-C. Yang, H. Kobori, C.-L. He, M.-H. Lin, H.-Y. Chen, C. Li, M. Kanehara, T. Teranishi, and S. Gwo, “Plasmon hybridization in individual gold nanocrystal dimers: Direct observation of bright and dark modes,” Nano Letters 10, 632–637 (2010), pMID: 20058898, https://doi.org/10.1021/nl903693v .
  * Gómez  _et al._ (2013) D. E. Gómez, Z. Q. Teo, M. Altissimo, T. J. Davis, S. Earl, and A. Roberts, “The dark side of plasmonics,” Nano Letters 13, 3722–3728 (2013), pMID: 23802620, https://doi.org/10.1021/nl401656e .
  * Jeyaram  _et al._ (2013) Y. Jeyaram, N. Verellen, X. Zheng, A. Silhanek, M. Hojeij, B. Terhalle, Y. Ekinci, V. Valev, G. Vandenbosch, and V. Moshchalkov, “Rendering dark modes bright by using asymmetric split ring resonators,” Optics express 21, 15464–15474 (2013).
  * Cao  _et al._ (2012) W. Cao, R. Singh, I. A. I. Al-Naib, M. He, A. J. Taylor, and W. Zhang, “Low-loss ultra-high-q dark mode plasmonic fano metamaterials,” Opt. Lett. 37, 3366–3368 (2012).
  * Cui  _et al._ (2011) Y. Cui, J. Xu, K. Hung Fung, Y. Jin, A. Kumar, S. He, and N. X. Fang, “A thin film broadband absorber based on multi-sized nanoantennas,” Applied Physics Letters 99, 253101 (2011).
  * Liu  _et al._ (2011) X. Liu, T. Tyler, T. Starr, A. F. Starr, N. M. Jokerst, and W. J. Padilla, “Taming the blackbody with infrared metamaterials as selective thermal emitters,” Physical review letters 107, 045901 (2011).
  * Bouchon  _et al._ (2012) P. Bouchon, C. Koechlin, F. Pardo, R. Haïdar, and J.-L. Pelouard, “Wideband omnidirectional infrared absorber with a patchwork of plasmonic nanoantennas,” Optics Letters 37, 1038–1040 (2012).
  * Fix  _et al._ (2017) B. Fix, J. Jaeck, P. Bouchon, S. Héron, B. Vest, and R. Haïdar, “High-quality-factor double fabry-perot plasmonic nanoresonator,” Opt. Lett. 42, 5062–5065 (2017).
  * Koechlin  _et al._ (2013) C. Koechlin, P. Bouchon, F. Pardo, J.-L. Pelouard, and R. Haïdar, “Analytical description of subwavelength plasmonic mim resonators and of their combination,” Optics Express 21, 7025–7032 (2013).
  * Seok  _et al._ (2011) T. J. Seok, A. Jamshidi, M. Kim, S. Dhuey, A. Lakhani, H. Choo, P. J. Schuck, S. Cabrini, A. M. Schwartzberg, J. Bokor, E. Yablonovitch, and M. C. Wu, “Radiation engineering of optical antennas for maximum field enhancement,” Nano Letters 11, 2606–2610 (2011), pMID: 21648393, https://doi.org/10.1021/nl2010862 .
  * Bouchon  _et al._ (2010) P. Bouchon, F. Pardo, R. Haïdar, and J. Pelouard, “Fast modal method for subwavelength gratings based on B-splines formulation,” Journal of the Optical Society of America A 27, 696–702 (2010).
  * Palik and Ghosh (1985) E. Palik and G. Ghosh,  _Handbook of optical constants of solids_ (Academic press, 1985).
  * Hugonin and Lalanne (2005) J. Hugonin and P. Lalanne, Reticolo software for grating analysis, Institut d’Optique, Palaiseau, France (2005).