# Multi-task problems are not multi-objective

Michael Ruchte
Representation Learning Lab
University of Freiburg
Freiburg, Germany
<EMAIL_ADDRESS>
&Josif Grabocka
Representation Learning Lab
University of Freiburg
Freiburg, Germany
<EMAIL_ADDRESS>

###### Abstract

Multi-objective optimization (MOO) aims at finding a set of optimal
configurations for a given set of objectives. A recent line of work applies
MOO methods to the typical Machine Learning (ML) setting, which becomes multi-
objective if a model should optimize more than one objective, for instance in
fair machine learning. These works also use Multi-Task Learning (MTL) problems
to benchmark MOO algorithms treating each task as independent objective.

In this work we show that MTL problems do not resemble the characteristics of
MOO problems. In particular, MTL losses are not competing in case of a
sufficiently expressive single model. As a consequence, a single model can
perform just as well as optimizing all objectives with independent models,
rendering MOO inapplicable. We provide evidence with extensive experiments on
the widely used Multi-Fashion-MNIST datasets. Our results call for new
benchmarks to evaluate MOO algorithms for ML. Our code is available at:
https://github.com/ruchtem/moo-mtl.

##  1 Introduction

Multi-Objective Optimization (MOO) is a large field with applications ranging
from car design to Adaptive Bitrate Control (Tanabe and Ishibuchi, 2020; Yang
et al., 2001; Mao et al., 2017). It is defined as finding the optimal trade-
offs between $J$ unweighted objectives:111In this work we ignore further
information like constraints on $x$ or $f_{j}(x)$ or some predefined ordering
of objectives.

$\displaystyle\min_{x\in\mathcal{X}}\left(f_{1}(x),f_{2}(x),\ldots,f_{J}(x)\right)$
(1)

With recent advances in gradient-based MOO (Fliege and Svaiter, 2000;
Désidéri, 2012) an emerging line of deep learning research focuses on applying
MOO algorithms for deep neural networks (Sener and Koltun, 2018; Lin et al.,
2019a; Navon et al., 2020). These works include casting Multi-Task Learning
(MTL) problems into the MOO setting or balancing accuracy and fairness thereby
optimizing neural networks for more than one loss. Different to fair machine
learning, MTL problems come with the availability of several large-scale
datasets and advanced tasks and thus is often used to benchmark MOO approaches
(Lin et al., 2019a; Navon et al., 2020; Ruchte and Grabocka, 2021).

In this work we demonstrate that MTL problems differ from MOO problems in a
significant way. In particular, the objectives are not necessarily competitive
as the network has dedicated parameters for every objective. As a result we
can optimize all objectives with just one model (one solution) rather than
trading off different objectives which renders MOO inapplicable. In fact, we
show that competition is caused by limited model capacity offering the simple
solution to just increase the parameters to solve MOO problems.

To confirm our claim we make extensive experiments on the frequently used
Multi-MNIST datasets (Lin et al., 2019a).

##  2 Related work

The foundations for gradient-based MOO were created by (Fliege and Svaiter,
2000) and (Désidéri, 2012). Particularly, Désidéri (2012) derived the Multiple
Gradient Descent Algorithm (MGDA) which makes a step towards the minimum norm
of the gradients with respect to each objective function, guaranteeing a
Pareto-optimal solution. Sener and Koltun (2018) applied this to MTL problems
and provided an approximate solution for neural networks which come with a
large parameter space. As they demonstrate, applying MOO methods to fit a
network for multiple tasks can be a great advantage. This is not to be
confused with our claim which questions benchmarking MOO algorithms on MTL
problems.

The first work to generate a full Pareto front was Lin et al. (2019a). They
proposed to train $k$ models independently, each bound to a region of the
objective space, naming their method Pareto Multi-Task Learning (PMTL). The
regions are defined by preference rays which reflect a practitioner’s trade-
off choices. Two works use hypernetworks (Ha et al., 2016) to condition the
parameters on any possible preference ray (Navon et al., 2020; Lin et al.,
2020). This way they can approximate the full Pareto front at inference time.
We consider Pareto Hypernetworks (PHN) (Navon et al., 2020) as baseline as its
code is available. Ruchte and Grabocka (2021) follow a similar idea but
condition the model in input feature space by concatenating randomly sampled
preference rays to the input features and add a penalty to ensure a well-
spread Pareto front.

All methods rely strongly on MTL problems to demonstrate their effectiveness,
mainly: Multi-Fashion-MNIST (Lin et al., 2019a), NYUv2 (Nathan Silberman and
Fergus, 2012), and CityScapes (Cordts et al., 2016). Navon et al. (2020) and
Ruchte and Grabocka (2021) additionally evaluate on fairness problems (e.g.
Adult (Kohavi et al., 1996)).

##  3 Multi-task problems are not multi-objective

In this section we compare the MTL setting with original MOO setting. For MTL
the model parameters $\theta$ depend not only on the objective functions
$\mathcal{L}$ but also on mini-batches of training data $D$. Thus with a
slight change towards common machine learning notation Equation 1 becomes

$\displaystyle\min_{\theta\in\Theta}\left(\mathcal{L}_{1}(\theta|D),\mathcal{L}_{2}(\theta|D),\ldots,\mathcal{L}_{J}(\theta|D)\right)$
(2)

This renders MOO a stochastic problem with (generally) non-convex loss
functions resulting in approximate solutions with many different local optima.
Maybe somewhat more challenging but not fundamentally different to Equation 1.
This is the problem to be solved in e.g. fair machine learning where a model
should be both accurate but also non-discriminative with respect to some
attributes like race or gender. Another example would be recommender systems
which should optimize both customer experience and revenue (Lin et al.,
2019b).

For MTL, however, each loss is independent of the others as each loss is tied
to a specific output of the model. Hence a perfectly valid though naive
approach is to optimize a different model for each output (Equation 3).

$\displaystyle\min_{\theta_{1}\in\Theta_{1},\theta_{2}\in\Theta_{2},\ldots,\theta_{J}\in\Theta_{J}}\left(\mathcal{L}_{1}(\theta_{1}|D),\mathcal{L}_{2}(\theta_{2}|D),\ldots,\mathcal{L}_{J}(\theta_{J}|D)\right)$
(3)

For efficiency reasons, however, it is common to share some parts of the
parameters across tasks (Ruder, 2017):

$\displaystyle\min_{\theta_{s}\in\Theta_{s},\theta_{1}\in\Theta_{1},\ldots,\theta_{J}\in\Theta_{J}}(\mathcal{L}_{1}(\theta_{s},\theta_{1}|D),\mathcal{L}_{2}(\theta_{s},\theta_{2}|D),\ldots,\mathcal{L}_{J}(\theta_{s},\theta_{J}|D))$
(4)

where $\theta_{s}$ are shared parameters across all tasks and $\theta_{j}$ are
similar to Equation 3 task-specific parameters. This is different to Equation
1 and Equation 2, where there exists no objective-specific part of $x$ or
$\theta$. With enough capacity in either $\theta_{s}$ or $\theta_{j}$, the
model should be able to optimize all objectives simultaneously as in Equation
3.

##  4 Experiments

We consider training on a single task (Single Task method) to obtain the
optimal score on any given MTL problem (i.e. Equation 3). Additionally, we
benchmark six methods: Uniform which minimizes the uniformly weighted sum of
all objectives $\min\sum^{J}_{j=1}\frac{1}{J}\mathcal{L}_{j}$; MGDA (Sener and
Koltun, 2018); PMTL (Lin et al., 2019a); PHN (linear scalarization if not
mentioned otherwise) (Navon et al., 2020); and COSMOS (Ruchte and Grabocka,
2021).

###  4.1 Experimental setup

We follow the experimental setup of the baselines with fixed 100 training
epochs, 256 batch size, and Adam (Kingma and Ba, 2014). As reference point for
Hypervolume (HV) (Zitzler et al., 2007) we use $(1,1)$. As model we use LeNet
following Lin et al. (2019a). We use the Multi-Fashion-MNIST dataset (Lin et
al., 2019a) with predefined splits: 108k training, 12k validation, and 20k
test images. In these datasets two instances are merged into one image with a
large overlap. The two tasks are to predict the class of each instance. We use
the validation images only for HPO and do not merge them for final training.
We report final results averaged over ten random seeds. We perform all
experiments on Nividia GeForce RTX 2080 GPUs.

Parameter | Range
---|---
Learning rate |  ${1}{\times}10^{-2},{7.5}{\times}10^{-3}$, ${5}{\times}10^{-3},{2.5}{\times}10^{-3}$, $\ldots,{1}{\times}10^{-4}$
Weight decay |  ${1}{\times}10^{-1},{7.5}{\times}10^{-2}$, ${5}{\times}10^{-2},{2.5}{\times}10^{-2}$, $\ldots,{1}{\times}10^{-4}$
Learning rate scheduler |  Cosine Annealing, Step-wise, None
Table 1: Method-independent hyperparameter search spaces

#### Hyperparameter optimization.

We perform random search over the hyperparameter search spaces defined in
Table 1. Learning rate scheduler specifics are ${1}{\times}10^{-6}$ as minimum
learning rate for Cosine Annealing (Loshchilov and Hutter, 2016) and Step-wise
with milestones at epoch 33 and 66 with multiplier 0.1. We sample 100
configurations uniformly at random, which covers approximately 30% of the
$9\times 13\times 3=351$ possible combinations. We select the best
configuration based on Hypervolume based on cross-entropy loss measured on the
validation set.

Some methods come with their own set of hyperparameters. For MGDA this is
gradient normalization
$\in\\{\text{l2},\text{loss},\text{loss+},\text{None}\\}$. For PHN this is the
internal solver $\in\\{\text{ls},\text{EPO}\\}$ and the Dirichlet sampling
parameter $\alpha\in\\{0.1,0.2,0.5,1,1.2,1.5\\}$. For COSMOS this is $\alpha$
and the cosine penalty weight $\lambda\in\\{1,2,4,8,16\\}$. The ranges are
chosen to cover the original values plus one value beyond. Method-specific HPO
covers a smaller part of the search space, approximately 7, 2.37, and 0.95
Percent for MGDA, PHN, and COSMOS, respectively. We denote these methods with
an asterix (∗). For PHN and COSMOS we provide also results not optimizing the
internal solver, $\alpha$, or $\lambda$ taking the original values reported by
the authors. See Section A in the supplementary for the parameters found by
HPO.

###  4.2 Results

