# CW-ERM: Improving Autonomous Driving Planning with
Closed-loop Weighted Empirical Risk Minimization

Eesha Kumar1,∗, Yiming Zhang2, Stefano Pini1, Simon Stent1,
Ana Ferreira2, Sergey Zagoruyko1, Christian S. Perone1,∗ 1Author is with Woven
Planet United Kingdom Limited, 114-116 Curtain Road, London, United Kingdom,
EC2A 3AH<EMAIL_ADDRESS>is with Woven Planet
North America, Inc., 900 Arastradero Rd, Palo Alto, CA, USA 94304.
<EMAIL_ADDRESS>contribution.

###### Abstract

The imitation learning of self-driving vehicle policies through behavioral
cloning is often carried out in an open-loop fashion, ignoring the effect of
actions to future states. Training such policies purely with Empirical Risk
Minimization (ERM) can be detrimental to real-world performance, as it biases
policy networks towards matching only open-loop behavior, showing poor results
when evaluated in closed-loop. In this work, we develop an efficient and
simple-to-implement principle called Closed-loop Weighted Empirical Risk
Minimization (CW-ERM), in which a closed-loop evaluation procedure is first
used to identify training data samples that are important for practical
driving performance and then we these samples to help debias the policy
network. We evaluate CW-ERM in a challenging urban driving dataset and show
that this procedure yields a significant reduction in collisions as well as
other non-differentiable closed-loop metrics.

##  I Introduction

Learning effective planning policies for self-driving vehicles (SDVs) from
data such as human demonstrations remains one of the major challenges in
robotics and machine learning. Since early works such as ALVINN [1], Imitation
Learning has seen major recent developments using modern Deep Neural Networks
(DNNs) [2, 3, 4, 5, 6, 7]. Imitation Learning (IL), and especially Behavioral
Cloning (BC), however, still face fundamental challenges [8], including causal
confusion [9] (later identified as a feedback-driven covariate shift[10]) and
dataset biases [8], to name a few.

There is one particular limitation of IL policies trained with BC that is,
however, often overlooked: the mismatch between training and inference-time
execution of the policy actions. Most of the time, BC policies are trained in
an open-loop fashion, predicting the next action given the immediate previous
action and optionally conditioned on recent past actions [2, 3, 4, 5, 7].
These policies, however, when executed in real-world, impact the future
states. Small prediction errors can drive covariate shift and make the network
predict in an out-of-distribution regime.

In this work, we address the mismatch between training and inference through
the development of a simple training principle. Using a closed-loop simulator,
we first identify and then reweight samples that are important for the closed-
loop performance of the planner. We call this approach CW-ERM (Closed-loop
Weighted Empirical Risk Minimization), since we use Weighted ERM [11] to
correct the training distribution in favour of closed-loop performance. We
extensively evaluate this principle on real-world urban driving data and show
that it can achieve significant improvements on planner metrics that matter
for real-world performance (e.g. collisions).

Our contributions are therefore the following:

  * •

We motivate and propose Closed-loop Weighted Empirical Risk Minimization (CW-
ERM), a technique that leverages closed-loop evaluation metrics acquired from
policy rollouts in a simulator to debias the policy network and reduce the
distributional differences between training (open-loop) and inference time
(closed-loop);

  * •

we evaluate CW-ERM experimentally on a challenging urban driving dataset in a
closed-loop fashion to show that our method, although simple to implement,
yields significant improvements in closed-loop performance without requiring
complex and computationally expensive closed-loop training methods;

  * •

we also show an important connection of our method to a family of methods that
addresses covariate shift through density ratio estimation.

In Section II, we detail the proposed CW-ERM and in Section IV we show the CW-
ERM experiments and compare them against ERM.

##  II Methodology

Figure 1: High-level overview of our proposed Closed-loop Weighted Empirical
Risk Minimization (CW-ERM) method. In steps (1-2) we train an identification
policy $\hat{\pi}_{\text{ERM}}$ using traditional ERM [12] on a set of
training data samples or driving "scenes". In step (3), we perform closed-loop
simulation of the policy $\hat{\pi}_{\text{ERM}}$ and collect metrics to
construct the error set in step (4). With the error set in hand, we upsample
scenes in the training set as shown in step (5). We train the final policy
$\hat{\pi}_{\text{CW-ERM}}$ using CW-ERM as shown in step (6) with the
upsampled $\mathcal{D}_{\text{up}}$ set.

###  II-A Problem Setup

The traditional formulation of supervised learning for imitation learning,
also called behavioral cloning (BC), can be formulated as finding the policy
$\hat{\pi}_{BC}$:

$\hat{\pi}_{BC}=\operatorname*{argmin}_{\pi\in\Pi}\mathbb{E}_{s\sim
d_{\pi^{*}},a\sim\pi^{*}(s)}[\ell(s,a,\pi)]$ (1)

where the state $s$ is sampled from the expert state distribution
$d_{\pi^{*}}$ induced when following the expert policy $\pi^{*}$. Actions $a$
are sampled from the expert policy $\pi^{*}(s)$. The loss $\ell$ is also known
as the surrogate loss that will find the policy $\hat{\pi}_{BC}$ that best
mimics the unknown expert policy $\pi^{*}(s)$. In practice, we only observe a
finite set of state-action pairs ${(s_{i},a^{*}_{i})}_{i=1}^{m}$, so the
optimization is only approximate and we then follow the Empirical Risk
Minimization (ERM) principle to find the policy $\pi$ from the policy class
$\Pi$.

If we let $\mathbb{E}_{s\sim
d_{\pi^{*}},a\sim\pi^{*}(s)}[\ell(s,a,\pi)]=\epsilon$, then it follows that
$J(\pi)\leq J(\pi^{*})+T^{2}\epsilon$ as shown by the proof in [13], where $J$
is the total cost and $T$ is the task horizon. As we can see, the total cost
can grow quadratically in $T$.

When the policy $\hat{\pi}_{BC}$ is deployed in the real-world, it will
eventually make mistakes and then induce a state distribution
$d_{\hat{\pi}_{BC}}$ different than the one it was trained on (
$d_{\pi^{*}}$). During closed-loop evaluation of driving policies, non-
imitative metrics such as collisions and comfort are also evaluated. However,
they are often ignored in the surrogate loss or only implicitly learned by
imitating the expert due to the difficulty of overcoming differentiability
requirements, as smooth approximations of these metrics are still different
than the non-differentiable counterparts often used. These policies can often
show good results in open-loop training, but perform poorly in closed-loop
evaluation or when deployed in a real SDV due to the differences between
$d_{\hat{\pi}_{BC}}$ and $d_{\pi^{*}}$, where the estimator is no longer
consistent.

###  II-B Closed-loop Weighted Empirical Risk Minimization

In our method, called “Closed-loop Weighted Empirical Risk Minimization” (CW-
ERM), we seek to debias a policy network from the open-loop performance
towards closed-loop performance, making the model rely on features that are
robust to closed-loop evaluation. Our method consists of three stages: the
training of an identification policy, the use of that policy in closed-loop
simulation to identify samples, and the training of a final policy network on
a reweighted data distribution. More explicitly:
Stage 1 (identification policy): train a traditional BC policy network in
open-loop using ERM, to yield $\hat{\pi}_{\text{ERM}}$.

Stage 2 (closed-loop simulation): perform rollouts of the
$\hat{\pi}_{\text{ERM}}$ policy in a closed-loop simulator, collect closed-
loop metrics and then identify the error set below:

$\displaystyle
E_{\hat{\pi}_{\text{ERM}}}=\\{(s_{i},a_{i})~{}\text{s.t.}~{}{C(s_{i},a_{i})}>0\\},$
(2)

where $\text{s}_{i}$ is a training data sample, or “scene” with a fixed number
of timesteps from the training set, $\text{a}_{i}$ is the action performed
during the roll-out and $C(\cdot)$ is a cost such as the number of collisions
found during closed-loop rollouts.

Stage 3 (final policy): train a new policy using weighted ERM where the scenes
belonging to the error set $E_{\hat{\pi}_{\text{ERM}}}$ are upweighted by a
factor $w(\cdot)$, yielding the policy $\hat{\pi}_{\text{CW-ERM}}$:

$\operatorname*{argmin}_{\pi\in\Pi}\mathbb{E}_{s\sim
d_{\pi^{*}},a\sim\pi^{*}(s)}[w(E_{\hat{\pi}_{\text{ERM}}},s)\ell(s,a,\pi)]$
(3)

As we can see, the CW-ERM policy in Equation 3 is very similar to the original
BC policy trained with ERM in Equation 1, with the key difference of a
weighting term based on the error set from closed-loop simulation in Stage 2.
In practice, although statistically equivalent, we upsample scenes by a fixed
factor rather than reweighting, as it is known to be more stable and robust
[14].