Figure 1: Results on Multi-MNIST dataset for different methods. Linear scalarization with uniform weights (Uniform) is a competetive baseline despite its simplicity given its hyperparameters are properly tuned. Method | Params | Multi MNIST | Multi Fashion | multi fashion+mnist
---|---|---|---|---
HV $\uparrow$ |  $\Delta$ ST $\downarrow$ | HV $\uparrow$ |  $\Delta$ ST $\downarrow$ | HV $\uparrow$ |  $\Delta$ ST $\downarrow$
PMTL | 211k | 0.8509 $\pm$ 0.0023 | 0.0058 | 0.7139 $\pm$ 0.0056 | 0.0048 | 0.8220 $\pm$ 0.0062 | 0.0254
PHN∗ | 3,243k | 0.8337 $\pm$ 0.0041 | 0.0229 | 0.6769 $\pm$ 0.0143 | 0.0417 | 0.7996 $\pm$ 0.0053 | 0.0478
PHN | 0.8368 $\pm$ 0.0046 | 0.0198 | 0.6769 $\pm$ 0.0143 | 0.0417 | 0.8073 $\pm$ 0.0051 | 0.0401
COSMOS∗ | 42k | 0.8467 $\pm$ 0.0039 | 0.0099 | 0.7023 $\pm$ 0.0087 | 0.0164 | 0.8132 $\pm$ 0.0026 | 0.0342
COSMOS | 0.8396 $\pm$ 0.0047 | 0.0170 | 0.7064 $\pm$ 0.0055 | 0.0123 | 0.7915 $\pm$ 0.0111 | 0.0559
MGDA∗ | 42k | 0.8482 $\pm$ 0.0060 | 0.0084 | 0.7091 $\pm$ 0.0040 | 0.0096 | 0.8222 $\pm$ 0.0046 | 0.0252
Uniform | 42k | 0.8484 $\pm$ 0.0091 | 0.0083 | 0.7065 $\pm$ 0.0111 | 0.0121 | 0.8189 $\pm$ 0.0067 | 0.0285
Single Task | 42k | 0.8566 $\pm$ 0.0023 | - | 0.7187 $\pm$ 0.0036 | - | 0.8474 $\pm$ 0.0028 | -
Table 2: Hypervolume (HV) based on misclassification rates of properly
optimized methods for the Multi-Mnist datasets (larger is better). $\Delta$ ST
denotes the distance of each method to Single Task, an indication whether the
problem is multi-objective. ∗ denotes optimizing method-specific
hyperparameters as well. Italic denotes the second-best method after Single
Task. Linear scalarization with uniform weights (Uniform) is a competetive
baseline despite its simplicity given its hyperparameters are properly tuned.
Surprisingly, almost all Pareto front generating models (upper section) are
not competitive even though their hyperparameters are also tuned to maximize
HV.

We show the results in terms of Misclassification Rate (MCR) in Figure 1 and
Table 2. We can see that Uniform despite being by far the most simple method
is competitive compared to the more complex methods and achieves nearly the
same performance as Single Task. The gap between Uniform and Single Task is
larger for Multi-Fashion indicating that here the two tasks are indeed
competing. We hypothesis, however, this is due to the limited capacity of the
model and results with different model size in Section 4.3. For scores on each
task and cross-entropy results see Section B in the supplementary.

Unsurprisingly, using HPO improves performance of all methods when comparing
to the numbers reported by the authors. For instance, PMTL reports Single Task
MCR on Multi-MNIST Top-Left of 0.91 and Bottom-Right of 0.88, whereas we get
0.93 and 0.91, respectively. Same is true for the other datasets and methods.
This underlines the importance of HPO for proper method evaluation.
Particularly, surpassing Single Task models in MTL settings should raise
concerns regarding overfitting and thus unfair comparisons.

Somewhat surprising is that even methods which can generate a Pareto front
(especially PHN and COSMOS), do not produce one when their learning rate,
weight decay and scheduler are optimized for validation HV. This problem was
already mentioned by Ruchte and Grabocka (2021) and their main motivation to
introduce the cosine penalty. We double checked and there were HPO
configurations sampled which resulted in a Pareto front but actually achieved
a much lower HV score. We see this as a further indication that we are not
looking at an MOO problem.

###  4.3 Ablations

Dataset | Single Task HV | Uniform HV
---|---|---
MM | 0.8551 $\pm$ 0.0033 | 0.8488 $\pm$ 0.0056
MF | 0.7234 $\pm$ 0.0028 | 0.7122 $\pm$ 0.0075
MFM | 0.8442 $\pm$ 0.0034 | 0.8200 $\pm$ 0.0068
Table 3: Grid seach HV results.

#### Grid search.

To quantify a possible negative impact of random search in HPO we also perform
grid search over all possible 351 hyperparameter combinations and report the
results for a subset of the methods in Table 3. Datasets are abbreviated as MM
(Multi-MNIST), MF (Multi-Fashion), and MFM (Multi-Fashion-MNIST). There are
almost no significant differences when comparing to Table 2 which demonstrates
that randomly searching about one third of the search space is a valid
approach in our case.

Figure 2: Results for Single Task and Uniform when varying the capacity of the
model. Errorbars denote standard deviation. We can see that the distance
between Single Task and Uniform (black arrows) as well as the absolute error
decrease significantly when increasing the capacity of the model.

#### Model capacity.

We ask the question whether capacity plays a crucial role in rendering MTL
problems multi-objective with tasks competing against each other. For this we
alter the capacity of the model by changing the number of channels with a
multiplier. Specifically, we change the model to have $10c$ channels in the
first layer and $20c$ channels in the second layer where
$c\in\\{0.5,1,10,50\\}$ is the channel multiplier (i.e. we only alter
$\theta_{s}$).

The results are presented in Figure 2 and Table 4. We can see that the
distance between Single Task and Uniform scores indeed depends on the capacity
of the models. Reducing capacity increases the distance, increasing capacity
decreases the distance up to a saturation point where error rates almost
match. Interestingly, for Multi-Fashion-MNIST there is still a considerable
gap between single task and uniform performance. This could be an indication
that not only capacity matters but also the similarity of the tasks.

The experiments support the hypothesis that for MTL problems it is more
beneficial to increase the capacity instead of treating it as a MOO problem
trying to balance presumed conflicts between losses. This is especially true
for methods which inherently increase the capacity as PHN does. Multiplying
the channels by 10 and using uniform scalarization yields a better Hypervolume
with about 25% of the parameters compared to PHN (0.03 HV improve on Multi-
MNIST).

Method | $c$ | Params | Multi Mnist | Multi fashion | Multi fashion+mnist
---|---|---|---|---|---
HV |  $\Delta$ ST | HV |  $\Delta$ ST | HV |  $\Delta$ ST
Single Task | 0.5x | 20k | 0.8340 $\pm$ 0.0028 | - | 0.6978 $\pm$ 0.0045 | - | 0.8313 $\pm$ 0.0025 | -
Uniform | 0.8238 $\pm$ 0.0065 | 0.0102 | 0.6767 $\pm$ 0.0143 | 0.0212 | 0.7885 $\pm$ 0.0095 | 0.0428
Single Task | 1x | 42k | 0.8566 $\pm$ 0.0023 | - | 0.7187 $\pm$ 0.0036 | - | 0.8474 $\pm$ 0.0028 | -
Uniform | 0.8484 $\pm$ 0.0091 | 0.0083 | 0.7065 $\pm$ 0.0111 | 0.0121 | 0.8189 $\pm$ 0.0067 | 0.0285
Single Task | 10x | 863k | 0.8677 $\pm$ 0.0022 | - | 0.7355 $\pm$ 0.0042 | - | 0.8576 $\pm$ 0.0025 | -
Uniform | 0.8661 $\pm$ 0.0023 | 0.0016 | 0.7288 $\pm$ 0.0107 | 0.0067 | 0.8391 $\pm$ 0.0026 | 0.0184
Single Task | 50x | 14,315k | 0.8680 $\pm$ 0.0012 | - | 0.7384 $\pm$ 0.0033 | - | 0.8571 $\pm$ 0.0023 | -
Uniform | 0.8667 $\pm$ 0.0032 | 0.0013 | 0.7317 $\pm$ 0.0076 | 0.0067 | 0.8341 $\pm$ 0.0042 | 0.0230
Table 4: Ablation results for different model sizes. $c$ denotes the channel
multiplier, $\Delta$ ST denotes the distance of each method to Single Task.
The results show that $\Delta$ ST decreases on overall when increasing the
model capacity.

##  5 Limitations

One limitation is that we did not do proper HPO for methods with additional
hyperparameters as this would result in unbearable computing costs. For the
same reason we also used only one seed for HPO, despite the results in Section
4.3 being aware that the obtained scores are noisy.

Also due to computing costs we decided to not repeat the analysis with larger
models on one of the larger benchmarks mentioned in Section 2. We, however, do
not see a reason why the results should not generalize to larger problems.
This could be a direction for future research, as well as an analysis on the
effect of altering the capacity of the task-specific heads, i.e. $\theta_{j}$.

##  6 Conclusion

In this work we showed that multi-task problems do not necessarily resemble
the characteristics of multi-objective problems. One important aspect in
rendering MTL tasks multi-objective is the capacity of the model and likely
the similarity between the tasks. The results call for a new set of benchmarks
which are truly multi-objective. Pointers could be Fairness (Navon et al.,
2020; Ruchte and Grabocka, 2021) and recommender systems (Lin et al., 2019b;
Milojkovic et al., 2019). Especially the latter come with large datasets and
require large models which is more suited for deep neural networks. Another
pointer is Deist et al. (2021) who evaluate on medical image segmentation with
possibly conflicting labeling.

It required roughly 3100 GPU hours to compute the results presented here,
where 2780 are due to the ablations.

## Acknowledgments

Michael Ruchte would like to thank Samuel Müller for the helpful discussions.

## References

  * Cordts et al. [2016] Marius Cordts, Mohamed Omran, Sebastian Ramos, Timo Rehfeld, Markus Enzweiler, Rodrigo Benenson, Uwe Franke, Stefan Roth, and Bernt Schiele.  The cityscapes dataset for semantic urban scene understanding.  In _Proc. of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016.
  * Deist et al. [2021] Timo M Deist, Monika Grewal, Frank JWM Dankers, Tanja Alderliesten, and Peter AN Bosman.  Multi-objective learning to predict pareto fronts using hypervolume maximization.  _arXiv preprint arXiv:2102.04523_ , 2021.
  * Désidéri [2012] Jean-Antoine Désidéri.  Mutiple-gradient descent algorithm for multiobjective optimization.  In _European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012)_ , 2012.
  * Fliege and Svaiter [2000] Jörg Fliege and Benar Fux Svaiter.  Steepest descent methods for multicriteria optimization.  _Mathematical methods of operations research_ , 51(3):479–494, 2000.
  * Ha et al. [2016] David Ha, Andrew Dai, and Quoc V Le.  Hypernetworks.  _arXiv preprint arXiv:1609.09106_ , 2016.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kohavi et al. [1996] Ron Kohavi et al.  Scaling up the accuracy of naive-bayes classifiers: A decision-tree hybrid.  In _Kdd_ , volume 96, pages 202–207, 1996.
  * Lin et al. [2019a] Xi Lin, Hui-Ling Zhen, Zhenhua Li, Qing-Fu Zhang, and Sam Kwong.  Pareto multi-task learning.  _Advances in neural information processing systems_ , 32:12060–12070, 2019a.
  * Lin et al. [2020] Xi Lin, Zhiyuan Yang, Qingfu Zhang, and Sam Kwong.  Controllable pareto multi-task learning.  _arXiv preprint arXiv:2010.06313_ , 2020.
  * Lin et al. [2019b] Xiao Lin, Hongjie Chen, Changhua Pei, Fei Sun, Xuanji Xiao, Hanxiao Sun, Yongfeng Zhang, Wenwu Ou, and Peng Jiang.  A pareto-efficient algorithm for multiple objective optimization in e-commerce recommendation.  In _Proceedings of the 13th ACM Conference on recommender systems_ , pages 20–28, 2019b.
  * Loshchilov and Hutter [2016] Ilya Loshchilov and Frank Hutter.  Sgdr: Stochastic gradient descent with warm restarts.  _arXiv preprint arXiv:1608.03983_ , 2016.
  * Mao et al. [2017] Hongzi Mao, Ravi Netravali, and Mohammad Alizadeh.  Neural adaptive video streaming with pensieve.  In _Proceedings of the Conference of the ACM Special Interest Group on Data Communication_ , pages 197–210, 2017.
  * Milojkovic et al. [2019] Nikola Milojkovic, Diego Antognini, Giancarlo Bergamin, Boi Faltings, and Claudiu Musat.  Multi-gradient descent for multi-objective recommender systems.  _arXiv preprint arXiv:2001.00846_ , 2019.
  * Nathan Silberman and Fergus [2012] Pushmeet Kohli Nathan Silberman, Derek Hoiem and Rob Fergus.  Indoor segmentation and support inference from rgbd images.  In _ECCV_ , 2012.
  * Navon et al. [2020] Aviv Navon, Aviv Shamsian, Ethan Fetaya, and Gal Chechik.  Learning the pareto front with hypernetworks.  In _International Conference on Learning Representations_ , 2020.
  * Ruchte and Grabocka [2021] Michael Ruchte and Josif Grabocka.  Scalable pareto front approximation for deep multi-objective learning.  _Proceedings of the IEEE International Conference on Data Mining (ICDM)_ , 2021.
  * Ruder [2017] Sebastian Ruder.  An overview of multi-task learning in deep neural networks.  _arXiv preprint arXiv:1706.05098_ , 2017.
  * Sener and Koltun [2018] Ozan Sener and Vladlen Koltun.  Multi-task learning as multi-objective optimization.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pages 525–536, 2018.
  * Tanabe and Ishibuchi [2020] Ryoji Tanabe and Hisao Ishibuchi.  An easy-to-use real-world multi-objective optimization problem suite.  _Applied Soft Computing_ , 89:106078, 2020.
  * Yang et al. [2001] RJ Yang, N Wang, CH Tho, JP Bobineau, and BP Wang.  Metamodeling development for vehicle frontal impact simulation.  In _International Design Engineering Technical Conferences and Computers and Information in Engineering Conference_ , volume 80227, pages 27–35. American Society of Mechanical Engineers, 2001.
  * Zitzler et al. [2007] Eckart Zitzler, Dimo Brockhoff, and Lothar Thiele.  The hypervolume indicator revisited: On the design of pareto-compliant indicators via weighted integration.  In _International Conference on Evolutionary Multi-Criterion Optimization_ , pages 862–876. Springer, 2007.

## Appendix

##  Appendix A HPO results

We show the hyperparameter values found by hyperparameter optimization for
Section 4.2 in Table 5 and for Section 4.3 in Table 6.

Parameter | Dataset | ST TL | ST BR | Uniform | MGDA | PMTL | PHN | PHN∗ | COSMOS | COSMOS∗
---|---|---|---|---|---|---|---|---|---|---
Lerning rate | MM | 0.0005 | 0.00075 | 0.0025 | 0.001 | 0.00025 | 0.00075 | 0.001 | 0.00075 | 0.001
MF | 0.0025 | 0.001 | 0.00075 | 0.0025 | 0.00075 | 0.00075 | 0.00075 | 0.001 | 0.00075
MFM | 0.0025 | 0.00075 | 0.00075 | 0.001 | 0.00075 | 0.0005 | 0.0025 | 0.0075 | 0.001
Weight decay | MM | 0.05 | 0.05 | 0.0075 | 0.01 | 0.075 | 0.05 | 0.075 | 0.0075 | 0.0075
MF | 0.01 | 0.01 | 0.01 | 0.005 | 0.05 | 0.075 | 0.075 | 0.0075 | 0.0075
MFM | 0.0075 | 0.01 | 0.01 | 0.01 | 0.05 | 0.075 | 0.005 | 0.0005 | 0.01
Learning rate scheduler | MM | CA | MS | CA | CA | CA | CA | CA | CA | CA
MF | MS | CA | CA | CA | CA | CA | CA | CA | CA
MFM | CA | MS | CA | CA | CA | CA | CA | CA | MS
Normal- ization | MM | - | - | - | loss | - | - | - | - | -
MF | - | - | - | loss+ | - | - | - | - | -
MFM | - | - | - | none | - | - | - | - | -
Internal solver | MM | - | - | - | - | - | linear† | epo | - | -
MF | - | - | - | - | - | linear† | linear | - | -
MFM | - | - | - | - | - | linear† | linear | - | -
$\alpha$ | MM | - | - | - | - | - | 0.2† | 0.1 | 1.2† | 0.5
MF | - | - | - | - | - | 0.2† | 0.2 | 1.2† | 1.0
MFM | - | - | - | - | - | 0.2† | 1.2 | 1.2† | 0.2
$\lambda$ | MM | - | - | - | - | - | - | - | 2† | 1
MF | - | - | - | - | - | - | - | 2† | 2
MFM | - | - | - | - | - | - | - | 8† | 1
Table 5: Hyperparameters found by HPO, †denotes the values from the original works, here just for reference. MM=Multi MNIST, MF=Multi Fashion, MFM=Multi Fashion+MNIST. Parameter | Dataset | Single Task TL | Single Task BR | Uniform
---|---|---|---|---
|  | 0.5x | 10x | 50x | 0.5x | 10x | 50x | 0.5x | 10x | 50x
Lerning rate | MM | 0.001 | 0.0001 | 0.0001 | 0.00075 | 0.00025 | 0.00025 | 0.0005 | 0.0005 | 0.0001
MF | 0.001 | 0.00075 | 0.0005 | 0.0025 | 0.001 | 0.0005 | 0.0025 | 0.00075 | 0.0005
MFM | 0.001 | 0.0005 | 0.00075 | 0.0025 | 0.00075 | 0.00025 | 0.0025 | 0.0005 | 0.0005
Weight decay | MM | 0.05 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.0075 | 0.0075 | 0.0075
MF | 0.0075 | 0.01 | 0.01 | 0.005 | 0.0075 | 0.01 | 0.005 | 0.0075 | 0.0075
MFM | 0.01 | 0.01 | 0.01 | 0.0075 | 0.01 | 0.01 | 0.0075 | 0.0075 | 0.005
Learning rate scheduler | MM | MS | MS | MS | MS | MS | CA | CA | MS | MS
MF | MS | MS | MS | MS | CA | CA | CA | CA | MS
MFM | MS | MS | CA | MS | MA | MS | CA | MS | CA
Table 6: Hyperparameters found by HPO for the models with different capacity.
The second row in the header refers to the channel multiplier (see Section
4.3). MM=Multi MNIST, MF=Multi Fashion, MFM=Multi Fashion+MNIST. Learning rate
scheduler abbreviations: MS=Step-wise (multi step), CA=Cosine Annealing.

##  Appendix B More results

We show the performance in terms of Misclassification rate for each task in
Table 7. The performance for each task in terms of cross entropy is shown in
Figure 3 and Table 9. An equivalent of Table 2 for Hypervolume calculated
based on cross entropy is presented in Table 8.

Method | Multi MNIST | Multi Fashion | multi fashion+mnist
---|---|---|---
TL | BR | TL | BR | TL | BR
Single Task | 0.0655 $\pm$ 0.0013 | 0.0834 $\pm$ 0.0023 | 0.1519 $\pm$ 0.0031 | 0.1526 $\pm$ 0.0020 | 0.0376 $\pm$ 0.0015 | 0.1195 $\pm$ 0.0023
Uniform | 0.0704 $\pm$ 0.0051 | 0.0874 $\pm$ 0.0051 | 0.1578 $\pm$ 0.0070 | 0.1611 $\pm$ 0.0066 | 0.0517 $\pm$ 0.0041 | 0.1365 $\pm$ 0.0036
MGDA | 0.0715 $\pm$ 0.0032 | 0.0865 $\pm$ 0.0036 | 0.1579 $\pm$ 0.0038 | 0.1580 $\pm$ 0.0030 | 0.0490 $\pm$ 0.0029 | 0.1354 $\pm$ 0.0035
ParetoMTL | 0.0806 $\pm$ 0.0030 | 0.0907 $\pm$ 0.0023 | 0.1646 $\pm$ 0.0064 | 0.1584 $\pm$ 0.0047 | 0.0732 $\pm$ 0.0044 | 0.1298 $\pm$ 0.0030
PHN∗ | 0.0765 $\pm$ 0.0020 | 0.0972 $\pm$ 0.0031 | 0.1746 $\pm$ 0.0092 | 0.1800 $\pm$ 0.0085 | 0.0583 $\pm$ 0.0038 | 0.1509 $\pm$ 0.0036
PHN | 0.0761 $\pm$ 0.0025 | 0.0943 $\pm$ 0.0031 | 0.1746 $\pm$ 0.0092 | 0.1800 $\pm$ 0.0085 | 0.0518 $\pm$ 0.0041 | 0.1486 $\pm$ 0.0030
COSMOS∗ | 0.0723 $\pm$ 0.0022 | 0.0873 $\pm$ 0.0025 | 0.1605 $\pm$ 0.0051 | 0.1634 $\pm$ 0.0058 | 0.0645 $\pm$ 0.0012 | 0.1308 $\pm$ 0.0023
COSMOS | 0.0767 $\pm$ 0.0028 | 0.0907 $\pm$ 0.0029 | 0.1586 $\pm$ 0.0026 | 0.1604 $\pm$ 0.0047 | 0.1081 $\pm$ 0.0069 | 0.1573 $\pm$ 0.0061
Table 7: Results for each task in terms of Misclassification rate (MRC) Figure 3: Results in terms of cross-entropy loss Method | Multi MNIST | Multi Fashion | multi fashion+mnist
---|---|---|---
HV |  $\Delta$ ST | HV |  $\Delta$ ST | HV |  $\Delta$ ST
Single Task | 0.6024 $\pm$ 0.0048 | 0 | 0.3272 $\pm$ 0.0043 | 0 | 0.5732 $\pm$ 0.0072 | 0
Uniform | 0.5777 $\pm$ 0.0191 | 0.0248 | 0.3116 $\pm$ 0.0202 | 0.0157 | 0.5147 $\pm$ 0.0148 | 0.0584
MGDA∗ | 0.5789 $\pm$ 0.0158 | 0.0235 | 0.3168 $\pm$ 0.0072 | 0.0104 | 0.5254 $\pm$ 0.0102 | 0.0478
PMTL | 0.5888 $\pm$ 0.0080 | 0.0136 | 0.3261 $\pm$ 0.0106 | 0.0011 | 0.5209 $\pm$ 0.0160 | 0.0523
PHN∗ | 0.5438 $\pm$ 0.0097 | 0.0586 | 0.2631 $\pm$ 0.0247 | 0.0641 | 0.4738 $\pm$ 0.0119 | 0.0994
PHN | 0.5539 $\pm$ 0.0110 | 0.0485 | 0.2631 $\pm$ 0.0247 | 0.0641 | 0.4905 $\pm$ 0.0099 | 0.0827
COSMOS∗ | 0.5772 $\pm$ 0.0102 | 0.0252 | 0.3082 $\pm$ 0.0172 | 0.0190 | 0.5011 $\pm$ 0.0069 | 0.0720
COSMOS | 0.5603 $\pm$ 0.0106 | 0.0421 | 0.3159 $\pm$ 0.0118 | 0.0113 | 0.4563 $\pm$ 0.0219 | 0.1169
Table 8: Hypervolume (HV) results based on cross-entropy loss Method | Multi MNIST | Multi Fashion | multi fashion+mnist
---|---|---|---
TL | BR | TL | BR | TL | BR
Single Task | 0.1993 $\pm$ 0.0035 | 0.2476 $\pm$ 0.0055 | 0.4151 $\pm$ 0.0059 | 0.4405 $\pm$ 0.0047 | 0.1219 $\pm$ 0.0050 | 0.3473 $\pm$ 0.0067
Uniform | 0.2174 $\pm$ 0.0139 | 0.2621 $\pm$ 0.0123 | 0.4318 $\pm$ 0.0190 | 0.4523 $\pm$ 0.0186 | 0.1634 $\pm$ 0.0122 | 0.3849 $\pm$ 0.0093
MGDA | 0.2189 $\pm$ 0.0105 | 0.2589 $\pm$ 0.0114 | 0.4316 $\pm$ 0.0079 | 0.4426 $\pm$ 0.0071 | 0.1558 $\pm$ 0.0070 | 0.3776 $\pm$ 0.0087
ParetoMTL | 0.2476 $\pm$ 0.0091 | 0.2706 $\pm$ 0.0062 | 0.4489 $\pm$ 0.0190 | 0.4380 $\pm$ 0.0155 | 0.2375 $\pm$ 0.0127 | 0.3644 $\pm$ 0.0066
PHN∗ | 0.2347 $\pm$ 0.0056 | 0.2895 $\pm$ 0.0079 | 0.4761 $\pm$ 0.0273 | 0.4989 $\pm$ 0.0239 | 0.1871 $\pm$ 0.0106 | 0.4172 $\pm$ 0.0104
PHN | 0.2308 $\pm$ 0.0074 | 0.2799 $\pm$ 0.0084 | 0.4761 $\pm$ 0.0273 | 0.4989 $\pm$ 0.0239 | 0.1681 $\pm$ 0.0102 | 0.4105 $\pm$ 0.0075
COSMOS∗ | 0.2208 $\pm$ 0.0065 | 0.2592 $\pm$ 0.0076 | 0.4360 $\pm$ 0.0144 | 0.4540 $\pm$ 0.0176 | 0.2101 $\pm$ 0.0040 | 0.3656 $\pm$ 0.0063
COSMOS | 0.2342 $\pm$ 0.0073 | 0.2684 $\pm$ 0.0080 | 0.4306 $\pm$ 0.0096 | 0.4453 $\pm$ 0.0127 | 0.3504 $\pm$ 0.0160 | 0.4289 $\pm$ 0.0124
Table 9: Results for each task in terms of cross-entropy loss