By training a policy using CW-ERM, we expect it to upsample scenes that
perform poorly in closed-loop evaluation, making the policy network robust to
the covariate shift seen during inference time while unrolling the policy.

We describe the complete CW-ERM training procedure in Algorithm 1 and in
Figure 1 we show a high-level overview of our method.

Input:Training set $\mathcal{D}$and hyperparameters $K$(number of epochs) and
$w$(upsampling factor). Stage 1: Identification policy1\. Train
$\hat{\pi}_{\text{ERM}}$on $\mathcal{D}$using ERM for $K$epochs (Equation 1).
Stage 2: Closed-loop simulation2\. Perform closed-loop simulation of the
policy $\hat{\pi}_{\text{ERM}}$in training scenes; 3\. Compute closed-loop
evaluation metrics; 4\. Build the error set $E_{\hat{\pi}_{\text{ERM}}}$of
training scenes from closed-loop metrics (Equation 2). Stage 3: Final
policy5\. Construct upsampled dataset $\mathcal{D}_{\text{up}}$by upsampling
the error set $E_{\hat{\pi}_{\text{ERM}}}$by $w$times; 6\. Train final model
$\hat{\pi}_{\text{CW-ERM}}$for $K$epochs on $\mathcal{D}_{\text{up}}$via CW-
ERM (Equation 3).

Algorithm 1 CW-ERM training procedure

###  II-C Relationship to covariate shift adaptation with density ratio
estimation

One important connection of our method is with covariate shift correction
using density ratio estimation [11]. To correct for the covariate shift, the
negative log-likelihood is often weighted by the density ratio $r(s)$:

$\operatorname*{argmin}_{\pi\in\Pi}\mathbb{E}_{s\sim
d_{\pi^{*}},a\sim\pi^{*}(s)}[r(s)\ell(s,a,\pi)]$ (4)

where $r(s)$ is defined as the density ratio between test and training
distributions:

$r(s)=\frac{p_{\text{test}}(s)}{p_{\text{train}}(s)}$ (5)

In practice, $r(s)$ is difficult to compute and is thus estimated. The density
ratio will be higher when the sample is more important for the test
distribution. In our method (CW-ERM), instead of using the density ratio to
weight training samples, we resample the training set based on an estimate of
each data point’s importance towards good closed-loop behaviours. Like the
density ratio, the weighting in our case will also be higher for when the
sample is important for the test distribution.

One key characteristic of the importance weighted estimator is that it can be
consistent even under covariate shift. We leave, however, the analysis of
theoretical properties of our approximation for future work.

##  III Related Work

Closely related to our work is the “Learning from Failure” method [15] (also
known as LfF), where the authors train two models at the same time with a
similar purpose of mitigating bias. The difference is that in CW-ERM we train
models sequentially and we use closed-loop evaluation metrics instead of the
loss to upsample, which permits the use of non-differentiable metrics. We also
do not use GCE (generalized cross-entropy) to bias the identification model.
Our method is simple from a training perspective, but unlike LfF, it does
assume the availability of a simulator.

A similar approach that has been successfully applied to computer vision and
natural language processing is the _Just Train Twice_ (JTT) algorithm [16].
JTT similarly first trains an identification model, then trains another model
which upweights samples misclassified by the identification model. Although
similar in the sense that two models are trained sequentially, in JTT the goal
is to deal with worst-group accuracy and not to improve robustness to closed-
loop behaviors of a planning model, as in our case.

Several works have attempted to address the covariate shift problem in
imitation learning. ChauffeurNet [2] and SafetyNet [7] add state perturbation
to the training data for improved generalization. Similarly, DAVE-2 [4] used
video captured from three different cameras as well as perturbation on the
captured images. Another common approach is to supplement training with on-
policy data [17, 18, 19], however, in practice, collecting on-policy data for
use during training can be extremely expensive and time-consuming.

Similar to our work, Urban Driver [20] also utilizes a closed-loop simulator,
but the simulator is used directly during training to generate unrolls while
using BPTT (backpropagation through time). Urban Driver needs a differentiable
simulator and does not scale well due to the need for rollouts during training
and the memory requirements of BPTT. In contrast, our work takes a much
simpler approach where the closed-loop simulator is only used to identify
which samples to up-weight and does not require a differentiable simulator,
while being able to directly identify scenes that are important for closed-
loop evaluation without having to change the training loss to add
differentiable collision losses as in Urban Driver [20]. In our work, any
closed-loop metric can be used, with no requirements for differentiability.

##  IV Experimental Evaluation

###  IV-A Policy network architecture

Our method is agnostic to model architecture choices. To evaluate our CW-ERM
approach, we adopt the recent network architecture of [7] to represent a
strong baseline performance for SDV planning. This model uses a transformer-
based [21] architecture with a vectorial input representation [22] to create
features for each element into vector sets. It consists of a PointNet-based
[23] module for local processing of vectorized inputs and a global graph using
a Transformer encoder for reasoning about interactions with agents and map
features. Differently from [7], we don’t use a safety layer, as we want to
evaluate the planner performance without external trajectory fallbacks. For
further details, please refer to [7].

###  IV-B Training

During training, we found that stopping training of the identification policy
before convergence, similar to what was done in JTT [16] and LfF [15], also
yielded better results. We limited the capacity of the identification policy
by training it until $K$ epochs. The insight is that important biases are
learned in early training phases [15] and limiting model capacity can avoid
overfitting and avoid depletion [16] of the error set used for the training of
the final policy. We train the final policy for 40 epochs.

For the ERM baseline, we compare our method against two different experiments
where we have the traditional BC trained with ERM with and without
perturbations (details can be found in Appendix -A).

###  IV-C Datasets

We train and test CW-ERM on a proprietary dataset. Our SDV data is collected
in challenging urban missions on San Francisco and Palo Alto roads. This
dataset is a collection of driving trajectories from our SDV and surrounding
agents, along with recorded HD Maps. Various types of behavioral scenarios in
urban driving such as stopping behind a lead vehicle, stopping at
intersections, and driving among dense cars, pedestrians, cyclists etc. are
captured. The majority of scenes in our dataset are between 11-13 seconds
long, with the longest lasting up to 30 seconds. The total data used during
training is 180 hours and we validate and test on 60 hours of driving data
each.

###  IV-D Evaluation framework

We compute the closed-loop evaluation metrics by doing rollouts of the policy
in the log-replayed scenes on a simulator 111Please refer to the
Reproducibility section for details on the open-sourcing of the simulator and
metrics used in this work.. During the unroll, trajectories are recorded. An
evaluation plan composed of a set of metrics and constraints is executed over
the recorded trajectories. We count every scene that violated a constraint
(e.g., a collision) and then compute the confidence intervals (CIs) for each
metric using a Binomial exact posterior estimation with a flat prior, which
gives similar results (up to rounding errors) to bootstrapping as recommended
in [24].

TABLE I: Experimental results from closed-loop evaluation in simulation. In
this table we show a baseline method of behavioral cloning (ERM) with and
without perturbations (details can be found in Appendix -A) together with the
results from single and multi-metric experiments. Lower is better for all
metrics.

Method |  Upsampled metric |  Perturbation |  Front Collisions |  Side Collisions |  Rear Collisions |  Dist. to ref. traj.
---|---|---|---|---|---|---
ERM (baseline) | (not applicable) |  | 67 (52.8, 85.0) | 97 (79.6, 118.2) | 114 (94.9, 136.8) | 75 (59.9, 93.9)
ERM (baseline) | (not applicable) | ✓ | 14 (8.4, 23.5) | 55 (42.3, 71.6) | 34 (24.4, 47.5) | 35 (25.2, 48.6)
CW-ERM (ours) | Front Collisions | ✓ |  9 (4.8, 17.1) | 52 (39.7, 68.1) | 38 (27.7, 52.1) | 39 (28.6, 53.2)
CW-ERM (ours) | Side Collisions | ✓ | 11 (6.2, 19.7) | 47 (35.4, 62.5) | 31 (21.9, 44.0) | 35 (25.2, 48.6)
CW-ERM (ours) | Dist. to Reference Trajectory | ✓ | 10 (5.5, 18.4) | 50 (37.9, 65.8) | 40 (29.4, 54.4) |  28 (19.4, 40.4)
CW-ERM (ours) | Front + Side Collisions | ✓ | 9 (4.8, 17.1) | 55 (42.3, 71.6) | 39 (28.5, 53.2) | 35 (25.2, 48.6)
CW-ERM (ours) | Front + Side + Rear Collisions | ✓ | 11 (6.2, 19.7) | 50 (37.9, 65.8) |  28 (19.4, 40.4) | 37 (26.9, 51.0)
CW-ERM (ours) | Front + Side + Dist. to ref. traj. | ✓ | 10 (5.5, 18.4) |  46 (34.5, 61.3) | 40 (29.4, 54.4) | 31 (21.9, 44.0)

Figure 2: Visual representation of the experimental results from closed-loop
evaluation in simulation shown in Table I. Confidence intervals (CIs) were
calculated using .95 interval from an exact Binomial posterior with a flat
prior. In this plot we only compare against the best baseline (ERM with
perturbation).

###  IV-E Metrics

Metrics computed in the closed-loop simulator are used to construct the error
set (Equation 2). In our evaluation we consider certain important metrics: the
number of front collisions, side collisions, rear collisions, and distance
from reference trajectory. The distance from reference trajectory considers
the entire target trajectory for the current simulated point. A failed scene
with respect to this metric is one where the distance of the simulated center
of the SDV to the closest point in the target trajectory is farther than four
meters.

In our evaluation, we perform two sets of experiments: _single metric_ and
_multi metric_. In single metric experiments we construct the error set using
only a single metric, while for multi metric we use scenes from multiple
metrics together.

###  IV-F Results

####  IV-F1 Single Metric

We show the results from single metric experiments in Table I. We can see that
the number of collisions significantly reduced for both side and front
collision experiments. We found improvements in the range of $\sim$35% on the
test set for some metrics when compared to the baseline.

We also found that the largest margin of improvements targeting single metrics
in isolation were seen when using single metric based error set, while a
balance was achieved when targeting multiple metrics, which suggests a Pareto
front of solutions when targeting multiple objectives.

Variance is also lower in some cases when compared to the baseline. We note
that while upsampling a certain metric, it shows noticeable improvements in
other related metrics. For example, in our single metric experiments, we see
that improving side collisions also improve rear collisions. This is evidence
that the model is not only getting better at side collisions but also becoming
less passive (as indicated by reduction in rear collisions, due to log-
replayed agents in simulation that are non-reactive).

Qualitative results during closed-loop unroll are shown in Figures 4 and 5.
Here, we show improved behavior when the scene is upsampled in CW-ERM - once
for front collisions Figure 4 and another for side collisions Figure 5. Here,
we see a better response in the CW-ERM model which avoids collisions by
waiting at intersection and slowing down next to a lead vehicle.

####  IV-F2 Multi Metric

In our multi-metric experiments, we combine two or more metrics - namely
$m_{1},m_{2}..m_{N}$ \- into a single upsampling experiment. The metrics are
equally weighted and hence scenes that fail due to any $m_{i}$ will be added
to the error set. While improvements are noticeable upon combining Front and
Side collisions or Front, Side and Distance to the reference trajectory in
Table I, considerable regression is observed when adding rear collisions. As
we can see from the experiments, this is clearly related to the amount of
false-positives (FPs) in rear collisions due to the lack of agent reactivity
during log playback in the simulator.

Figure 3: A sample plot showing the effect of upsampling factor on the
performance based on single cost upsampling. Here, the scenes are upsampled
based on side collisions by factor $w$ along X-axis. The resulting side
collisions obtained while evaluating on the validation dataset is obtained on
the Y-axis. It is evident from the plot that the performance improves until an
upsampling factor of $w=50$ after which number of side collisions begin to
increase.

###  IV-G Hyper-parameter Tuning / Upsampling Experiments

We evaluate the performance of the upsampled training set using various
identification models $K\in\\{10,20\\}$ with various upsampling factors
$w\in\\{10,20,30,50\\}$ on the validation dataset. Single metric upsampling
experiments responded better to error set extracted from the training where
$K=10$, while using $K=20$ performed better for multi-metric experiments. We
find that the size of the resulting upsampled error set influences
performance. As seen from Table I and Figure 3, there exists a limit beyond
which performance does not improve during the upsampling experiments. A
similar observation was noted in JTT [16], where they also found an upsampling
factor for which beyond it, worst-group accuracy could not be improved.

Figure 4: A scene from the test dataset showing the behavior of ego in
Baseline (ERM) (with perturbation) and CW-ERM front collisions upsampled. The
blue dots are the target trajectory and the yellow dots are the predicted
trajectory. The ego is the box in pink and blue-green boxes are other agents.
Here, we see that the baseline moves ahead at intersection ignoring the car
from the right resulting in a front collision. In contrast, for the CW-ERM
policy, it waits at the intersection. Figure 5: A scene from the test dataset
showing the behavior of ego in Baseline (ERM) (with perturbation) and CW-ERM
side collisions upsampled. The blue dots are the target trajectory and the
yellow dots are the predicted trajectory. The ego is the box in pink and blue-
green boxes are other agents. Here, we see that the baseline deviates from
target trajectory and collides with the bus ahead. The CW-ERM policy slows
down behind the bus and prevents a collision.

##  V Limitations

Although our method is efficient, easy to implement and showed significant
improvements, it also comes with limitations that are important to highlight.
In this work, we use a proprietary dataset for evaluation, primarily due to
the current lack of available closed-loop evaluation benchmarks. Most public
datasets available today are focused on agent prediction tasks and on open-
loop metrics (e.g. [25, 26]). Recently, the closed-loop planning benchmark
nuPlan [27] was released, but is still under active development and requires a
special license for industrial labs to use.

Our method also introduces two new hyper-parameters: $K$ (number of epochs for
the identification model) and $w$ (upsampling factor), however, we found the
improvement to be robust to different parameterizations of these parameters,
similar to past observations [16].

We performed log-replay of agents in simulation, which can produce false
positives for the rear collision metric. We leave further analysis and the
usage of a reactive simulator as future direction.

We have not yet deployed our policy in a real-world SDV, however, we evaluated
it on a closed-loop evaluation framework that is known to be closely
representative of real-world performance. Deployment and testing of a policy
in a real-world SDV on public roads requires further safety evaluations that
we leave as future work.

##  VI Reproducibility

We make available our closed-loop simulator and the closed-loop metrics used
in this work in the following open-source repository:

https://woven-planet.github.io/l5kit

##  VII Discussion

Most recent improvements in imitation learning are based on improving the
asymptotic performance of algorithms. In this work we showed a different
direction that tackles the problem by directly addressing the mismatch between
training and inference without requiring an extra human oracle or adding extra
complexity during training. Our method is as simple as upsampling scenes by
leveraging any existing simulator and training two models, yet it showed that
there is still room for significant improvements without having to deal with
human-in-the-loop, training rollouts or impacting the policy inference
latency. We also described an important potential connection of our method
with density ratio estimation for covariate shift correction [11], which we
believe is an exciting future research direction that could provide better
theoretical understanding of the improvements seen in our experiments.

###  -A Perturbation

In our experiments, we employed similar perturbation techniques as in [7]. We
randomly perturb the ego’s current state to shift from the current trajectory
in some of the training examples. To put it more concretely, for the perturbed
states, we add a zero-mean Gaussian noise to the ego’s current position and
heading. We perturb the ego’s speed by $av+|b|$ where $v$ is the current
speed, $a$ and $b$ are speed multiplier and bias terms generated by a zero-
mean Gaussian distribution. Taking the absolute value of the bias term is to
ensure that the perturbed speed is always non-negative. Additionally, we
perform collision checks for every perturbed state and we do not include any
perturbed states which includes collisions.

###  -B Policy network hyper-parameters

To train the baselines and our policy we employ a distributed training with a
local batch size of 64 for each replica (with an effective batch size of 4096
when using 64 GPU replicas). We use a learning rate of 0.001 that is annealed
with cosine schedule during training for 40 epochs (except for the
identification models as described in the Section IV). We also used MAE (mean
absolute error) loss and the Adam [28] optimizer with default PyTorch [29]
parameters for $\beta_{1}$ and $\beta_{2}$.

###  -C Scene distribution

We select training, validation and test data subsets to be balanced. The
scenarios selected for the datasets are diverse, and a variety of complex
urban scenarios have been curated such that it is possible to evaluate closed-
loop performance, even on minority scenario groups. In Figure 6, a detailed
split of the various scenarios is provided.

Figure 6: Scenario distribution of the training dataset. Scenes can be of
single or multiple scenarios - in the case of multiple scenario tags, the
scenario names are comma separated.

## ACKNOWLEDGMENT

We would like to thank Kenta Miyahara, Nobuhiro Ogawa and Ezequiel Castellano
for the review of this work and everyone from the UK Research Team and the ML
Planning team who supported this work through the ecosystem needed for all
experiments, and for the fruitful discussions.