footnote  ††thanks: These authors have equally contributed to this work.
††thanks: These authors have equally contributed to this work.

# J-coupling NMR Spectroscopy with Nitrogen Vacancy Centers at High Fields

P. Alsina-Bolívar  Department of Physical Chemistry, University of the Basque
Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain  EHU Quantum Center,
University of the Basque Country UPV/EHU, Leioa, Spain  A. Biteri-Uribarren
Department of Physical Chemistry, University of the Basque Country UPV/EHU,
Apartado 644, 48080 Bilbao, Spain  EHU Quantum Center, University of the
Basque Country UPV/EHU, Leioa, Spain  C. Munuera-Javaloy  Department of
Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644,
48080 Bilbao, Spain  EHU Quantum Center, University of the Basque Country
UPV/EHU, Leioa, Spain  J. Casanova  Department of Physical Chemistry,
University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain
EHU Quantum Center, University of the Basque Country UPV/EHU, Leioa, Spain
IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5, 48009 Bilbao,
Spain

###### Abstract

A diamond-based sensor utilizing nitrogen-vacancy (NV) center ensembles
permits the analysis of micron-sized samples through NMR techniques at room
temperature. Current efforts are directed towards extending the operating
range of NV centers into high magnetic fields, driven by the potential for
larger nuclear spin polarization of the target sample and the presence of
enhanced chemical shifts. Especially interesting is the access to J-couplings
as they carry information of chemical connectivity inside molecules. In this
work, we present a protocol to access J-couplings in both homonuclear and
heteronuclear cases with NV centers at high magnetic fields. Our protocol
leads to a clear spectrum exclusively containing J-coupling features with high
resolution. This resolution is limited primarily by the decoherence of the
target sample, which is mitigated by the noise filtering capacities of our
method.

_Introduction.–_ Nuclear Magnetic Resonance (NMR) spectroscopy has established
itself as a fundamental tool for exploring materials since its inception in
the 1950s Slichter78 ; levitt2013spin ; tampieri2022a . This technique is
essential for analyzing molecular structure and function. It finds extensive
applications in the examination of macroscopic biological tissues and
functional studies in medical imaging chary2008nmr ; bottomley1984nmr ;
abragam1961principles . However, the inherent low sensitivity of standard NMR
spectroscopy restricts its application to bulky samples allert2022advances .
In this scenario, newly developed solid-state quantum sensors, open up new
avenues for spectroscopic measurements with unprecedented sensitivity and
spatial resolution, enabling studies at the cellular neuling2023prospects  and
single-molecular level sushkov2014all ; shi2015single ; lovchinsky2016nuclear
. In particular, quantum sensors consisting on NV ensembles have demonstrated
their utility in conducting NMR spectroscopy on microscopic samples at room
temperature glenn2018high ; bucher2020hyperpolarisation ; arunkumar2021micron
, while the pursuit of extending the operating range of NV-based sensors into
the high magnetic field regime is driven by the advantages offered by elevated
fields. Namely, a larger thermal nuclear spin polarization and enhanced energy
shifts that simplify data extraction.

The ability to access J-couplings – also known as scalar or indirect
couplings– is of particular interest in the actual spectroscopic framework.
J-couplings are indirect intramolecular dipole-dipole interactions mediated
through chemical bonds. They are responsible for complex splittings in NMR
spectral lines, and the examination of J-coupling patterns allows scientists
to obtain information regarding nuclear connectivity within a molecule. Hence,
accessing J-couplings at the microscale regime would open new avenues in mass-
limited scenarios for applications in, e.g., combinatorial chemistry, as well
as in cell and membrane biology liu2017combinatorial ; liang2016nmr .
Nevertheless, J-couplings are not trivial to detect since (i) They are
independent of the externally applied field. Then, at high magnetic fields
other field-dependent terms dominate the dynamics, and (ii) Long coherence
times are required to resolve the spectral peaks corresponding to J-couplings
owing to their weak value (typically, on the order of a few Hertz).

In this Letter, we present J-INSECT (J-coupling Induced Nuclear Signal with
Extended Coherence Time). This protocol enables the detection of homonuclear
and heteronuclear J-couplings with NV ensembles (thus, accessing microscale
samples) in the regime of high magnetic fields. J-INSECT integrates heterodyne
techniques schmitt2017submillihertz ; munuera2023high ; meinel2023high ;
boss2017quantm  and exerts noise filtering over the sample effectively
enhancing its coherence. This results into highly precise estimations of
J-couplings. We investigate the performance of J-INSECT in regimes including
both a reduced and a large number of measurements, and analyze the impact of
distinct noise sources (including the presence of always-on J-couplings)
demonstrating the feasibility of the method in common experimental conditions.

Figure 1: Scheme of J-INSECT with RF and MW controls. Each channel (black bold
horizontal line) denotes the radiation pattern over each system constituent.
The upper RF channel corresponds to the pulses over non-hydrogen nuclear
species. In this respect: to detect heteronuclear J-couplings, simultaneous
$\pi$-pulses over H and another isotope are delivered. The central channel
pictures the RF that must be applied to H nuclei, i.e. an initial
$(\pi/2)_{y}$-pulse, followed by an alternation of $(\pi)_{x}$-pulses and two
full Rabi oscillations along $Y$ for detection, with winding-unwinding
mechanism for robustness as originally proposed in munuera2023high  and
further developed in daly2023nutation . The resultant magnetization of
hydrogen nuclei, $\vec{M}\propto\sum_{i}^{N_{H}}\left[\langle
S_{i}^{x}\rangle,\langle S_{i}^{y}\rangle,\langle S_{i}^{z}\rangle\right]$
suppmat , at different instants is represented with red arrows in the spheres:
dashed arrows indicate the magnetization at the beginning of each protocol-
step (pulses or free evolutions), while the solid arrow indicates the
magnetization at the end of the step. The evolution of the magnetization
outside the xy-plane (induced by RF drivings) is depicted with blue arrows
indicating the magnetization’s trajectory. During the detection stage, the
$Y\bar{Y}$ pulse induces a magnetic field ($B_{z}$) on the NV ensemble along
the NV axis that precess at a controllable speed depending on the RF Rabi
frequency. The amplitude of $B_{z}$ at the $n$-th detection period is
proportional to the amplitude of the red solid arrow (i.e. to the
magnetization at the end of the $n$-th encoding period). The $Y\bar{Y}$ pulses
bring the magnetization back to the state at the start of the detection stage
as they induce two $2\pi$ rotations. $B_{z}$ is captured by the NV ensemble
via a MW control sequence (e.g., the XY4 in the lower channel), while green
pulses indicate the standard initialization and readout of the NVs.

_Molecular target Hamiltonian.–_ To exemplify the method we consider an
isotropic liquid sample comment  containing a target molecular ensemble with
two types of atoms, i.e. each molecule comprises $N_{\text{H}}$ hydrogens (H
from now on) and a number $N_{\text{A}}$ of nuclei of the “A” species. The
resulting Hamiltonian is

$\displaystyle H/\hbar$
$\displaystyle=\omega_{\text{H}}\sum_{i=1}^{N_{\text{H}}}S_{i}^{z}\>+\>\sum_{i=1}^{N_{\text{H}}}\delta_{i}^{\text{H}}\>S_{i}^{z}+\sum_{i<j}^{N_{\text{H}}}J_{i,j}^{\text{H}}\>\vec{S}_{i}\cdot\vec{S}_{j}$
$\displaystyle+\omega_{\text{A}}\sum_{i}^{N_{\text{A}}}\>I_{i}^{z}\>+\>\sum_{i=1}^{N_{\text{A}}}\delta_{i}^{\text{A}}\>I_{i}^{z}+\sum_{i<j}^{N_{\text{A}}}J_{i,j}^{\text{A}}\>\vec{I}_{i}\cdot\vec{I}_{j}$
$\displaystyle+\sum_{i}^{N_{\text{H}}}\sum_{j}^{N_{\text{A}}}J_{i,j}^{\text{H}-\text{A}}\>\vec{S}_{i}\cdot\vec{I}_{j}$
$\displaystyle+\sum_{i}^{N_{\text{H}}}\Omega_{i}^{\text{H}}(t)S_{i}^{x}\cos{\left(\omega_{1}\>t\right)}+\sum_{i}^{N_{\text{A}}}\Omega_{i}^{\text{A}}(t)I_{i}^{x}\cos{(\omega_{2}\>t)},$
(1)

where we denote the spin operators of H nuclei as $S_{k}^{x,y,z}$, and those
of the A species as $I_{k}^{x,y,z}$. On the first line we have: Larmor terms
of the $N_{\text{H}}$ H (with
$\omega_{\text{H}}=|\gamma_{\text{H}}|B_{\text{ext}}$), their chemical shifts,
and the homonuclear J-couplings of magnitude $J_{i,j}^{\rm{H}}$. The terms on
the second line are the equivalent interactions for the A species in the
target molecule. The third line comprises heteronuclear J-couplings between H
and A nuclei ($J_{i,j}^{\text{H}-\text{A}}$), while the last one contains the
RF controls over H and A nuclear spins. Due to significant differences in
gyromagnetic ratios among distinct spin species, crosstalk effects are safely
neglected using the RWA, especially in high-field scenarios. By setting the RF
fields such that $\omega_{1}\equiv\omega_{\text{H}}$ and
$\omega_{2}\equiv\omega_{\text{A}}$, in a rotating frame with respect to
(w.r.t) Larmor terms, Eq. (1) simplifies to

$H/\hbar=H_{shift}+\sum_{\begin{subarray}{c}i<j\\\
\text{neq}\end{subarray}}^{N_{\text{H}}}J_{i,j}^{\text{H}}\>S_{i}^{z}S_{j}^{z}+\sum_{i<j}^{N_{\text{A}}}J_{i,j}^{\text{A}}\>I_{i}^{z}I_{j}^{z}+\sum_{i}^{N_{\text{H}}}\sum_{j}^{N_{\text{A}}}{J}_{i,j}^{\text{H}-\text{A}}\>S_{i}^{z}I_{j}^{z}+H_{c}$
(2)