## References

  * [1] D. A. Pomerleau, “Alvinn: An autonomous land vehicle in a neural network,” in _Advances in Neural Information Processing Systems_ , vol. 1, 1989.
  * [2] M. Bansal, A. Krizhevsky, and A. Ogale, “Chauffeurnet: Learning to drive by imitating the best and synthesizing the worst,” in _Proceedings of Robotics: Science and Systems_ , June 2019.
  * [3] H. Xu, Y. Gao, F. Yu, and T. Darrell, “End-to-end learning of driving models from large-scale video datasets,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 3530–3538.
  * [4] M. Bojarski, D. Del Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang, X. Zhang, J. Zhao, and K. Zieba, “End to end learning for self-driving cars,” 2016.
  * [5] F. Codevilla, M. Müller, A. López, V. Koltun, and A. Dosovitskiy, “End-to-end driving via conditional imitation learning,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 4693–4700.
  * [6] A. Kuefler, J. Morton, T. Wheeler, and M. Kochenderfer, “Imitating driver behavior with generative adversarial networks,” in _2017 IEEE Intelligent Vehicles Symposium (IV)_ , 2017, pp. 204–211.
  * [7] M. Vitelli, Y. Chang, Y. Ye, A. Ferreira, M. Wołczyk, B. Osiński, M. Niendorf, H. Grimmett, Q. Huang, A. Jain, _et al._ , “Safetynet: Safe planning for real-world self-driving vehicles using machine-learned policies,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 897–904.
  * [8] F. Codevilla, E. Santana, A. M. López, and A. Gaidon, “Exploring the limitations of behavior cloning for autonomous driving,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 9329–9338.
  * [9] P. de Haan, D. Jayaraman, and S. Levine, “Causal confusion in imitation learning,” in _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds., vol. 32. Curran Associates, Inc., 2019.
  * [10] J. C. Spencer, S. Choudhury, A. Venkatraman, B. D. Ziebart, and J. A. Bagnell, “Feedback in imitation learning: The three regimes of covariate shift,” _CoRR_ , vol. abs/2102.02872, 2021.
  * [11] H. Shimodaira, “Improving predictive inference under covariate shift by weighting the log-likelihood function,” _Journal of Statistical Planning and Inference_ , vol. 90, no. 2, pp. 227–244, 2000.
  * [12] V. Vapnik, “Principles of risk minimization for learning theory,” in _Advances in Neural Information Processing Systems_ , 1991.
  * [13] S. Ross and D. Bagnell, “Efficient reductions for imitation learning,” in _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, Y. W. Teh and M. Titterington, Eds., vol. 9. Chia Laguna Resort, Sardinia, Italy: PMLR, 13–15 May 2010, pp. 661–668.
  * [14] J. An, L. Ying, and Y. Zhu, “Why resampling outperforms reweighting for correcting sampling bias with stochastic gradients,” in _9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021_. OpenReview.net, 2021.
  * [15] J. Nam, H. Cha, S. Ahn, J. Lee, and J. Shin, “Learning from failure: De-biasing classifier from biased classifier,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 20 673–20 684, 2020.
  * [16] E. Z. Liu, B. Haghgoo, A. S. Chen, A. Raghunathan, P. W. Koh, S. Sagawa, P. Liang, and C. Finn, “Just train twice: Improving group robustness without training group information,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 6781–6792.
  * [17] S. Ross, G. Gordon, and D. Bagnell, “A reduction of imitation learning and structured prediction to no-regret online learning,” in _Proceedings of the fourteenth international conference on artificial intelligence and statistics_. JMLR Workshop and Conference Proceedings, 2011, pp. 627–635.
  * [18] Y. Pan, C.-A. Cheng, K. Saigol, K. Lee, X. Yan, E. A. Theodorou, and B. Boots, “Imitation learning for agile autonomous driving,” _The International Journal of Robotics Research_ , vol. 39, no. 2-3, pp. 286–302, 2020.
  * [19] A. Prakash, A. Behl, E. Ohn-Bar, K. Chitta, and A. Geiger, “Exploring data aggregation in policy learning for vision-based urban autonomous driving,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11 763–11 773.
  * [20] O. Scheel, L. Bergamini, M. Wolczyk, B. Osiński, and P. Ondruska, “Urban driver: Learning to drive from real-world demonstrations using policy gradients,” in _Conference on Robot Learning_. PMLR, 2022, pp. 718–728.
  * [21] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in Neural Information Processing Systems_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds., vol. 30. Curran Associates, Inc., 2017\.
  * [22] J. Gao, C. Sun, H. Zhao, Y. Shen, D. Anguelov, C. Li, and C. Schmid, “Vectornet: Encoding hd maps and agent dynamics from vectorized representation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11 525–11 533.
  * [23] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in _CVPR_ , 2017.
  * [24] R. Agarwal, M. Schwarzer, P. S. Castro, A. C. Courville, and M. Bellemare, “Deep reinforcement learning at the edge of the statistical precipice,” _Advances in Neural Information Processing Systems_ , vol. 34, 2021.
  * [25] M.-F. Chang, J. Lambert, P. Sangkloy, J. Singh, S. Bak, A. Hartnett, D. Wang, P. Carr, S. Lucey, D. Ramanan, and J. Hays, “Argoverse: 3d tracking and forecasting with rich maps,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019.
  * [26] S. Ettinger, S. Cheng, B. Caine, C. Liu, H. Zhao, S. Pradhan, Y. Chai, B. Sapp, C. R. Qi, Y. Zhou, Z. Yang, A. Chouard, P. Sun, J. Ngiam, V. Vasudevan, A. McCauley, J. Shlens, and D. Anguelov, “Large scale interactive motion forecasting for autonomous driving: The waymo open motion dataset,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , October 2021, pp. 9710–9719.
  * [27] K. T. e. a. H. Caesar, J. Kabzan, “Nuplan: A closed-loop ml-based planning benchmark for autonomous vehicles,” in _CVPR ADP3 workshop_ , 2021.
  * [28] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _International Conference on Learning Representations (ICLR)_ , 2015.
  * [29] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in _Advances in Neural Information Processing Systems 32_. Curran Associates, Inc., 2019, pp. 8024–8035.


# Is solid copper oxalate a spin chain or a mixture of entangled spin pairs?

Pavel Pokhilkoa, Dmitry S. Bezrukovb,c, and Anna I. Krylova
a Department of Chemistry, University of Southern California, Los Angeles,
California 90089-0482
b Department of Chemistry, M.V. Lomonosov Moscow State University, Moscow
119991, Russia
c Skolkovo Institute of Science and Technology, Skolkovo Innovation Center,
Nobel str. 3, Moscow 143026, Russia

###### Abstract

Macroscopic assemblies of interacting spins give rise to a broad spectrum of
behaviors determined by the spatial arrangement of the magnetic sites and the
electronic interactions between them. Compounds of copper (II), in which each
copper carries spin $\frac{1}{2}$, exhibit a vast variety of physical
properties. For antiferromagnetically coupled spin sites, there are two
limiting scenarios: spin chains in which the spins can exhibit a long-range
order or a mixture of dimers in which the spins within each pair are entangled
but do not communicate with the spins from other dimers. In principle, the two
types can be distinguished on the basis of experimental observations and
modeling using empirically parameterized effective Hamiltonians, but in
practice, ambiguity may persist for decades, as is the case for copper
oxalate. Here we use high-level _ab initio_ calculations to establish the
validity of the nearest-site Heisenberg model and to predict the interaction
strength between the magnetic sites. The computed magnetic susceptibility
provides an unambiguous interpretation of magnetic experiments performed
throughout half a century, clearly supporting the infinite spin-chain behavior
of solid copper oxalate.

##  I Introduction

Figure 1: An assembly of interacting spins is described by the Heisenberg
Hamiltonian. Depending on the structure and strength of interactions, such
system can exhibit vastly different physics. For one-dimensional systems, the
two limiting cases are spin chains in which the spins can exhibit a long-range
order or a mixture of dimers in which the spins within each pair are entangled
but do not communicate with the spins from other dimers.

Strongly correlated materials exhibit unusual properties, which are exploited
in emerging applications such as quantum information science and
spintronicsWunderlich2016:review:antif_spintronics . Physical properties
depend on the spatial arrangement and the strengths of electronic couplings
between the magnetic sites, as illustrated in Fig. 1. For example, chains of
antiferromagnetically coupled spins exhibit long-range spin correlations and a
zero spectral gap at zero magnetic field. In contrast, an assembly of weakly
coupled spin pairs has no long-range order and has a gap at zero field. Both
cases give rise to a bell-shape magnetic susceptibility curve $\chi(T)$. In
both types of systems, the equilibrium state can be converted to a fully spin-
polarized state above saturation field, roughly equal to the exchange
coupling. The differences between the two cases can be seen at low
temperatures, e.g., spin chains have finite $\chi(T)$, in contrast to dimers.
The spectral gap and long-range order can also be probed by neutron scattering
experiments, e.g., dynamic structural factors for the spin-chains at zero
field start at zero energy and show characteristic sinusoidal shape, following
the spinon dispersion relation $\epsilon(k)=\frac{\pi}{2}|J|\sin k$, as was
illustrated for copper sulfateMourigal:CopperSulfate:spinon:2013 ; whereas the
dimers are expected to show flat bands at finite energy. In addition to their
technological relevanceWunderlich2016:review:antif_spintronics , spin chains
are also fundamentally interesting—for example, they give rise to quantum
criticality around the saturation fieldSchofield:QuanumCriticality:05 ;
Lorenz:2017:Cu:spin_chain .

Copper (II) salts can behave both as antiferromagnetic spin
chainsLorenz:2017:Cu:spin_chain ; Mourigal:CopperSulfate:spinon:2013  and as
spin-paired dimersBowers:mag_sus:fitting:1952 ; CuAc:Neese:2011 . When
experimental data is limited to finite-temperature magnetic susceptibility, as
in the case of copper oxalate, one cannot confidently distinguish between the
two regimes without first-principle calculations. Here we report the first
fully ab initio determination of the effective Hamiltonian and macroscopic
magnetic susceptibility for copper oxalate. The results provide unambiguous
interpretation of magnetic experiments performed throughout half a century and
clearly support the infinite spin chain behavior of the solid copper oxalate.

As shown in Fig. 1, an infinite chain of spins can be treated with the XXX
Heisenberg Hamiltonian

$\displaystyle H=-J\sum_{i}\mathbf{S}_{i}\mathbf{S}_{i+1},$ (1)

where the $\mathbf{S}_{i}$ are local spins and $J$ is the effective exchange
constant. If the effective exchange constant is negative, the system adopts
the antiferromagnetic singlet ground state with opposite spin orientation of
the adjacent magnetic centers. The case of $S=1/2$ is one of the few known
quantum integrable models. Its exact solution is given by the Bethe
ansatzBethe:ansatz:1931 , which facilitated the development of solid-state
physics and mathematics of integrable systems. The thermodynamic properties of
this model are well-known and have been used to explain experimental
observables of real materials, containing, for instance,
copperLorenz:2017:Cu:spin_chain ; Uchida:1996:Cu:spin_chain ,
vanadiumJacobson:1987:V:spin_chain , and magnesiumMorosin:1973:Mn:spin_chain
magnetic centers.

Despite its simplicity, the magnetic structure of copper oxalate has not been
settled. On the basis of susceptibility measurements, EPR, and EXAFS, both a
dimerCu:oxalate:1966 ; Dubicki:Cu:chains:1966 ; Grivel:CuC2O4:cryst  and an
infinite spin chain Heisenberg modelsSoos:CuC2O4:EPR:1976 ;
Goulon:Cu:EXAFS:1979  have been proposed and used to fit the experimental
data. Even after the determination of the crystal structure in 2014, the
choice of magnetic model remained openGrivel:CuC2O4:cryst .

Even when the model is known, the first-principle determination of the
effective exchange constant is not trivial. The strongly correlated nature of
magnetic systems combined with their large sizes, pose a formidable challenge
for quantum chemistry. Here we use state-of-the-art wave-function methods
based on equation-of-motion coupled-cluster (EOM-CC) theoryKrylov:EOMRev:07 .
In contrast to popular broken-symmetry density functional theory, our approach
has no empirical or system-dependent parameters and does not rely on
unphysical spin-scrambled solutions. We begin with the full-electron treatment
of model systems and use the resulting wave functions to build effective
HamiltoniansMarlieu:MagnetRev:2014 ; Pokhilko:EffH:2020 , which afford a
coarse-grained description of the electronic structure and allow extrapolation
to infinite systems. This is the first application of the spin-flip EOM-CC
method to describe a periodic system, opening a new route in the treatment of
periodic strongly correlated systems. Our theoretical results provide
unambiguous interpretation of the magnetic measurementsCu:oxalate:1966 ;
Dubicki:Cu:chains:1966 ; Soos:CuC2O4:EPR:1976 ; Goulon:Cu:EXAFS:1979 ;
Grivel:CuC2O4:cryst  of copper oxalate performed throughout several decades.

##  II Results and discussion

Figure 2: Top, left: Crystal structure of copper oxalate (CuC2O4). Top, right:
fragments with 2, 3, and 4 copper atoms from the crystal. Middle: Localized
orbitals defining the open-shell model spaces. Bottom: Des Cloizeaux effective
Hamiltonians built from the EOM-SF wave functions in the basis of open-shell
determinants. The energies are shifted to produce a zero trace. The numbers in
blue show the couplings between nearest neighbors. The numbers in red show the
couplings between distant centers.

The crystal structureGrivel:CuC2O4:cryst  of copper oxalate, shown in Figure
2, reveals a regular pattern of copper oxalate chains with a bidentate
orientation of the oxalate ligands. We begin with EOM-SF-CCSD/cc-pVDZ
calculations for the fragments of increasing length (containing 2, 3, and 4
copper centers, as shown in Figure 2, and capped with hydrogen atoms). We use
the resulting energies and wave functions to build effective
HamiltoniansMarlieu:MagnetRev:2014 ; Pokhilko:EffH:2020 , which afford a
coarse-grained description of the electronic structure.

The effective Heisenberg Hamiltonians for fragments built from these highly
accurate EOM-SF-CCSD calculations are shown in the middle panel of Figure 2.
Bloch’s and des Cloizeaux’s versions of the effective Hamiltonians are
consistent with each other and yield nearly the same values of effective
exchange constants $J$, as shown in Table S2 in the SI. This procedure
provides a rigorous mapping between the full-electron description and
effective interactions between the magnetic centers. In particular, the
resulting effective exchange constant includes not only the contribution from
the exchange integral, but also screening and correlation effects, folded in
by means of the Bloch wave operator. We then use the Heisenberg Hamiltonian to
obtain the states that are not reachable by a single spin-flip, as was
pioneered by Mayhall and Head-Gordon for single-molecule
magnetsMayhall:1SF:2015 .

Note that in these calculations we constructed the full Heisenberg Hamiltonian
(as shown in the pink box in Fig. 1), i.e., without assuming nearest-neighbor
approximation or constraining all neighboring $J$ to be the same. The analysis
of these ab initio constructed effective Hamiltonians clearly reveals
antiferromagnetic spin interactions. The off-diagonal elements decay rapidly
with the distance between the copper centers, such that only the nearest-
neighbor spin couplings are significant. Non-nearest-neighbor couplings are
about 1 cm-1 or less, which is comparable with the thresholds used in the
calculations. To further quantify the impact of neglecting the non-nearest
neighbors contributions, we carried out two additional numeric tests. First,
we computed eigenvalues of the Heisenberg Hamiltonian in which we zeroed the
contributions from non-nearest neighbors and compared them with the
eigenvalues of the full Hamiltonian; the results are shown in Table S3 of the
SI. As one can see, the individual eigenvalues change by less than 0.5 cm-1
and the effect on the spectral gap is of the same magnitude—for example, the
difference between the highest and the lowest eigenstate in the 4-copper
Hamiltonian changes by 0.44 cm-1 due to the zeroing out non-nearest neighbor
elements, which is indeed small compared to the value of the gap (301.07
cm-1). Second, we tri-diagonalized the full effective Hamiltonian for the 3
and 4 copper centers, so that the effect of non-nearest neighbor contributions
is folded into the remaining matrix elements. The results, which are shown in
Fig. S2 in the SI, show that this operation changes the value of the key
matrix element (between the two middle coppers) by only 0.63 cm-1. Therefore,
we conclude that the nearest-neighbor approximation is indeed fully justified
for this system, such that the full Heisenberg Hamiltonian for the fragments
can be replaced by the XXX Heisenberg Hamiltonian for the infinite spin chain,
Eq. (1), which we use to extrapolate the Hamiltonian to the limit of the
infinite number of copper centers. The effective exchange constants computed
for the fragments of increasing sizes (collected in Table S2) show little
variation and converge rapidly with respect to the model system size. Hence,
for calculations of thermodynamic properties we use the $|J|$ constant from
the middle of the 4-copper fragment (-178.1 cm-1).

To account for weak correlations beyond EOM-CCSD, we calculated a perturbative
triples correction using the (fT) model for a model dimer. Surprisingly,
inclusion of triple excitations leads to a more than 20% increase of $|J|$
(42.7 cm-1), illustrating that weak correlations are important for
quantitative results. The direction of the change can be explained by the
Pauli repulsion principle. Pauli repulsion, which is already built-in in the
wave-function at the mean-field level by virtue of using Slater determinants,
allows the electrons in the triplet state to avoid each other. The resulting
Pauli hole covers the errors due to an incomplete description of the Coulomb
hole, which requires electron correlation. Consequently, the correlation
effects are always smaller for triplet statessfpaper  and improving
correlation treatment results in stabilizing singlets relative to triplets.
Previous calculations of dicopper single-molecule magnets suggested that the
cc-pVDZ basis is nearly sufficientOrms:magnets:17 ; Pavel:OSFNO:2019 . To
account for a basis-set effect beyond cc-pVDZ, we performed EOM-SF-CCSD
calculations for the dimer with the cc-pVTZ basis set. These calculations show
that this improvement of the basis-set reduces the $|J|$ value by 8.5 cm-1.