where
$H_{shift}=\sum_{i=1}^{N_{\text{H}}}\delta_{i}^{\text{H}}\>S_{i}^{z}+\sum_{i=1}^{N_{\text{A}}}\delta_{i}^{\text{A}}\>I_{i}^{z}$
and
$H_{c}=\sum_{i}^{N_{\text{H}}}\Omega_{i}^{\text{H}}S_{i}^{x}+\sum_{i}^{N_{\text{A}}}\Omega_{i}^{\text{A}}I_{i}^{x}$.
For the sake of simplicity in the presentation of the method, Eq. (2) assumes
the elimination of the transversal spin operators ($I_{i}^{x,y},I_{j}^{x,y}$,
$S_{i}^{x,y}$ and $S_{j}^{x,y}$). This applies in an scenario such that
$|\omega_{\text{H}}-\omega_{i}|>>J_{i,j}^{\text{H}-\text{A}}$,
$|\delta^{\text{H}}_{i}-\delta^{\text{H}}_{j}|>>J^{H}_{i,j}$, and
$|\delta^{\text{A}}_{i}-\delta^{\text{A}}_{j}|>>J^{\text{A}}_{i,j}$ (note, the
validity of these conditions is stronger at high magnetic fields). In
addition, the label “neq” in the first summation indicates that the
homonuclear J-couplings among magnetically equivalent nuclei (i.e., those with
identical chemical shifts) are not included in Eq. (2), since they do not have
an impact in the dynamics as confirmed via numerical simulations. As a
summary, Eq. (2) features only ZZ type interactions and serves as a model to
explain the basics of our protocol. We remark that our numerical simulations
consider deviations w.r.t. Eq. (2) due to shortcomings in the application of
prior approximations and the consideration of all J-couplings throughout the
entire protocol. This encompasses instances during the application of finite-
width RF pulses, where deviations in the Rabi frequency are also taken into
account.

_The protocol.–_ J-INSECT is illustrated in Fig. 1. This starts with an
initial $\frac{\pi}{2}$-pulse over the H nuclei that transfer their thermal
polarization to the x-axis, and continues by alternating encoding and
detection stages.

During free evolutions (of length $\tau/2$) at the encoding stages, Eq. (2)
simplifies to

$H/\hbar=H_{shift}+\sum_{\begin{subarray}{c}i<j\\\
\text{neq}\end{subarray}}^{N_{\text{H}}}J_{i,j}^{\text{H}}\>S_{i}^{z}S_{j}^{z}+\sum_{i}^{N_{\text{H}}}\sum_{j}^{N_{\text{A}}}{J}_{i,j}^{\text{H}-\text{A}}\>S_{i}^{z}I_{j}^{z},$
(3)

as $\sum_{i<j}^{N_{\text{A}}}J_{i,j}^{\text{A}}\>I_{i}^{z}I_{j}^{z}$ commutes
with every other term in Eq. (3) and does not impact the dynamics of H nuclei
–the emitters in J-INSECT– (check readout stage). $\pi$-pulses are delivered
in the middle of the encoding stage with a twofold objective. On the one hand,
they cancel $H_{shift}$ mitigating broadening caused by spatial
inhomogeneities in the external magnetic field. In other words, the refocusing
$\pi$-pulses lengthen the nuclear coherence time from $T_{2}^{*}$ to $T_{2}$.
This enables a longer scanning duration of the nuclear sample, which proves
especially advantageous for estimating J-couplings due to their relatively
weak values. On the other hand, $\pi$-pulses over H (and A) enable the
selective encoding of $J_{i,j}^{\text{H}}$ (and $J_{i,j}^{\text{H}-\text{A}}$)
in the system dynamics. This becomes evident upon calculating the effective
propagator at the end of the encoding stage after delivering $\pi$-pulses over
H nuclei, resulting in

$U_{1}(t)=\exp{\left[-it\sum_{\begin{subarray}{c}i<j\\\
\text{neq}\end{subarray}}^{N_{\text{H}}}J_{i,j}^{\text{H}}\>S_{i}^{z}S_{j}^{z}\right]}.$
(4)

Alternatively, in the scenario where synchronized $\pi$-pulses are applied to
H and A, resembling a double nucleus-nucleus resonance, one gets

$\displaystyle U_{2}(t)=\exp{\left[-it\left(\sum_{\begin{subarray}{c}i<j\\\
\text{neq}\end{subarray}}^{N_{\text{H}}}J_{i,j}^{\text{H}}\>S_{i}^{z}S_{j}^{z}+\sum_{i}^{N_{\text{H}}}\sum_{j}^{N_{\text{A}}}J_{i,j}^{\text{H}-\text{A}}\>S_{i}^{z}I_{j}^{z}\right)\right]}.$
(5)

Besides, the influence of additional nuclear species (such as B) is negated
within the dynamics due to the refocusing introduced by J-INSECT. This is, the
unavoidable presence of extra nuclear spin isotopes in the target sample would
not alter the system dynamics unless a $\pi$-pulse that synchronously flips
the B species is introduced on purpose (see later).

Then, the initial $\frac{\pi}{2}$-pulse along the y-axis takes the thermal
state of hydrogen nuclei to the x-axis (see Fig. 1) while ZZ-interactions that
govern the encoding stage –i.e. those in Eqs. (4) or (5)– encode J-couplings
in the dynamics of $\langle S_{i}^{x}\rangle$ and $\langle S_{i}^{y}\rangle$
of H nuclei. Importantly, owing to the initialization of the thermal
polarization along the x-axis, it can be demonstrated that following the
encoding stage, $\langle S_{i}^{x}\rangle$ exhibits significantly greater
amplitude than $\langle S_{i}^{y}\rangle$ (approximately $10^{5}$ times
larger), as detailed in SM suppmat . In this manner, J-INSECT is tuned to
target $\langle S_{i}^{x}\rangle$ in the detection stages.

At a generic $n$-th detection stage, owing to a rotation of the H nuclei along
the y-axis (see Fig. 1.), the magnetic field $B_{z}^{n}$ induced by the sample
over the NV ensemble is along the NV axis and reads meriles2010imaging ;
suppmat

$B_{z}^{n}(t)=B_{0}(n\tau)\sin(\Omega^{\text{H}}t),$ (6)

where $\Omega^{\text{H}}$ is the Rabi frequency of the $Y$ and $\bar{Y}$
pulses (see Fig. 1) while, importantly:

$B_{0}(n\tau)\propto\frac{2}{N^{\text{H}}}\sum{\langle
S_{i}^{x}\rangle_{t=n\tau}}.$ (7)

In conclusion, the subsequent magnetic field $B^{n}_{z}(t)$ has an amplitude
proportional to $\langle S_{i}^{x}\rangle$ of H nuclei, and oscillates with a
controllable frequency $\Omega^{\text{H}}$. This is, irrespective of the
involved Larmor frequencies, what enables to operate at high magnetic fields.
Note that, at 2 T hydrogen Larmors reach $\approx(2\pi)\times 85$ MHz, posing
a technically challenging task for their tracking. On the contrary,
$\Omega^{\text{H}}$ can be tuned, for instance to tens of kHz
herb2020broadband ; yudilevich2023coherent , such that the oscillations of
$B^{n}_{z}(t)$ are easily followed by the NV ensemble by applying, e.g., a
standard XY4 sequence. Recording the resultant phase acquisition of the NV
ensemble after every detection period, leads to a direct reconstruction of the
evolution of $\langle S_{i}^{x}\rangle$ of H nuclei, thus of the J-couplings
encoded in $\langle S_{i}^{x}\rangle$. Finally, it is worth remarking that the
NV centers are active only during the detection periods; this implies a high
spectral resolution constrained not by the sensor, but by the
$T_{2}(>T_{2}^{*})$ of H nuclei.

_Numerical results.–_ We consider a picoliter sized liquid sample on top of a
diamond hosting an NV ensemble at micrometer depth glenn2018high  and a static
magnetic field of $B_{\text{ext}}=2$ T aligned with the NV axis.

To exemplify J-INSECT we address a sample consisting on fluoromethanol
molecules ($\rm{CH_{2}FOH}$) that contains five spin-1/2 nuclei of three
different species: H, carbon and fluorine (13C and F from now on, note oxygen
has no spin) in a thermal state, see scheme of the molecule on inset of Fig.
2. Our protocol enables accessing heteronuclear J-couplings among H and the
13C nucleus, as well as with the F. Moreover, it gives access to the
homonuclear coupling between non-equivalent H (the two H nuclei in Fig. 2
labeled as H${}_{\text{a}}$ are magnetically equivalent, whereas
H${}_{\text{b}}$ is magnetically unequivalent with respect to the former). We
choose values falling within the typical range found in the literature
levitt2013spin , namely:
$J_{\rm{\text{H}_{\text{a}}\text{H}_{\text{b}}}}\equiv J=(2\pi)\times
8\>\rm{Hz}$, $J_{\rm{\text{H}_{\text{a}}C}}\equiv J_{1}=(2\pi)\times
130\>\rm{Hz}$, $J_{\rm{H_{b}C}}\equiv J_{2}=(2\pi)\times 6\>\rm{Hz}$, while
interactions with the F nuclei are $J_{\rm{H_{a}F}}\equiv
J^{F}_{1}=(2\pi)\times 80\>\rm{Hz}$ and $J_{\rm{H_{b}F}}\equiv
J^{F}_{2}=(2\pi)\times 4\>\rm{Hz}$, whereas $J_{CF}=(2\pi)\times 160$.
Furthermore, for a magnetic field $B_{\text{ext}}=2$ T, the chemical shifts
are $\delta_{\rm{H_{a}}}=(2\pi)\times 512\rm{Hz}$,
$\delta_{\rm{H_{b}}}=(2\pi)\times 236\rm{Hz}$, $\delta_{\rm{C}}=(2\pi)\times
85\rm{Hz}$ and $\delta_{\rm{F}}=(2\pi)\times 450\rm{Hz}$.

Our numerical simulations start with Hamiltonian (1) in an interaction picture
w.r.t the Larmor terms
$\omega_{H}\sum_{i=1}^{N_{\text{H}}}S_{i}^{z}+\omega_{\text{A}}\sum_{i}^{N_{\text{A}}}I_{i}^{z}$.
In this manner, the only approximation made is simplifying the heteronuclear
coupling term
$\sum_{i}^{N_{\text{H}}}\sum_{j}^{N_{\text{A}}}J_{i,j}^{\text{H}-\text{A}}\>\vec{S}_{i}\cdot\vec{I}_{j}$
in Eq. (1) to only ZZ interactions (i.e. to
$\sum_{i}^{N_{\text{H}}}\sum_{j}^{N_{\text{A}}}J_{i,j}^{\text{H}-\text{A}}\>S_{i}^{z}I_{j}^{z}$,
like in Eq. (2)). This can be done since
$|\omega_{\text{H}}-\omega_{i}|\sim(2\pi)\times 10\text{
MHz}>>J_{i,j}^{\text{H}-\text{A}}\sim(2\pi)\times 100\text{ Hz}$, thus, the
rotating wave approximation safely applies. The reader can find the full
simulated Hamiltonian in SM suppmat .

Figure 2: Spectra obtained from 1000 repetitions of the protocol with 600
_encoding_ \+ _detection_ stages. Due to the dephasing time $T_{2}=0.6$ s, the
FWHM is approximately 1 Hz. Within each figure, there is an illustrative
depiction of the molecule, highlighting the targeted J-couplings. Colored
nuclei signify the species upon which $\pi$ pulses are applied. Additionally,
each figure includes a boxed scheme detailing the RF pulse train applied on
each case. (a) Targeting only H and 13C, five peaks arise, which encode $J$,
$J_{1}$ and $J_{2}$. b) Targeting the H, the 13C and the F, ten peaks arise,
encoding: $J$, $J_{1}$, $J_{2}$, $J_{1}^{\rm{F}}$ and $J_{2}^{\rm{F}}$.