To account for spin–orbit interaction, which can split and shift energy
levels, we computed spin–orbit couplings (for the model dimer) using the
Breit–Pauli Hamiltonian within spin–orbit mean-field approximation and non-
relativistic EOM-CC states. In these calculations we used the EOM-DIP-CCSD
ansatz, which can describe not only the lowest singlet and triplet states
(i.e., the Heisenberg states), but the entire manifold of the low-lying states
derived from different occupations of coppers’ $d$-orbitals. Although these
states are relatively far in energy from the lowest singlet/triplet pair
(about 1.6 eV above the lowest singlet), they are expected to be important for
spin–orbit coupling by virtue of El-Sayed’s ruleElSayedRule-orig ; ElSayedRule
; Pavel:SOCNTOs:2019 , which says that states with different orbital
occupancies give rise to large spin–orbit interactions.

The EOM-DIP-CCSD ansatz, which can handle strong correlation well but
underestimates weak correlation and orbital relaxation, overestimates the
energy gap relative to EOM-SF-CCSD ($J=-388.6$ cm-1 for the dimer). To correct
for this limitation, we combined the EOM-SF-CCSD energy gap with EOM-DIP-CCSD
spin–orbit couplings to evaluate the effective $|J|$. Table S4 illustrates the
convergence of the $J$ constant with respect to the the lowest singlet number
of interacting states. We observe that both the lowest singlet and the lowest
triplet states interact strongly through the spin–orbit operator with higher
excited states. This interaction shifts both states, largely canceling the
effect of the spin–orbit coupling on the singlet-triplet gap. Nevertheless,
the singlet interacts more strongly with the excited states than the triplet
does, which results in a substantial (13 %) increase of antiferromagnetic
$|J|$ by $32$ cm-1. We note that these results indicate that neglect of the
spin–orbit interaction in previous studies of di-copper single-molecule
magnets might be responsible for the reported systematic discrepancies between
the theoretical and computed $J$-valuesOrms:magnets:17 .

Table 1 summarizes the various contributions to the effective $J$ and gives
our best estimate of $J=-244$ cm${}^{-1}=30.25$ meV. Comparing to other copper
spin chains, this value is 33 times higher than that of copper
pyrazineLorenz:2017:Cu:spin_chain  and 118 times higher than of copper
sulfateMourigal:CopperSulfate:spinon:2013 , implying that a much stronger
magnetic field (on the order of 400-500 T) would be required to bring this
system to the fully polarized spin state and into the quantum criticality
regime.

Table 1: Contributions to the effective exchange constant $J$, cm-1.  $-178.1$ | Converged effective Hamiltonian
---|---
$-42.7$ | Perturbative triples
$-32$ | Spin–orbit shift
$+8.5$ | Basis-set effects
$-244$ | Total

To compare the computed effective exchange constants with the experimental
observable, one needs to compute macroscopic magnetic susceptibility
$\chi(T)$. Even with the solvable Bethe ansatz, such calculations require
additional approximations. Here we consider two approaches: analytic classical
and numeric quantum treatment of the temperature dependence of $\chi(T)$. In
the classical treatment, the quantum Heisenberg Hamiltonian is mapped onto the
classical Heisenberg HamiltonianFisher:class:spin_chain:1964 . The connection
between the quantum and classical parameters in the case of $S=1/2$ is given
by

$\displaystyle J^{quant}/2=J^{cl},$ (2) $\displaystyle g^{quant}=g^{cl}.$ (3)

The analytic solution of the $\chi(T)$ exists for the case of a classical
linear spin chain with nearest-neighbor approximation under the assumption of
uniformness, meaning that $J_{i,i+1}$ and $g_{i}$ are the same for all local
spins. The solution has the following formFisher:class:spin_chain:1964 :

$\displaystyle\chi(T)=\frac{N(g^{cl})^{2}\mu_{B}^{2}}{12kT}\frac{1+u(K)}{1-u(K)},$
(4) $\displaystyle u(K)=\coth K-(1/K),$ (5) $\displaystyle
K=\frac{J^{cl}/2}{kT},$ (6)

where $N$ is the number of centers, $g^{cl}$ is a classical electron
$g$-factor (we assumed $g^{cl}=g^{quant}=2$), $\mu_{B}$ is Bohr’s magneton,
and $k$ is the Boltzmann constant. This expression predicts a maximum on the
susceptibility curve at:

$\displaystyle kT^{\text{cl}}_{max}\approx 0.2382|J^{cl}|.$ (7)

The analytic expression for quantum susceptibility is not known, but there are
Padé approximantsJohnston:pade_approx:2000  fitted to numerical curves with
high quality. We used such a curve (Fit 1 in Table I in the Ref.
Johnston:pade_approx:2000 ) for calculation of quantum magnetic
susceptibility. The susceptibility curve of a quantum infinite spin chain has
a maximum at

$\displaystyle kT_{max}\approx 0.640851|J|.$ (8)

One can immediately see that for the same value of $J^{quant}$, the classical
model yields Tmax, which is about 5.4 times lower that of the quantum model.
Conversely, using the classical model to fit the experimental data would yield
a 5.4 times larger effective exchange constant.

It is instructive to compare the predictions of the Heisenberg model with the
dimer model, which gives the well-known Bleaney–Bowers expression for magnetic
susceptibilityBowers:mag_sus:fitting:1952

$\displaystyle\chi^{\text{dim}}(T)=\frac{N_{A}g^{2}\mu^{2}_{B}}{kT}\frac{2}{3+e^{-\frac{J}{kT}}}$
(9)

The susceptibility of the dimer reaches maximum at

$\displaystyle kT^{\text{dim}}_{max}\approx 0.624|J|.$ (10)

Thus, although based on very different physics, the dimer model yields very
similar value of Tmax, i.e., within 3% from the quantum susceptibility derived
from the Heisenberg model.

Figure 3: Left: Classical and quantum susceptibilities of the XXX Heisenberg
model and the susceptibility for a dimer model (Bleaney–Bowers expression)
computed with $J=-244$ cm-1. Right: Quantum susceptibility computed with the
effective exchange constants from different approximations. $J_{0}$ is the
converged value of the exchange coupling constant from the 4-center fragment.
$J_{\text{tot}}=-244$ cm-1 is the overall value of the coupling that includes
$J_{0}$, perturbative triples, spin–orbit shift, and basis-set extrapolation.
The experimental data set 1 and 2 are from the Refs. Dubicki:Cu:chains:1966
and Zelentsov:Cu:oxalate:1964 , respectively.

The results of the classical and quantum Heisenberg models and the dimer model
are shown in Fig. 3. The susceptibility of the classical spin model differs
dramatically from the quantum spin model (Figure 3, left), making the
classical approximation not satisfactory in the interpretation of the magnetic
data. This is expected, since the susceptibility of the classical model
approaches the susceptibility of the quantum model only for the large spin on
each centerXiang:Spin_chains_large_S:1998 . Classical, quantum, and
experimental susceptibility curves reach maximum at 42 K, 225 K, and 260 K,
respectively. The quantum dimer model gives the maximum at 219 K. The
difference between the quantum infinite chain and the dimer Heisenberg models
in terms of the maximum position of the $\chi(T)$ curve is not large, but the
absolute values are vastly different (the dimer model yields absolute values
that are about 4 times higher). In the SI (Fig. S4), we show the comparison
between the quantum $\chi(T)$ and the scaled as well as scaled and shifted
dimer’s $\chi(T)$ curves. As one can see, if the dimer $\chi(T)$ is scaled to
match the value of the quantum $\chi(T)$ at its maximum, then the difference
in the slope at low and high temperatures is still clearly visible. One can
scale and shifts dimer’s $\chi(T)$ to match both the maximum and the high-
energy part of the curve (as was done in Ref. Grivel:CuC2O4:cryst ), but then
the slope in the low-energy part remains visibly different. This analysis
clearly shows that the infinite chain model reproduces not only the absolute
value of $\chi(T)$ but also the overall shape of the curve much better than
the scaled and shifted dimer model.

We also see that the $\chi(T)$ curve is broad and rather flat, which makes it
difficult to precisely extract $T_{max}$ from noisy experimental data. Hence,
the comparison on the basis of the maximum may be misleading, and one should
look at the overall shape. Near zero-temperature results would be important to
further highlight the discrepancy, but, unfortunately, there are no reliable
data, as the study that attempted such measurements clearly shows
contamination of the sampleGrivel:CuC2O4:cryst .

The right panel of Figure 3 shows the impact of different effects on the
overall magnetic susceptibility of the quantum spin chain model. The value of
the effective exchange constant from the 4-center fragment results in a
susceptibility curve that is far from the experimental one. Inclusion of
spin–orbit effects and weak correlation through perturbative triples improves
agreement. The impact of the basis set is relatively small but it improves
agreement even further. The more effects that are included, the closer the
agreement with the experiment becomes.

Inelastic neutron scattering allows one to probe the manifold of spin states
of a magnetic system. The intensity of inelastic neutron scattering is
proportional to the dynamic structure factor, given by

$\displaystyle
S^{\sigma\sigma^{\prime}}(k,\omega)=\sum_{i}\braket{\Psi_{GS}}{S_{k}^{\sigma}}{\Psi_{i}}\braket{\Psi_{i}}{S_{k}^{\sigma^{\prime}}}{\Psi_{GS}}|\delta(\omega-\omega_{i}),$
(11)

where $S^{\sigma}_{k}=\frac{1}{\sqrt{N}}\sum_{j}^{N}e^{iqj}S_{j}^{\sigma}$ are
the spin operators in the periodic representation, $\Psi_{GS}$ and $\Psi_{i}$
are the ground and excited states, and $\omega_{i}$ is the excitation energy
of excited state $\Psi_{i}$. Mixtures of weakly interacting antiferromagnetic
dimers would have a sharp intensity peak at the transfer energy matching the
singlet–triplet gapFurrer:INS:dimer:1977 ; Zheludev:dimer:INS:1996 . The
energy position of the peak does not depend on the momentum transfer $k$. In
contrast, extended spin chains show a characteristic dispersion law, where the
intensity is large at the distribution of energies and wave-
vectorsCaux:spin_chain:dyn_str_fact:2005 ; Kohno:INS:chain:2009 . Thus, future
neutron scattering experiments for copper oxalate could further confirm the
infinite spin-chain model and tighten the experimental error bar on the
effective exchange constant.

The rigorously derived ab initio effective Hamiltonian for copper oxalate
corresponds to the infinite Heisenberg spin chain model and validates the
nearest-neighbor approximation. This is the first purely theoretical
validation of the XXX Heisenberg model that is based on non-empirical all-
electron many-body calculations and not relying on spin-symmetry broken
solutions. Our results show that both strong and weak correlation contribute
to the effective exchange constant, whereas the dependence on the basis set is
relatively weak. The spin–orbit interaction introduces a noticeable additional
shift that further changes the effective exchange constant. Only after
accounting for all these effects, does the magnetic susceptibility of the
quantum spin model agree well with the experimental
valuesZelentsov:Cu:oxalate:1964 ; Dubicki:Cu:chains:1966 . This is the first
application of the EOM-SF-CC method to describe a periodic system, opening a
new route in the treatment of strongly correlated periodic systems. These
results provide a solid basis for future quantitative modeling of spin chains,
including magnetic field effects.

##  III Methods

Cartesian coordinates were extracted from the crystal unit cell parameters of
copper oxalate (CuC2O4) reported in Ref. 13. The crystal structure is composed
of the polymer CuC2O4 (Figure 2). The unit cell was repeated and the CuC2O4
stripes were extracted. Hydrogen atoms were added to the oxygen ends of the
chain with the corresponding to oxalic acid bond lengths and angles. The
Cartesian coordinates for all structures are given in the SI.

All calculations were performed with the Q-Chem quantum chemistry
packageqchem_2014 ; qchem_feature . We used Dunning’s cc-pVDZ and cc-pVTZ
basis sets for all atomsDunning:ccpvxz:1989 ; Duninng:ccpvxz:Sc-Zn ;
Dunning:ccpvxz:Ca ; Dunning:ccpvxz:Ga-Kr . Unrestricted orbitals were used in
calculations with open-shell references. To accelerate electronic structure
calculations, we used single precisionPavel:SP:2018  for CCSD, intermediates,
and EOM-SF-CCSD, combined with the libxm library for tensor
contractionsKaliman:LibXM:16 .

We used the open-shell frozen natural orbital truncation of the virtual
spacePavel:OSFNO:2019  with the truncation threshold corresponding to 99% of
the total preserved occupation in the virtual space. Core electrons were
frozen in all calculations. Triples corrections were computed using (fT)
perturbative triplessd3paper  for EOM-SF-CCSD. In the (fT) calculations for a
model dimer, the 40 lowest occupied orbitals were frozen for the dimer
fragment.

Calculations of spin–orbit couplingsPokhilko:SOC:19  were performed using EOM-
DIP-CCSD wave functionsWang:DIP:SOC:2020 .

Bloch’s and des Cloizeaux’s effective Hamiltonians were built for the open-
shell model spaces of the determinants expressed in the open-shell localized
orbitals following the procedure from Ref. Pokhilko:EffH:2020 . Localized
orbitals were determined using the Foster–Boys criterion of
localizationBoys:LMO:60 .

## Acknowledgments

We are grateful to Professor Devin Matthews from the South Methodist
University for his generous help with testing the effects of correlation
beyond perturbative triples correction. This work is supported by the
Department of Energy through the DE-SC0018910 grant.

P.P. current address is: Department of Chemistry, University of Michigan, Ann
Arbor, Michigan 48109, USA

The authors declare the following competing financial interest(s): A.I.K. is a
member of the Board of Directors and a part-owner of Q-Chem, Inc.