In order to maintain the Hamiltonian close to Eq. (2) the condition
$\tau>1/(|\delta_{\rm{\text{H}_{\text{a}}}}-\delta_{\rm{\text{H}_{\text{a}}}}|)$
must hold, such that the homonuclear term of H spins approaches to a ZZ
interaction. The previous condition imposes a lower bound on $\tau$, limiting
the amount of measurements that can be carried out for a given experimental
time. In our case, we set
$\tau=1.2/(|\delta_{\rm{\text{H}_{\text{a}}}}-\delta_{\rm{\text{H}_{\text{a}}}}|)=4.3$
ms (note later we explore the J-INSECT performance with a larger number of
measurements, i.e., departing from this condition). In addition, the amplitude
of the RF targeting H is $\Omega^{\text{H}}=(2\pi)\times 50$ kHz while, for
robustness purposes, $\pi$-pulses over H are implemented via a CORPSE
architecture Cummins2003tackling  which reduce the impact of detunings during
pulse execution. We observe no discernible impact on the system dynamics when
standard top-hat pulses are applied to 13C and F, thus we exclusively deploy
CORPSE pulses over Hs. This results in $t_{\pi}^{CORPSE}=43.3\ \mu$s,
$t_{\pi}^{\rm{C}}=39.7\ \mu$s and $t_{\pi}^{\rm{F}}=10.6\ \mu$s. Our
simulations firstly consider 600 measurements (this is, 600 encoding and
detection stages in a single experimental run) hence, the duration J-INSECT
per experimental run is
$T_{total}=600\>\left(\tau+t_{\pi}^{CORPSE}+2/\Omega^{\rm{\text{H}}}\right)\approx
2.66\>\rm{s}$, where the term $2/\Omega^{\rm{\text{H}}}$ corresponds to the
duration of two full RF Rabi oscillations during the readout stage.

Regarding error and imperfections, our numerical simulations comprise: (1)
Always on chemical shifts and J-couplings, (2) RF Cross-talk effects during RF
irradiation, (3) Errors on the amplitude of the RF drivings (these are
simulated through an Ornstein-Uhlenbeck process uhlenbeck1930on ; wang1945on ;
gillespie1996exact  with a 1% amplitude shift and a noise correlation time of
1 ms) and (4) A nuclear dephasing time $T_{2}=0.6\>\rm{s}$ glenn2018high ,
incorporated via an exponential decay in the field emitted by the sample.

We consider two cases: When J-INSECT targets H and 13C nuclei (Case 1) we get
the spectrum in Fig. 2.a. Note that in this case, F does not actively
participate, since we do not apply $\pi$-pulses over it, although it is
included in our numerical model. Case 1 exhibits five resonances at different
frequencies $f_{i}$ encoding heteronuclear J-couplings between H and 13C and
homonuclear J-coupling between H nuclei. The assignment of J-couplings based
on the obtained $f_{i}$ can be accomplished either by analytically solving the
evolution in Eq. (5) or by employing standard spin multiplicity techniques
suppmat ; multiplicity . Conversely, simultaneous targeting of the F nuclei
with an additional $\pi$-pulse (Case 2) results in the emergence of ten
resonances, as illustrated in Fig. 2.b. Note, in this second scenario the
heteronuclear coupling between the F and H nuclei is also encoded in
$B_{0}(n\tau)$. The resulting estimations for the J-couplings in Cases 1 and 2
are collected in Table 1.

| $\langle J\rangle$ | $\langle J_{1}\rangle$ | $\langle J_{2}\rangle$ | $\langle J_{1}^{\rm{F}}\rangle$ | $\langle J_{1}^{\rm{F}}\rangle$
---|---|---|---|---|---
Case 1 | $8.1\pm 0.5$ | $130.1\pm 0.8$ | $6.1\pm 0.9$ | - | -
Case 2 | $8.1\pm 0.5$ | $130.1\pm 0.9$ | $6\pm 1$ | $80.0\pm 0.9$ | $4\pm 1$
Table 1: J-coupling values obtained from the spectra shown in Fig.2.a (first
row) and Fig.2.b (second row). Values and their uncertainty are given in
Hertz.

In summary, J-INSECT effectively discloses J-couplings with large spectral
resolution resulting in a full-width-at-half-maximum (FWHM) of, at most, 1 Hz,
corresponding to the chosen value of $T_{2}=0.6$ s. Regarding sensitivities,
the use of elevated magnetic fields allows achieving a signal-to-noise ratio
(SNR) of 30 with approximately 18,000 repetitions of the protocol (estimation
carried out for Case 1 with photon statistics of NV experiments, assuming a 7%
of luminescence contrast. Additional details concerning Case 1 and the
development of Case 2, which yields similar results, are available in SM
suppmat ). This implies an experimental time of $\sim$13h. Hence, J-INSECT
facilitates a high-fidelity retrieval of J-couplings within a reasonable
signal-collecting time, with the potential to further reduce the experimental
duration by incorporating complementary methods, such as prior
hyperpolarization of the sample reported in Refs. bucher2020hyperpolarisation
; arunkumar2021micron , or advanced readout techniques such as repetitive
readout photoncount ; repetitivereadout .

An interesting final aside to consider is enhancing the rate of measurements,
as an strategy to extra reduce experimental time. To exemplify the procedure
we inspect the simpler Case 1. Enhancing the measurement rate implies reducing
$\tau$ in the encoding stages, hence departing from the condition
$\tau>1/(|\delta_{\rm{H_{a}}}-\delta_{\rm{H_{a}}}|)$. In this scenario the
spectrum is no longer given by the usual multiplicity rules suppmat , thus we
estimate the values of the J-couplings by performing Bayesian inference. In
particular, we choose to extend the number of measurements on each run (i.e.
the number of detection stages) from 600 to 6553 and reduce the overall
experimental time four times (from $\sim$ 13 h to $\sim$ 3h 30’). In this
case, the inference gives $\langle J\rangle=8.2\pm 0.5$, $\langle
J_{1}\rangle=130.1\pm 0.6$ and $\langle J_{2}\rangle=5.7\pm 0.4$ suppmat .

_Conclusion.–_ We introduced J-INSECT, a protocol that enables the detection
of homonuclear and heteronuclear J-couplings with NV ensembles at high
magnetic fields. J-INSECT incorporates heterodyne techniques and exerts noise
filtering, enhancing the sample’s coherence thus enabling highly precise
estimations of J-couplings limited by $T_{2}>T_{2}^{*}$. We explore J-INSECT’s
performance in two scenarios, considering both reduced and enhanced data
harvesting. We assess the influence of different noise sources, including the
persistent presence of always-on J-couplings, showcasing the method’s
feasibility in experimental conditions.

_Acknowledgements.–_ P. A. B. and A. B. U. acknowledge the financial support
of the IKUR STRATEGY (IKUR-IKA-23/22) and (IKUR-IKA-23/04), respectively. C.
M. J. acknowledges the predoctoral MICINN Grant No. PRE2019-088519. J. C.
acknowledges the Ramón y Cajal (RYC2018-025197-I) research fellowship, the
financial support from Spanish Government via the Nanoscale NMR and complex
systems (PID2021-126694NB-C21) project, the ELKARTEK project Dispositivos en
Tecnologías Cuánticas (KK-2022/00062), and the Basque Government grant
IT1470-22.

## References

  * (1) C. P. Slichter, Principles of Magnetic Resonance (Springer-Verlag, Berlin Heidelberg, 1978).
  * (2) M. H. Levitt, Spin Dynamics: basics of Nuclear Magnetic Resonance, 2nd ed. (Wiley, Chichester, 2008).
  * (3) A. Tampieri, M. Szabó, F. Medina and H. Gulyás Physical Sciences Review 6, 20190086 (2022).
  * (4) K.V.R. Chary, G. Govil,NMR in Biological Systems. From Molecules to Humans, Vol. 6 (Springer Science and Business Media, Dordrecht, 2008).
  * (5) P.A. Bottomley, Computerized radiology 8, 57 (1984)
  * (6) A. Abragam, The Principles of Nuclear Magnetism, (Oxford University Press, London, 1961).
  * (7) R. D. Allert, K. D. Breigel, and D. B. Bucher, Chemical Communications 58, 8165 (2022).
  * (8) N. R. Neuling, R.D. Allert and D. B. Bucher, Current Opinion in Biotechnology 83, 102975 (2023).
  * (9) A. O. Sushkov, N. Chisholm, I. Lovchinsky, M. Kubo, P. K. Lo, S. D. Bennett, D. Hunger, A. Akimov, R. L. Walsworth, H. Park and M. D. Lukin, Nano Letters 14, 6443 (2014)
  * (10) F. Shi, P. Wang, H. Sun, J. Wang, X. Rong, M. Chen, C. Ju, F. Reinhard, H. Chen, J. Wrachtrup, J. Wang and J. Du, Science 347, 1135 (2015)
  * (11) I. Lovchinsky, A. O. Sushkov, E. Urbach, N. P. de Leon, S. Choi, K. De Greve, R. Evans, R. Gertner, E. Bersin, C. Müller, L. McGuiness, F. Jelezko, R. L. Walsworth, H. Park and M. D. Lukin, Science 351, 836 (2016)
  * (12) D. R. Glenn, D. B. Bucher, J. Lee, M. D. Lukin, H. Park, and R. L. Walsworth, Nature, 555, 351 (2018).
  * (13) D. B. Bucher, D. R. Glenn, H. Park, M. D. Lukin and R. L. Walsworth, Physical Review X 10, 021053 (2020).
  * (14) N. Arunkumar, D. B. Bucher, M. J. Turner, P. TomHon, D. Glenn, S. Lehmkuhl, M. D. Lukin, H. Park, M. S. Rosen, T. Thesis and R. L. Walsworth, PRX Quantum 2, 010305 (2021).
  * (15) R. Liu, X. Li and K. S. Lam Current Opinion in Chemical Biology 38, 117 (2017).
  * (16) B. Liang and L K. Nature Structural & Molecular Biology 23, 468 (2016).
  * (17) S. Schmitt, T. Gefen, F. M. Stürner, T. Unden, G. Wolff, C. Müller, J. Scheuer, B. Naydenov, M. Markham, S. Pezzagna, J. Meijer, I. Schwarz, M. Plenio, A. Retzker, L. P. McGuinness and Fedor Jelezko, Science 356, 832 (2017)
  * (18) C. Munuera-Javaloy, A. Tobalina and J. Casanova, Physical Review Letters 130, 133603 (2023)
  * (19) J. Meinel, M. S. Kwon, D. Dasari, H. Sumiya, S. Onoda, J. Isoya, V. Vorobyov and J. Wrachtrup, Communication Physics 6, 302 (2023)
  * (20) J. M. Boss, K. S. Kujia, J. Zopes and C. L. Degen, Science 356, 837 (2017)
  * (21) The full form of the J-coupling is a tensor, although in an isotropic liquid –this is, in common liquids– this tensor is averaged out by rapid molecular tumbling, leading to a scalar interaction. Besides, for the same reason, the direct dipole-dipole coupling can be neglected.
  * (22) D. Daly, E. Huckestein, J. W. Blanchard and R. L. Walsworth, arXiv preprint arXiv:2310.08499 (2023)
  * (23) Supplemental Material
  * (24) C.A. Meriles, L. Jiang, J. S. Hodges, J. Maze, M. D. Lukin, P. Cappellaro, The Journal of Chemical Physics 133, 124105 (2010).
  * (25) K. Herb, J. Zopes, K. S. Cujia and C. L. Degen , Review of Scientific Instruments 91, 113106 (2020)
  * (26) D. Yudilevich, A. Salhov, I. Schaefer, K. Herb, A. Retzker and A. Finkler, New Journal of Physics 25, 113042 (2023)
  * (27) H. K. Cummins, G. Llewellyn and J.A. Jones, Physical Review A 67, 042308 (2003)
  * (28) G. E. Uhlenbeck and L.S. Ornstein, Physical Review 36, 823 (1930)
  * (29) M. C. Wang and G. E. Uhlenbeck, Reviews of Modern Physics 17, 323 (1945)
  * (30) D. T. Gillespie, Physical Review E 54, 2084 (1996)
  * (31) L. Jiang, J. S. Hidges, J. R. Maze, P. Maurer, J. M. Taylor, D. G. Cory, P. R. Hemmer, R.L. Walsworth, A. Yacoby, A. S. Zibrov, M.D. Lukin, Science 326, 267 (2009).
  * (32) E. Breitmaier, Stucture Elucidation by NMR in Organic Chemistry, 3rd ed. (John Wiley and Sons, 2003).
  * (33) N. Arunkumar, K. S. Olsson, J. T. Oon, C. A. Hart, D. B. Bucher, D. R. Glenn, M. D. Lukin, H. Park, D. Ham, and R. L. Walsworth, Physical Review Letters 131, 100801 (2023).