## References

  * (1) T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich,  Antiferromagnetic spintronics,  Nat. Nanotechnol. 11, 231 (2016).
  * (2) M. Mourigal, M. Enderle, A. Klöpperpieper, J.-S. Caux, A. Stunault, and H. M. Rønnow,  Fractional spinon excitations in the quantum Heisenberg antiferromagnetic chain,  Nat. Phys. 9, 435 (2013).
  * (3) P. Coleman and A. J. Schofield,  Quantum criticality,  Nature 433, 226 (2005).
  * (4) O. Breunig, M. Garst, A. Klümper, J. Rohrkamp, M. M. Turnbull, and T. Lorenz,  Quantum criticality in the spin-1/2 Heisenberg chain system copper pyrazine dinitrate,  Science 3 (2017).
  * (5) B. Bleaney and K. D. Bowers,  Anomalous paramagnetism of copper acetate,  Proc. R. Soc. London A 214, 451 (1952).
  * (6) R. Maurice, K. Sivalingam, D. Ganyushin, N. Guihéry, C. de Graaf, and F. Neese,  Theoretical determination of the zero-field splitting in copper acetate monohydrate,  Inorg. Chem. 50, 6229 (2011).
  * (7) H. Bethe,  Zur Theorie der Metalle,  Z. Phys. 71, 205 (1931).
  * (8) N. Motoyama, H. Eisaki, and S. Uchida,  Magnetic susceptibility of ideal spin $1/2$ Heisenberg antiferromagnetic chain systems, Sr2CuO3 and SrCuO2,  Phys. Rev. Lett. 76, 3212 (1996).
  * (9) D. C. Johnston, J. W. Johnson, D. P. Goshorn, and A. J. Jacobson,  Magnetic susceptibility of (VO)2P2O7: A one-dimensional spin-1/2 Heisenberg antiferromagnet with a ladder spin configuration and a singlet ground state,  Phys. Rev. B 35, 219 (1987).
  * (10) P. M. Richards, R. K. Quinn, and B. Morosin,  Magnetic and structural properties of manganese pyridine linear chain salts,  J. Chem. Phys. 59, 4474 (1973).
  * (11) B. N. Figgis and D. J. Martin,  The magnetic properties and structure of the cupric $\alpha$,$\omega$-dicarboxylates,  Inorg. Chem. 5, 100 (1966).
  * (12) L. Dubicki, C. M. Harris, E. Kokot, and R. L. Martin,  Magnetic studies with copper (II) salts. VII. The structure of copper (II) $\alpha$, $\omega$-dicarboxylates and their amine derivatives,  Inorg. Chem. 5, 93 (1966).
  * (13) A. N. Christensen, B. Lebech, N. H. Andersen, and J.-C. Grivel,  The crystal structure of paramagnetic copper(II) oxalate (CuC2O4): formation and thermal decomposition of randomly stacked anisotropic nano-sized crystallites,  Dalton Trans. 43, 16754 (2014).
  * (14) T. Kenneth K. T. McGregor and Z. G. Soos,  Anisotropic exchange and temperature-dependent electron paramagnetic resonance line width in one-dimensional copper(II) complexes. 2. Magnetic properties of copper(II) oxalate-1/3-water,  Inorg. Chem. 15, 2159 (1976).
  * (15) A. Michalowicz, J. J. Girerd, and J. Goulon,  EXAFS [extended x-ray absorption fine structure] determination of the copper oxalate structure. Relation between structure and magnetic properties,  Inorg. Chem. 18, 3004 (1979).
  * (16) A. I. Krylov,  Equation-of-motion coupled-cluster methods for open-shell and electronically excited species: The hitchhiker’s guide to Fock space,  Annu. Rev. Phys. Chem. 59, 433 (2008).
  * (17) J. P. Malrieu, R. Caballol, C. J. Calzado, C. de Graaf, and N. Guihéry,  Magnetic interactions in molecules and highly correlated materials: Physical content, analytical derivation, and rigorous extraction of magnetic Hamiltonians,  Chem. Rev. 114, 429 (2013).
  * (18) P. Pokhilko and A. I. Krylov,  Effective Hamiltonians derived from equation-of-motion coupled-cluster wave-functions: Theory and application to the Hubbard and Heisenberg Hamiltonians,  J. Chem. Phys. 152, 094108 (2020).
  * (19) N. J. Mayhall and M. Head-Gordon,  Computational quantum chemistry for multiple-site Heisenberg spin couplings made simple: Still only one spin-flip required,  J. Phys. Chem. Lett. 6, 1982 (2015).
  * (20) A. I. Krylov,  Size-consistent wave functions for bond-breaking: The equation-of-motion spin-flip model,  Chem. Phys. Lett. 338, 375 (2001).
  * (21) N. Orms and A. I. Krylov,  Singlet-triplet energy gaps and the degree of diradical character in binuclear copper molecular magnets characterized by spin-flip density functional theory,  Phys. Chem. Chem. Phys. 20, 13127 (2018).
  * (22) P. Pokhilko, D. Izmodenov, and A. I. Krylov,  Extension of frozen natural orbital approximation to open-shell references: Theory, implementation, and application to single-molecule magnets,  J. Chem. Phys. 152, 034105 (2020).
  * (23) M. A. El-Sayed,  Spin–orbit coupling and the radiationless processes in nitrogen heterocycles,  J. Chem. Phys. 38, 2834 (1963).
  * (24) M. A. El-Sayed,  Triplet state: Its radiative and non-radiative properties,  Acc. Chem. Res. 1, 8 (1968).
  * (25) P. Pokhilko and A. I. Krylov,  Quantitative El-Sayed rules for many-body wavefunctions from spinless transition density matrices,  J. Phys. Chem. Lett. 10, 4857 (2019).
  * (26) M. E. Fisher,  Magnetism in one-dimensional systems—the Heisenberg model for infinite spin,  Am. J. Phys. 32, 343 (1964).
  * (27) D. C. Johnston, R. K. Kremer, M. Troyer, X. Wang, A. Klümper, S. L. Bud’ko, A. F. Panchula, and P. C. Canfield,  Thermodynamics of spin ${S}=1/2$ antiferromagnetic uniform and alternating-exchange Heisenberg chains,  Phys. Rev. B 61, 9558 (2000).
  * (28) V. V. Zelentsov and T. G. Aminov,  Magnetic susceptibility of cupric oxalate and cupric succinate,  Dokl. Akad. Nauk SSSR 158, 1393 (1964).
  * (29) T. Xiang,  Thermodynamics of quantum Heisenberg spin chains,  Phys. Rev. B 58, 9142 (1998).
  * (30) A. Furrer and H. U. Güdel,  Interference effects in neutron scattering from magnetic clusters,  Phys. Rev. Lett. 39, 657 (1977).
  * (31) A. Zheludev, G. Shirane, Y. Sasago, M. Hase, and K. Uchinokura,  Dimerized ground state and magnetic excitations in CaCuGe2O6,  Phys. Rev. B 53, 11642 (1996).
  * (32) J.-S. Caux, R. Hagemans, and J. M. Maillet,  Computation of dynamical correlation functions of Heisenberg chains: the gapless anisotropic regime,  J. Stat. Mech. Theory Exp 2005, P09003 (2005).
  * (33) M. Kohno,  Dynamically dominant excitations of string solutions in the spin-$1/2$ antiferromagnetic Heisenberg chain in a magnetic field,  Phys. Rev. Lett. 102, 037203 (2009).
  * (34) Y. Shao, Z. Gan, E. Epifanovsky, A.T.B. Gilbert, M. Wormit, J. Kussmann, A.W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P.R. Horn, L.D. Jacobson, I. Kaliman, R.Z. Khaliullin, T. Kus, A. Landau, J. Liu, E.I. Proynov, Y.M. Rhee, R.M. Richard, M.A. Rohrdanz, R.P. Steele, E.J. Sundstrom, H.L. Woodcock III, P.M. Zimmerman, D. Zuev, B. Albrecht, E. Alguires, B. Austin, G.J.O. Beran, Y.A. Bernard, E. Berquist, K. Brandhorst, K.B. Bravaya, S.T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S.H. Chien, K.D. Closser, D.L. Crittenden, M. Diedenhofen, R.A. DiStasio Jr., H. Do, A.D. Dutoi, R.G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M.W.D. Hanson-Heine, P.H.P. Harbach, A.W. Hauser, E.G. Hohenstein, Z.C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R.A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C.M. Krauter, K.U. Laog, A. Laurent, K.V. Lawler, S.V. Levchenko, C.Y. Lin, F. Liu, E. Livshits, R.C. Lochan, A. Luenser, P. Manohar, S.F. Manzer, S.-P. Mao, N. Mardirossian, A.V. Marenich, S.A. Maurer, N.J. Mayhall, C.M. Oana, R. Olivares-Amaya, D.P. O’Neill, J.A. Parkhill, T.M. Perrine, R. Peverati, P.A. Pieniazek, A. Prociuk, D.R. Rehn, E. Rosta, N.J. Russ, N. Sergueev, S.M. Sharada, S. Sharmaa, D.W. Small, A. Sodt, T. Stein, D. Stuck, Y.-C. Su, A.J.W. Thom, T. Tsuchimochi, L. Vogt, O. Vydrov, T. Wang, M.A. Watson, J. Wenzel, A. White, C.F. Williams, V. Vanovschi, S. Yeganeh, S.R. Yost, Z.-Q. You, I.Y. Zhang, X. Zhang, Y. Zhou, B.R. Brooks, G.K.L. Chan, D.M. Chipman, C.J. Cramer, W.A. Goddard III, M.S. Gordon, W.J. Hehre, A. Klamt, H.F. Schaefer III, M.W. Schmidt, C.D. Sherrill, D.G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A.T. Bell, N.A. Besley, J.-D. Chai, A. Dreuw, B.D. Dunietz, T.R. Furlani, S.R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D.S. Lambrecht, W.Z. Liang, C. Ochsenfeld, V.A. Rassolov, L.V. Slipchenko, J.E. Subotnik, T. Van Voorhis, J.M. Herbert, A.I. Krylov, P.M.W. Gill, and M. Head-Gordon,  Advances in molecular quantum chemistry contained in the Q-Chem 4 program package,  Mol. Phys. 113, 184 (2015).
  * (35) A. I. Krylov and P. M. W. Gill,  Q-Chem: An engine for innovation,  WIREs: Comput. Mol. Sci. 3, 317 (2013).
  * (36) T. H. Dunning, Jr.,  Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen,  J. Chem. Phys. 90, 1007 (1989).
  * (37) N. B. Balabanov and K. A. Peterson,  Systematically convergent basis sets for transition metals. I. All-electron correlation consistent basis sets for the 3d elements Sc-Zn,  J. Chem. Phys. 123, 064107 (2005).
  * (38) J. Koput and K. A. Peterson,  Ab initio potential energy surface and vibrational–rotational energy levels of X${}^{2}\Sigma$+ CaOH,  J. Phys. Chem. A 106, 9595 (2002).
  * (39) A. K. Wilson, D. E. Woon, K. A. Peterson, and T. H. Dunning,  Gaussian basis sets for use in correlated molecular calculations. ix. the atoms gallium through krypton,  J. Chem. Phys. 110, 7667 (1999).
  * (40) P. Pokhilko, E. Epifanovskii, and A. I. Krylov,  Double precision is not needed for many-body calculations: Emergent conventional wisdom,  J. Chem. Theory Comput. 14, 4088 (2018).
  * (41) I. Kaliman and A. I. Krylov,  New algorithm for tensor contractions on multi-core CPUs, GPUs, and accelerators enables CCSD and EOM-CCSD calculations with over 1000 basis functions on a single compute node,  J. Comput. Chem. 38, 842 (2017).
  * (42) P. U. Manohar and A. I. Krylov,  A non-iterative perturbative triples correction for the spin-flipping and spin-conserving equation-of-motion coupled-cluster methods with single and double substitutions,  J. Chem. Phys. 129, 194105 (2008).
  * (43) P. Pokhilko, E. Epifanovsky, and A. I. Krylov,  General framework for calculating spin–orbit couplings using spinless one-particle density matrices: theory and application to the equation-of-motion coupled-cluster wave functions,  J. Chem. Phys. 151, 034106 (2019).
  * (44) H. Zhao, Z. Wang, M. Guo, and F. Wang,  Splittings of d8 configurations of late-transition metals with eom-dip-ccsd and fsccsd methods,  J. Chem. Phys. 152, 134105 (2020).
  * (45) S. F. Boys,  Construction of some molecular orbitals to be approximately invariant for changes from one molecule to another,  Rev. Mod. Phys. 32, 296 (1960).