Supplemental Material:
J-coupling NMR Spectroscopy with Nitrogen Vacancy Centers at High Fields

##  I Spin dynamics of the sample

In this section, we provide further details regarding the nuclear spin
dynamics in the sample.

###  I.1 Amplitudes of the $\langle\sigma_{i}^{x}\rangle$ and
$\langle\sigma_{i}^{y}\rangle$ components

We consider an evolution governed by the next propagator (note this is similar
to those in Eqs. (4,5) in the main text)

$U_{t}=\Pi_{i,j}e^{-i\frac{\phi_{i,j}}{2}\sigma_{i}^{z}\sigma_{j}^{z}},$ (S1)

over an initial state such that

$(\rho_{1}^{x}...\rho_{k}^{x})(\rho_{k+1}^{z}...\rho_{N}^{z}),$ (S2)

where $\rho_{j}^{x,z}$ are thermal states polarized along $x$ and $z$
directions. This is, $\rho_{j}^{z}=\frac{1}{2}I+\frac{B_{\rm n}}{4}\sigma_{z}$
while
$\rho_{j}^{x}=e^{i\frac{\pi}{4}\sigma_{j}^{y}}\rho_{j}^{z}e^{-i\frac{\pi}{4}\sigma_{j}^{y}}=\frac{1}{2}I-\frac{B_{\rm
H}}{4}\sigma_{x}$, with $B_{H}=\frac{\hbar\gamma_{h}B_{\rm{ext}}}{K_{b}T}$ and
$B_{n}=\frac{\hbar\gamma_{n}B_{\rm{ext}}}{K_{b}T}$ being small quantities at
room temperature (on the order of $10^{-5}$ for a magnetic field of 2 T). In
the previous expressions, $B_{\rm ext}=2\>\rm{T}$ is the external magnetic
field, $\hbar=1.054\cdot 10^{-34}\>\rm{J\cdot s}$ is the Planck constant
divided by $2\pi$, $\gamma_{H}=(2\pi)\times 42.57\>\rm{MHz/T}$ is the
gyromagnetic ratio of the H and $\gamma_{n}$ the corresponding one of other
nuclei, $K_{B}=1.38\cdot 10^{-23}\>\rm{J/K}$ the Boltzmann constant, and the
temperature $T=300$ K. Here we differentiate among thermal states
corresponding to the hydrogen (H) nuclei, i.e. the nuclear spin which is
rotated to the X axis at the start of the protocol (denoted by
$\rho_{j}^{x}$), and the other generic nuclear spin, which is represented with
$\rho_{j}^{z}$.

On the one hand, one can demonstrate that

$\langle\sigma_{1}^{x}\rangle={\rm
Tr}[U_{t}(\rho_{1}^{x}...\rho_{k}^{x})(\rho_{1}^{z}...\rho_{k}^{z})U^{\dagger}_{t}\sigma_{1}^{x}]={\rm
Tr_{\neq 1}}[\langle
e|\rho_{1}^{x}|g\rangle...\rho_{k}^{x}\rho_{k+1}^{z}...\rho_{N}^{z}\Pi_{j>1}e^{-i\phi_{1,j}\sigma_{j}^{z}}+\langle
g|\rho_{1}^{x}|e\rangle...\rho_{k}^{x}\rho_{k+1}^{z}...\rho_{N}^{z}\Pi_{j>1}e^{i\phi_{1,j}\sigma_{j}^{z}}],$
(S3)

where $\rm{Tr}_{\neq 1}[...]$ represents the trace over all spins excepting
the first one. Using $\langle e|\rho_{1}^{x}|g\rangle=\langle
g|\rho_{1}^{x}|e\rangle=-\frac{B_{\rm H}}{4}$ one reaches

$\displaystyle\langle\sigma_{1}^{x}\rangle=-\frac{B_{\rm
H}}{4}\Pi_{j=1}^{k}\cos{\phi_{1,j}}\
\Pi_{j=k+1}^{N}\bigg{(}\cos{\phi_{1,j}}-\frac{iB_{\rm
n}}{2}\sin{\phi_{1,j}}\bigg{)}-\frac{B_{\rm
H}}{4}\Pi_{j=1}^{k}\cos{\phi_{1,j}}\
\Pi_{j=k+1}^{N}\bigg{(}\cos{\phi_{1,j}}+\frac{iB_{\rm
n}}{2}\sin{\phi_{1,j}}\bigg{)}.$ (S4)

On the other hand, one can demonstrate that $\langle\sigma^{y}_{1}\rangle$
leads to

$\displaystyle\langle\sigma_{1}^{y}\rangle=-\frac{iB_{\rm
H}}{4}\Pi_{j=1}^{k}\cos{\phi_{1,j}}\
\Pi_{j=k+1}^{N}\bigg{(}\cos{\phi_{1,j}}-\frac{iB_{\rm
n}}{2}\sin{\phi_{1,j}}\bigg{)}+\frac{iB_{\rm
H}}{4}\Pi_{j=1}^{k}\cos{\phi_{1,j}}\
\Pi_{j=k+1}^{N}\bigg{(}\cos{\phi_{1,j}}+\frac{iB_{\rm
n}}{2}\sin{\phi_{1,j}}\bigg{)}.$ (S5)

Note the above expressions hold for any H, thus for any quantity
$\langle\sigma_{i}^{x}\rangle$ and $\langle\sigma_{i}^{y}\rangle$ that
involves a H spin in the target sample.

Notice how the leading order for each expectation value differs by five order
of magnitude:
$\langle\sigma_{1}^{x}\rangle\propto\mathcal{O}\left(B_{H}\right)$ and
$\langle\sigma_{1}^{y}\rangle\propto\mathcal{O}\left(B_{H}B_{n}\right)$. Thus,
$\langle\sigma_{i}^{x}\rangle$ is $\sim 10^{5}$ times greater than
$\langle\sigma_{i}^{y}\rangle$.

###  I.2 Fluoromethanol resonances

We compute the resonance positions for the specific case of the fluoromethanol
molecule. For the sake of simplicity we consider a polarised H, and the other
nuclei in the maximally mixed state, this is
$\rho_{i}^{x}=\frac{1}{2}(I+\sigma_{i}^{x})$ for H and
$\rho_{i}^{m}=\frac{1}{2}I$ for other nuclei. Note, these considerations only
alter the amplitude of the signal, not the resonance position. However, the
results in the main text consider thermal states for all nuclei in the
fluoromethanol molecule. Thus, following Eq. (S3), we have

$\displaystyle\langle
S_{i}^{x}\rangle=\frac{\langle\sigma_{i}^{x}\rangle}{2}=\frac{1}{2}\Pi_{j=1}^{N}\cos{\phi_{i,j}}.$
(S6)

where $\phi_{1,i}$ corresponds to the phase accumulated due to the interaction
between the i-th and j-th spins during a time t, that is
$\phi_{i,j}(\text{t})=\frac{J_{i,j}}{2}\text{t}$.

When H nuclei and 13C are targeted in the fluoromethanol molecule
(corresponding to Case 1 in the main text), the resulting signal is

$\displaystyle\langle\sum_{i=1}^{3}{S^{x}_{i}}\rangle(\text{t})=\frac{1}{6}\cos{\left(\frac{J}{2}\text{t}\right)}\left[2\cos{\left(\frac{J_{1}}{2}\text{t}\right)}+\cos{\left(\frac{J}{2}\text{t}\right)}\cos{\left(\frac{J_{2}}{2}\text{t}\right)}\right],$
(S7)

which leads to the resonance peaks in Fig.2.a of the main text, $f_{1}=\pm
J_{2}/2$, $f_{2}=\pm(J_{2}-2J)/2$, $f_{3}=\pm(J_{2}+2J)/2$,
$f_{4}=\pm(J_{1}-J)/2$, and $f_{5}=\pm(J_{1}+J)$.

In the case where the fluorine (F) is also targeted with $\pi$-pulses
(corresponding to Case 2 in the main text), one finds

$\displaystyle\langle\sum_{i=1}^{3}{S^{x}_{i}}\rangle(\text{t})=\frac{1}{6}\cos{\left(\frac{J}{2}\text{t}\right)}\left[2\cos{\left(\frac{J_{1}}{2}\text{t}\right)}\cos{\left(\frac{J_{1}^{\rm{F}}}{2}\text{t}\right)}+\cos{\left(\frac{J}{2}\text{t}\right)}\cos{\left(\frac{J_{2}}{2}\text{t}\right)}\cos{\left(\frac{J_{2}^{\rm{F}}}{2}\text{t}\right)}\right],$
(S8)

with the corresponding resonances depicted in Fig.2.b of the main text. These
are: $f_{1}=\pm\left(J_{2}-J_{2}^{\rm{F}}\right)/2$,
$f_{2}=\pm\left(+2J-J_{2}-J_{2}^{\rm{F}}\right)/2$,
$f_{3}=\pm\left(J_{2}+J_{2}^{\rm{F}}\right)/2$,
$f_{4}=\pm\left(2J+J_{2}^{\rm{F}}-J_{2}\right)/2$,
$f_{5}=\pm\left(2J+J_{2}-J_{2}^{\rm{F}}\right)/2$,
$f_{6}=\pm\left(2J+J_{1}+J_{1}^{\rm{F}}\right)/2$,
$f_{7}=\pm\left(J_{1}-J_{1}^{\rm{F}}-J\right)/2$,
$f_{8}=\pm\left(J+J_{1}-J_{1}^{\rm{F}}\right)/2$,
$f_{9}=\pm\left(J_{1}+J_{1}^{\rm{F}}-J\right)/2$ and
$f_{10}=\pm\left(J+J_{1}+J_{1}^{\rm{F}}\right)/2$.

Alternatively, the resonant frequencies can be obtained via standard NMR
multiplicity rules smultiplicity . This is, when a nuclei A is coupled with
$n_{x}$ spin-$\frac{1}{2}$ nuclei of an isotope X, the A signal splits into
$n_{x}+1$ lines. If A is also coupled to an element different to X (lets say
Y) this induces additional $n_{y}+1$ splittings in each previous peak.

In our case, since the broadening is mitigated, all J-coupling peaks tend to
be centered at 0 Hz. Thus, when H and 13C are targeted, the signal of $H_{a}$
(recall the $H_{a}$ and $H_{b}$ are magnetically inequivalent hydrogen groups
in the molecule, refer to Fig. 2 of the main text) is split into $2\times 2=4$
peaks, due to couplings with $H_{b}$ and the 13C, while the signal of $H_{b}$
is split into $3\times 2=6$ peaks, 3 due to the pair of equivalent $H_{a}$,
and 2 due to the 13C. Then, we are left with $6+4=10$ peaks around 0 Hz. If J
coupling interactions with F are also excited, the signal of $H_{a}$ is split
into $2\times 2\times 2=8$ (2 per every other element), while the signal of
$H_{b}$ splits into $3\times 2\times 2=12$ peaks, leading to a total of 20
peaks, see, e.g. sabragam1961principles  for more information.

##  II From the sample to the NV centers

The oscillation of hydrogens in the detection stage of J-INSECT impacts on the
NV ensemble (which is considered at a micrometer depth from diamond surface)
as a classical magnetic field. In this respect, see e.g. sglenn2018high ;
smeriles2010imaging ). In the high field regime, the magnetic signal induced
by the sample in the orthogonal directions (i.e. x and y) with respect to the
external magnetic field (oriented with the NV axis, i.e. $z$) can be neglected
as it rapidly oscillate at a speed $\gamma_{N}B_{z}$ with $B_{z}$ being
several Tesla. The resultant field emanated from the sample during a n-th
detection period ($n\in\mathbb{N}$) can be well approximated by sglenn2018high
; smeriles2010imaging

$\displaystyle B_{z}(t)=B_{0}^{n}\sin(\Omega^{H}t)$ (S9)

where $\Omega^{H}=\gamma_{H}B_{RF}$ is the amplitude of the RF driving. On the
other hand, the amplitude of the magnetic field is

$B_{0}^{n}=\frac{(2\pi)^{2}(\hbar\>\gamma_{H})^{2}\mu_{0}\>\rho_{H}\>B_{z}F_{3}}{16\pi
K_{B}T}M_{x}$ (S10)

where $M_{x}$ is the x component of the normalized magnetisation, defined as
$\vec{M}=\frac{2}{N^{H}}\langle\sum_{i}{\vec{S}_{i}}\rangle(n\tau)$. Note that
the quantity $\langle\sum_{i}{\vec{S}_{i}}\rangle(n\tau)$ is computed via
numerically evolving the next molecular Hamiltonian:

$\displaystyle H/\hbar=$
$\displaystyle\delta^{H_{A}}(S_{1}^{z}+S_{2}^{z})+\delta^{H_{B}}S_{3}^{z}+\delta^{C}I^{z}+\delta^{F}P^{z}+J_{H_{a}H_{a}}\vec{S_{1}}\cdot\vec{S_{2}}+J_{H_{a}H_{b}}(\vec{S_{1}}\cdot\vec{S_{3}}+\vec{S_{2}}\cdot\vec{S_{3}})$
$\displaystyle+J_{H_{a}C}(S_{1}^{z}+S_{2}^{z})I^{z}+J_{H_{b}C}S_{3}^{z}I^{z}+J_{H_{a}F}(S_{1}^{z}+S_{2}^{z})P^{z}+J_{H_{b}F}S_{3}^{z}P^{z}+J_{CF}I^{z}P^{z}$
(S11)

In addition, $\hbar=1.054\cdot 10^{-34}\>\rm{J\cdot s}$ is the Planck constant
divided by $2\pi$, $\gamma_{H}=(2\pi)\times 42.57\>\rm{MHz/T}$ is the
gyromagnetic ratio of the H, $\mu_{0}=4\pi\cdot 10^{-7}\>\rm{H/m}$ the
magnetic permeability of free space, $\rho_{H}\approx 6.6\times
10^{28}\>m^{-3}$ the H density of the sample, $B_{z}=2\>\rm{T}$ the external
magnetic field, $K_{B}=1.38\cdot 10^{-23}\>\rm{J/K}$ the Boltzmann constant,
$T=300K$ the temperature and $F_{3}=4.1$ characterizes the geometry of the
sample sglenn2018high .

The Hamiltonian of a generic NV in the ensemble under the influence of
$B_{z}(t)$ is

$\displaystyle H_{NV}=\frac{\sigma_{z}\>\gamma_{e}}{2}B_{z}(t)+H_{c},$ (S12)

where $H_{c}$ accounts for the MW driving on the NV ensemble. In particular,
applying a XY4 control sequence during the detection stage over the NV-
ensemble with $t_{RF}=2\frac{1}{\Omega^{H}}$ and the spacing between the MW
pulses $T_{MW}=t_{RF}/4$, yields:

$\displaystyle\langle\sigma_{y}\rangle_{NV}\approx\frac{2\gamma_{e}\>t_{RF}}{\pi}B_{0}^{n},$
(S13)

provided that $\frac{2\>\gamma_{e}t_{RF}}{\pi}\>B_{0}^{n}<<1$.

##  III READOUT CONSIDERATIONS

In order to evaluate the impact of the measurement rate of the protocol under
the high external magnetic field paradigm, we make an estimation of the
required signal averaging time leading to a SNR of 30.

In sphotoncount  they record the mean number of photon-count per measurement
when a single NV is in $\ket{0}$, that is $n_{0}=0.016$, and a photon count
difference between the ground and excited spin state of $n_{0}-n_{1}$=0.005,
yielding a $30\%$ of luminescence contrast between both states. In this case,
since we are working with NV ensembles, we considered a fluorescence contrast
of $7\%$, which was achieved by lowering the photon emission $n_{0}$.
Employing the latter photon statistics, and taking as a basis the signal with
no readout error (in particular, the inverse Fourier transform of the spectra
in Fig. 2.a, note that this signal carries the rest of errors described in the
main text), we have reproduced an experimental outcome.

In particular, we assume a configuration similar to that in Ref.
sglenn2018high . This is an NV center density of $0.8\times
10^{23}\text{m}^{3}$, a beam diameter of $20\ \mu$m, and a $10\ \mu$m thick
layer of NV centers sensitive to the thermal spin signal. This leads to a
sensor containing $\approx 2.5\times 10^{8}$ active NV centers. The protocol
is repeated $18000$ times. This is, we average the emitted signal across
$18000$ repetitions of the protocol, which leads to a duration of $18000\times
2.66\text{s}\approx\text{13\text{ h }}20^{\prime}$.

Specifically, in Case 1, the noise in the resulting averaged signal has a
standard deviation of approximately $6\times 10^{-5}$ while the signal’s
amplitude is $\langle\sigma_{NV}^{z}\rangle\sim 1\times 10^{-3}$, see inset of
Fig. S1.a. This leads to a SNR of about 30 in the spectrum by taking the the
height of the smallest peak. For Case 2, the noise standard deviation is
similar. However, the SNR is approximately 15. This is because even though
both signals (from Case 1 and Case 2) have similar amplitudes, the spectrum of
Case 2 exhibits more peaks which results in peaks of smaller amplitude as it
can be seen in Fig. S1 for clarification.

In conclusion, the thermal polarization induced by the high external magnetic
field is enough to balance out the relatively low measurement rate, yielding a
detectable signal (thus of clean J-coupling spectra) under reasonable
experimental time in both cases.

Figure S1: Simulating an experiment. Figure a) corresponds to Case 1, while b)
holds for Case 2. Insets: in blue are the signals without considering
luminescence issues, while red squares are the outcome that results from
averaging $18000$ repetitions of the respective protocol with $2.5\times
10^{8}$ NV centers assuming a contrast of 7%. Main figures: In blue are the
spectra that correspond to the signals without luminescence issues, while the
red squares are the FT of averaged signal in the inset (thus, carrying readout
noise). Both spectra are normalized by a common factor, ensuring that the
maximum peak in Case 1 has amplitude 1. In Case 1, referencing the amplitude
of the smallest peak, the corresponding SNR is just above 30. In Case 2 under
the same reference criterion, it is slightly below 15.

##  IV Beyond $\tau>1.2/(|\delta_{\rm{H_{a}}}-\delta_{\rm{H_{a}}}|)$

To increase the number of measurements on each experimental run, we reduce
$\tau$ in encoding stages from $\tau=4.3\text{ ms}$ to $\tau=75\>\mu s$ and
perform 6553 measurements instead of 600. We consider imperfect pulses (with
1% error on the amplitude driving) as well as always-on J-couplings and
readout errors (see previous section). This increase in the number of
measurements leads to a reduction on the total experimental repetitions to be
performed. In this case the protocol includes an extra ($\pi$)-pulse during
the encoding period to further reduce the impact of imperfect drivings (see
the scheme of the pulse sequence in Fig. S2).

On the other hand, the reduction on $\tau$ leads to an intricate spectrum
without clear resonance peaks, see Fig. S2.a). Thus, to infer the values of
J-couplings we perform Bayesian analysis. As an input to the Bayesian
inference we employ the previously described noisy signal with imperfect
pulses and readout errors. Regarding the parameter set, as the prior knowledge
we assume a uniform distribution. To sum up, via Bayes inference, simplicity
of the spectra can be traded for shorter experimental times.

As a figure of merit, we simulate the resulting signal performing four times
less averages (thus a four-times shorter experiment). This results into a
readout error twice greater. The posterior distributions are depicted in Fig
S2.b, from which the following values and uncertainties arise (respectively):
$\langle J\rangle=8.2\pm 0.5$ Hz, $\langle J_{1}\rangle=130.1\pm 0.6$ Hz and
$\langle J_{2}\rangle=5.7\pm 0.4$ Hz.

Figure S2:  a) Top, scheme of the RF sequence employed for the high
measurement-rate protocol (with $(\pi)_{x}$ and $(\pi)_{-x}$ pulses altogether
with CORPSE structures over Hs during encoding for enhanced robustness) and
the corresponding spectrum (bottom). This shows the result obtained with 6553
measurements (solid blue) while, for comparison, we include the spectrum
obtained with 600 measurements (dashed red), which corresponds to the same
spectrum shown in Fig. 2.a of the main text. b) Posterior distributions for
each scanned parameter via Bayesian inference. I: Posterior distribution
corresponding to J, II: to $\rm J_{1}$ and III: to $\rm J_{2}$. The exact
value of each parameter is marked by a vertical dotted line in I, II, and III.
The uncertainties are the standard deviation $\sigma_{\theta}$, given by
$\sigma_{\theta}^{2}=\langle(\theta-\langle\theta\rangle)^{2}\rangle=\int(\theta-\langle\theta\rangle)^{2}\>P(\theta|\rm{D})\>d\theta$,
where $\theta$ stands for J, $\rm J_{1}$ and $\rm J_{2}$.

## References

  * (1) E. Breitmaier, Stucture Elucidation by NMR in Organic Chemistry, 3rd ed. (John Wiley and Sons, 2003).
  * (2) A. Abragam, The Principles of Nuclear Magnetism, (Oxford University Press, London, 1961).
  * (3) D. R. Glenn, D. B. Bucher, J. Lee, M. D. Lukin, H. Park, and R. L. Walsworth, Nature, 555, 351 (2018).
  * (4) C.A. Meriles, L. Jiang, J. S. Hodges, J. Maze, M. D. Lukin, P. Cappellaro, The Journal of Chemical Physics 133, 124105 (2010).
  * (5) L. Jiang, J. S. Hidges, J. R. Maze, P. Maurer, J. M. Taylor, D. G. Cory, P. R. Hemmer, R.L. Walsworth, A. Yacoby, A. S. Zibrov, M.D. Lukin, Science 326, 267 (2009).