# CoTAR: Chain-of-Thought Attribution Reasoning with Multi-level Granularity

Moshe Berchansky Daniel Fleischer Moshe Wasserblat Peter Izsak
Intel Labs
<EMAIL_ADDRESS>

###### Abstract

State-of-the-art performance in QA tasks is currently achieved by systems
employing Large Language Models (LLMs), however these models tend to
hallucinate information in their responses. One approach focuses on enhancing
the generation process by incorporating attribution from the given input to
the output. However, the challenge of identifying appropriate attributions and
verifying their accuracy against a source is a complex task that requires
significant improvements in assessing such systems. We introduce an
attribution-oriented Chain-of-Thought reasoning method to enhance the accuracy
of attributions. This approach focuses the reasoning process on generating an
attribution-centric output. Evaluations on two context-enhanced question-
answering datasets using GPT-4 demonstrate improved accuracy and correctness
of attributions. In addition, the combination of our method with finetuning
enhances the response and attribution accuracy of two smaller LLMs, showing
their potential to outperform GPT-4 in some cases.

CoTAR: Chain-of-Thought Attribution Reasoning with Multi-level Granularity

Moshe Berchansky Daniel Fleischer Moshe Wasserblat Peter Izsak Intel Labs
<EMAIL_ADDRESS>

##  1 Introduction

Text generation from sources, such as information obtained through retrieval,
is a key aspect of grounded question answering. One way to ensure the
credibility of these models is through attributed text generation. This method
pairs the generated text with supporting evidence, enhancing the model’s
trustworthiness and allowing for easier detection of errors.

The pairing process can be done in several levels. Firstly, we can pair each
sentence in our response with a set of passage identifiers Gao et al. (2023).
On a model detailed level, individual spans in the answer can correspond to a
specific passage, hence attributing specific portions to various passages
Schuster et al. (2023). These spans are designed to be direct copies from the
passages, assuring that the content provided is fully supported. Since not all
the spans are copied, the model has a sufficient degree of flexibility to
produce both a coherent and factually correct answer.

Figure 1: Usage of CoT for attribution-based answers.

However, these approaches fail at times in two primary ways to cite
accurately. On one hand, the model might focus only a very specific portion of
the input, therefore missing relevant sections when it constructs the answer.
On the other hand, it might cite too many passages, some of which might not be
relevant to the answer.

To solve these, we propose utilizing a Chain-of-Thought approach (CoT) Wei et
al. (2022b), which we denote as CoTAR, which allows the model to perform
reasoning over the input passages before generating the output. In our
approach, we instruct the model to extract relevant information from the
passages as well as specify the form of attribution to generate. The
information can be in the span, sentence, or passage level. We inspect how the
various CoT methods affect the model’s output, in each of the citation levels.
In addition to few-shot instruction, we also finetune various models using our
method, resulting in models that are competitive with and even outperform
GPT-4 in some cases.

We can summarize our contributions in this paper as follows:

  * •

We perform rigorous measurements of both the answer quality and citation
quality, across multiple models and citation levels.

  * •

We show that utilizing CoT reasoning improves the ability of an LLM to produce
both better quality answers, and more precise and faithful citations from the
source, demonstrated on multiple models.

  * •

We demonstrate that using by finetuning, smaller models can be competitive
with or outperform GPT-4 in some cases in answer and citation quality metrics.

##  2 Method

###  2.1 Attribution-Oriented Question Answering

Attribution-oriented question answering can be defined as a task where, given
a query and a set of relevant contextual resources, the objective is to
accurately answer the question while attributing specific portions or the
entire answer to the appropriate contextual sources.

We can identify three levels of attribution with different granularity levels:
Span, Sentence and Passage. Attribution levels refer to how citations are
displayed within an answer, indicating whether complete passages, individual
sentences, or specific text segments are credited back to their original
sources. Detailed examples for each level are available in Table 4 in the
appendix.

###  2.2 Chain-of-Thought Attribution Reasoning

We propose a multi-step CoT Wei et al. (2022b) reasoning scheme with varying
levels of attribution, similar to the three levels presented in Section 2,
hypothesizing that this could encourage the model to generate more accurate
answers. The process involves identifying the most crucial aspects of the
given context for answering the question, by incorporating direct citations to
the referenced parts. We denote the process of granular attribution as CoT
attribution guidance, or CoT method for short. The three levels are as
follows:

  * •

Span Guidance: Produce the relevant spans of information per passage.

  * •

Sentence Guidance: For every passage, write sentences that summarize how the
passage answers the question.

  * •

Passage Guidance: State which passages are relevant for the question.

We explore all combinations of CoT methods with citation levels, and asses how
each one effects the generated answers. A example of all three levels of
attribution-guidance is in Table 5 in the appendix.

Level CoT Answer Quality Citation Quality BERT HEM RL Avg. ALCE F1 CSCA DOC F1
SEM-F1 Avg. Span None 91.71 67.58 66.08 75.12 77.06 81.15 85.33 67.57 77.78
Span 92.06 66.95 68.55 75.85 78.83 89.88 86.54 72.23 81.87 Sent. 91.87 66.88
68.76 75.84 84.76 69.69 86.38 70.52 77.84 Pass. 91.91 66.89 67.67 75.49 79.21
86.94 86.49 69.87 80.63 Sent. None 91.47 62.48 62.18 72.04 82.79 - 85.32 59.20
75.77 Span 91.92 62.06 64.76 72.91 82.39 - 86.71 61.74 76.94 Sent. 91.96 64.59
67.36 74.64 86.42 - 86.51 64.44 79.12 Pass. 91.69 61.19 64.42 72.43 84.25 -
85.92 61.79 77.32 Pass. None 91.55 61.13 62.76 71.81 75.38 - 85.19 44.31 68.29
Span 91.90 61.66 64.70 72.75 71.60 - 86.17 43.41 67.06 Sent. 91.84 62.43 65.32
73.20 70.63 - 86.41 42.59 66.54 Pass. 91.69 60.40 63.88 71.99 73.28 - 85.64
44.12 67.68

Table 1: Results of GPT-4 on the QuoteSUM dataset. The results are categorized
by response attribution (level), and for each level, every CoT approach is
demonstrated. The best results per level are marked in bold.

##  3 Evaluation Metrics

In order to evaluate the answers, we specify metrics of two types. Answer
Quality, measuring similarity between the predicted answer and the gold
answer. Citation Quality, assessing the similarity between the cited text, the
cited passage, and the citations present in the gold answer.

For answer quality , we use the n-grams based method ROUGE-L (RL), and the
semantic method BERTScore (BERT) Zhang* et al. (2020). For hallucination
evaluation, we use the HEM111vectara/hallucination_evaluation_model model,
which was finetuned on NLI datasets, and predict whether two texts are
factually consistent. For citation quality, we measure the quality of the
cited content; we propose specific metrics for each of the citation levels. We
note that span citations must match the source passages precisely, while
sentence and passage attributions are not required to do so. We use the SEM-F1
metric proposed by Schuster et al. (2023) for n-gram token level similarity
between the cited content. In addition, we include the Citation
Precision/Recall introduced in Gao et al. (2023), and combine them into ALCE
F1, by using the harmonic mean. We also measure the F1 score over passage
indices between the cited passages and the passages cited by the expected
answer, denoted as DOC F1, and the Correct Span Citation Attribution (CSCA)
which indicates whether a predicted span is a direct span from the attributed
passage.

##  4 Experimental Setup

###  4.1 Datasets

We use the QuoteSum (QSUM) Schuster et al. (2023) dataset, which contains
semi-extractive answers written by humans, for natural questions and generated
ones. In addition, we include MS MARCO
(MS)222https://huggingface.co/datasets/ms_marco Campos et al. (2016), which is
a crowd sourced dataset of responses to queries from Bing. Each example
include 10 passages which can be relevant or irrelevant to answering the
question.

###  4.2 Span Dataset Preparation

In order to include span citations in datasets such as MSMARCO, we can create
span-based citation labels out of passage-based labels. We do so by finding
the longest common substrings between the input passages and the answer, and
choosing those that contain named entities, thus ensuring the semantic
relevance of the cited spans. The answer is then formatted as proposed by
Schuster et al. (2023). For more details about the process, we refer the
reader to Algorithm 1.

###  4.3 Models

We use GPT-4 Achiam et al. (2023) for evaluating our CoT method. We use the
Mistral 7B Jiang et al. (2023) decoder model and the Flan-T5 XXL Chung et al.
(2022) encoder-decoder as representative small-sized models in our further
finetuning evaluations.

###  4.4 Method Comparison

For each of the tasks described in Section 2.1, we utilize every CoT method
described in Section 2.2. We run each combination with GPT-4, a small decoder-
only model, and an encoder-decoder model. All results are presented on the
test sets of the datasets mentioned in Section 4.1. For full implementation
details of the training and inference, we refer the reader to section A in the
Appendix.

##  5 Experimental Results

To examine how CoT affects the model, we inspect the GPT-4 case, and provide a
detailed overview of the metrics in Table 1. In the answer metrics, the
differences between the various CoT methods and the standard run (shown as
"None") are not significant, aside from the sentence level, where the sentence
CoT method produces better scores across all metrics. For the citation quality
metrics, the span and sentence level produce better scores than their
counterparts when using the appropriate CoT methods. In particular, the CSCA
metric in the span level case showcases a significant advantage, mainly since
the model is guided to produce the spans relevant for the answer in advance,
thus allowing the model to accurately reference them. At the passage level,
the SEM-F1 and ALCE-F1 scores showcase an advantage for the standard run. We
hypothesize that the addition of CoT for such a broad citation task is less
helpful than in the more granular sentence and span level citations. We
observe that the GPT-4 model performed best in the span level when using the
span CoT method, and the same applies for the sentence level with sentence
CoT. In the passage level, while the CoT methods improve the answer quality,
they do not aid in the citation quality.

We showcase similar results for all the models, across the various citation
levels and CoT methods in Tables 7 and 8. For each setting, we display all
metrics and category averages for answer and citation quality, respectively.

We present an overview of the various combinations of citation levels and CoT
methods in Table 2. For smaller decoder-only models (i.e. Mistral 7B), the
answer scores perform best in general when the citation level is used with its
appropriate CoT method. In the citation metrics, results are different between
the QSUM and MS. While the CoT methods behave well in the MS dataset, they are
proving less helpful in QSUM, with the main difference being a higher ALCE-F1
score for the standard run. The difference can be mostly attributed to the
higher complexity of the QSUM dataset, as opposed to the shorter and simpler
MS answers.

For encoder-decoder models, the CoT methods produce better results than the
standard run, both in the answer and citation metrics. However, aside from the
MS citation quality results, the advantage of using each citation level task
with its corresponding CoT method is less prominent. The sentence CoT method
is highly dominant in most of the metrics in the QSUM dataset, and passage CoT
method being prominent in the MS dataset.

Finally, we observe in Table 2 that the smaller models are competitive with
and even outperform GPT-4. This conclusion is similar to the one reached by
Schuster et al. (2023), when comparing to other LLMs without finetuning.

| CoT | Answer Quality | Citation Quality
---|---|---|---
Span | Sent. | Pass. | Span | Sent. | Pass.
| QuoteSUM
GPT-4 | None | 75.1 | 72.0 | 71.8 | 77.8 | 75.8 | 68.3
Span | 75.9 | 72.9 | 72.8 | 81.9 | 76.9 | 67.1
Sent. | 75.8 | 74.6 | 73.2 | 77.8 | 79.1 | 66.5
Pass. | 75.5 | 72.4 | 72.0 | 80.6 | 77.3 | 67.7
Mistral 7B | None | 74.2 | 72.6 | 74.4 | 83.0 | 72.7 | 66.4
Span | 76.2 | 75.6 | 75.1 | 82.2 | 72.3 | 64.7
Sent. | 75.2 | 75.6 | 75.1 | 79.9 | 73.4 | 65.0
Pass. | 75.6 | 74.8 | 75.8 | 82.0 | 73.5 | 66.1
Flan-T5 | None | 76.6 | 73.8 | 74.0 | 84.7 | 71.4 | 64.3
Span | 77.3 | 71.8 | 70.6 | 84.1 | 69.7 | 62.1
Sent. | 78.4 | 78.2 | 76.7 | 84.7 | 79.7 | 66.5
Pass. | 78.5 | 74.4 | 73.1 | 85.6 | 72.7 | 61.2
| MSMARCO
Mistral 7B | None | 69.6 | 69.1 | 70.2 | 68.4 | 66.1 | 66.0
Span | 70.8 | 70.5 | 71.5 | 74.8 | 67.8 | 68.9
Sent. | 71.0 | 71.6 | 71.3 | 72.4 | 69.1 | 68.4
Pass. | 70.8 | 71.6 | 71.6 | 74.4 | 69.2 | 69.1
Flan-T5 | None | 70.4 | 71.5 | 72.4 | 74.3 | 67.7 | 68.5
Span | 70.4 | 70.3 | 68.4 | 74.7 | 68.7 | 68.6
Sent. | 71.3 | 71.6 | 71.9 | 73.2 | 69.0 | 69.3
Pass. | 69.0 | 72.7 | 72.6 | 73.9 | 69.0 | 70.0
Table 2: Average answer and citation quality results, per citation level, per
CoT method, for QuoteSUM and MSMARCO. Rows: Citation levels. Columns: CoT
method used and metrics averages per level.

##  6 Related Work

Several recent studies have suggested different techniques for text generation
including citation attribution, with varying degrees of detail. Attributing
relevant information to web pages was done in Thoppilan et al. (2022); Liu et
al. (2023). Menick et al. (2022) attributes the response to specific snippets
from the source passages. In particular, Gao et al. (2023) assigns attribution
for each produced sentence to one or several input passages. Our sentence
level citation approach and the evaluation metrics are based on this work.

Substantial work has been done to reduce hallucination in the generated text
produced by the model Li et al. (2024); Rawte et al. (2023), as citation has
the ability to increase the reliability of the generated text. Some work has
been done on correlating answer to text sections in the retrieved passages Gao
et al. (2023); Bohnet et al. (2023). However, these techniques require
additional processing.

The usage of CoT has shown significant improvements in a wide variety of tasks
Wei et al. (2022b). Additionally, in-context learning can be added to the
input, by providing examples of how to solve the task at hand Wei et al.
(2022a). We explore the effects of using CoT with the various citation level
tasks, in addition to the usage of in-context learning examples in every
scenario. Luo et al. (2023) instruct and finetune models to produce claims,
based on retrieved passages, differentiating between informative and
distracting input segments. We specify the text granularity level we expect
the model should cite and phrase its answers. Slobodkin et al. (2024) focus on
the sentence citation task, splitting up the generation task into multiple
reasoning steps, including CoT as well. Our approach focuses on using a single
call, producing the answer at-once, with a CoT approach.

##  7 Conclusions

Through extensive measurements applied across various models and citation
levels, we have been able to ascertain the quality of both answers and
citations. Our findings indicate that the use of our CoTAR reasoning
significantly enhances the capacity of a model to generate superior quality
answers and more accurate, faithful citations from the source. This
improvement is evident in both decoder-only and encoder-decoder models.
Furthermore, our research shows that by using finetuning, smaller models can
compete with or even surpass the performance of GPT-4 across a wide range of
answer and citation metrics. This is particularly evident when considering the
diversity of citation levels and CoT methods.

##  8 Limitations

In our experimentation, we did not compare our approach directly to the work
done in Slobodkin et al. (2024), due to lack of access to the publicly
available code-base. In our GPT-4 runs, we ran the experimentation solely on
the QuoteSUM dataset, in order to compare GPT-4 to our models on an previously
established dataset. When considering datasets for our work, we primarily
focused on the QuoteSUM dataset from Schuster et al. (2023), and also the
MSMARCO dataset, due to the extensive annotation of the answers and the
accompanying passages. For future work, we will add similar datasets in
structure and content.

##  9 Potential Risks

Since this project is primarily focused on assessing current models using
datasets that are openly accessible, we do not foresee any possible adverse
effects.

## References

  * Achiam et al. (2023) Josh Achiam, Steven Adler, Sandhini Agarwal, Lama Ahmad, Ilge Akkaya, Florencia Leoni Aleman, Diogo Almeida, Janko Altenschmidt, Sam Altman, Shyamal Anadkat, et al. 2023.  Gpt-4 technical report.  _arXiv preprint arXiv:2303.08774_.
  * Bohnet et al. (2023) Bernd Bohnet, Vinh Q. Tran, Pat Verga, Roee Aharoni, Daniel Andor, Livio Baldini Soares, Massimiliano Ciaramita, Jacob Eisenstein, Kuzman Ganchev, Jonathan Herzig, Kai Hui, Tom Kwiatkowski, Ji Ma, Jianmo Ni, Lierni Sestorain Saralegui, Tal Schuster, William W. Cohen, Michael Collins, Dipanjan Das, Donald Metzler, Slav Petrov, and Kellie Webster. 2023.  Attributed question answering: Evaluation and modeling for attributed large language models.  _Preprint_ , arXiv:2212.08037.
  * Campos et al. (2016) Daniel Fernando Campos, Tri Nguyen, Mir Rosenberg, Xia Song, Jianfeng Gao, Saurabh Tiwary, Rangan Majumder, Li Deng, and Bhaskar Mitra. 2016.  Ms marco: A human generated machine reading comprehension dataset.  _ArXiv_ , abs/1611.09268.
  * Chung et al. (2022) Hyung Won Chung, Le Hou, Shayne Longpre, Barret Zoph, Yi Tay, William Fedus, Yunxuan Li, Xuezhi Wang, Mostafa Dehghani, Siddhartha Brahma, Albert Webson, Shixiang Shane Gu, Zhuyun Dai, Mirac Suzgun, Xinyun Chen, Aakanksha Chowdhery, Alex Castro-Ros, Marie Pellat, Kevin Robinson, Dasha Valter, Sharan Narang, Gaurav Mishra, Adams Yu, Vincent Zhao, Yanping Huang, Andrew Dai, Hongkun Yu, Slav Petrov, Ed H. Chi, Jeff Dean, Jacob Devlin, Adam Roberts, Denny Zhou, Quoc V. Le, and Jason Wei. 2022.  Scaling instruction-finetuned language models.  _Preprint_ , arXiv:2210.11416.
  * Gao et al. (2023) Tianyu Gao, Howard Yen, Jiatong Yu, and Danqi Chen. 2023.  Enabling large language models to generate text with citations.  _Preprint_ , arXiv:2305.14627.
  * Honnibal and Montani (2017) Matthew Honnibal and Ines Montani. 2017.  spaCy 2: Natural language understanding with Bloom embeddings, convolutional neural networks and incremental parsing.  To appear.
  * Jiang et al. (2023) Albert Q. Jiang, Alexandre Sablayrolles, Arthur Mensch, Chris Bamford, Devendra Singh Chaplot, Diego de las Casas, Florian Bressand, Gianna Lengyel, Guillaume Lample, Lucile Saulnier, Lélio Renard Lavaud, Marie-Anne Lachaux, Pierre Stock, Teven Le Scao, Thibaut Lavril, Thomas Wang, Timothée Lacroix, and William El Sayed. 2023.  Mistral 7b.  _Preprint_ , arXiv:2310.06825.
  * Li et al. (2024) Junyi Li, Jie Chen, Ruiyang Ren, Xiaoxue Cheng, Wayne Xin Zhao, Jian-Yun Nie, and Ji-Rong Wen. 2024.  The dawn after the dark: An empirical study on factuality hallucination in large language models.  _Preprint_ , arXiv:2401.03205.
  * Liu et al. (2023) Nelson F. Liu, Tianyi Zhang, and Percy Liang. 2023.  Evaluating Verifiability in Generative Search Engines.  _arXiv preprint_.  ArXiv: 2304.09848 [cs].
  * Luo et al. (2023) Hongyin Luo, Yung-Sung Chuang, Yuan Gong, Tianhua Zhang, Yoon Kim, Xixin Wu, Danny Fox, Helen Meng, and James Glass. 2023.  Sail: Search-augmented instruction learning.  _Preprint_ , arXiv:2305.15225.
  * Menick et al. (2022) Jacob Menick, Maja Trebacz, Vladimir Mikulik, John Aslanides, Francis Song, Martin Chadwick, Mia Glaese, Susannah Young, Lucy Campbell-Gillingham, Geoffrey Irving, and Nat McAleese. 2022.  Teaching language models to support answers with verified quotes.  _Preprint_ , arXiv:2203.11147.
  * Rawte et al. (2023) Vipula Rawte, Prachi Priya, S. M Towhidul Islam Tonmoy, S M Mehedi Zaman, Amit Sheth, and Amitava Das. 2023.  Exploring the relationship between llm hallucinations and prompt linguistic nuances: Readability, formality, and concreteness.  _Preprint_ , arXiv:2309.11064.
  * Schuster et al. (2023) Tal Schuster, Ádám Dániel Lelkes, Haitian Sun, Jai Gupta, Jonathan Berant, William W. Cohen, and Donald Metzler. 2023.  Semqa: Semi-extractive multi-source question answering.  _ArXiv_ , abs/2311.04886.
  * Slobodkin et al. (2024) Aviv Slobodkin, Eran Hirsch, Arie Cattan, Tal Schuster, and Ido Dagan. 2024.  Attribute first, then generate: Locally-attributable grounded text generation.  _Preprint_ , arXiv:2403.17104.
  * Thoppilan et al. (2022) Romal Thoppilan, Daniel De Freitas, Jamie Hall, Noam Shazeer, Apoorv Kulshreshtha, Heng-Tze Cheng, Alicia Jin, Taylor Bos, Leslie Baker, Yu Du, YaGuang Li, Hongrae Lee, Huaixiu Steven Zheng, Amin Ghafouri, Marcelo Menegali, Yanping Huang, Maxim Krikun, Dmitry Lepikhin, James Qin, Dehao Chen, Yuanzhong Xu, Zhifeng Chen, Adam Roberts, Maarten Bosma, Vincent Zhao, Yanqi Zhou, Chung-Ching Chang, Igor Krivokon, Will Rusch, Marc Pickett, Pranesh Srinivasan, Laichee Man, Kathleen Meier-Hellstern, Meredith Ringel Morris, Tulsee Doshi, Renelito Delos Santos, Toju Duke, Johnny Soraker, Ben Zevenbergen, Vinodkumar Prabhakaran, Mark Diaz, Ben Hutchinson, Kristen Olson, Alejandra Molina, Erin Hoffman-John, Josh Lee, Lora Aroyo, Ravi Rajakumar, Alena Butryna, Matthew Lamm, Viktoriya Kuzmina, Joe Fenton, Aaron Cohen, Rachel Bernstein, Ray Kurzweil, Blaise Aguera-Arcas, Claire Cui, Marian Croak, Ed Chi, and Quoc Le. 2022.  Lamda: Language models for dialog applications.  _Preprint_ , arXiv:2201.08239.
  * Wei et al. (2022a) Jason Wei, Yi Tay, Rishi Bommasani, Colin Raffel, Barret Zoph, Sebastian Borgeaud, Dani Yogatama, Maarten Bosma, Denny Zhou, Donald Metzler, Ed H. Chi, Tatsunori Hashimoto, Oriol Vinyals, Percy Liang, Jeff Dean, and William Fedus. 2022a.  Emergent abilities of large language models.  _Preprint_ , arXiv:2206.07682.
  * Wei et al. (2022b) Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, Ed Huai hsin Chi, F. Xia, Quoc Le, and Denny Zhou. 2022b.  Chain of thought prompting elicits reasoning in large language models.  _ArXiv_ , abs/2201.11903.
  * Zhang* et al. (2020) Tianyi Zhang*, Varsha Kishore*, Felix Wu*, Kilian Q. Weinberger, and Yoav Artzi. 2020.  Bertscore: Evaluating text generation with bert.  In _International Conference on Learning Representations_.

##  Appendix A Implementation Details

The prompt is built using 4 fewshot examples, each composed of a question, a
list of passages, and the expected answer in the task’s format. As done by
Schuster et al. (2023), the prompt for the models is the same as they have
specified, with an example in Table 6, and the fewshot examples are chosen
with the question similarity between the test question and the training
example questions. As stated above, the expected answers cite their passages
of origin in one of the levels described in Table 4. For the CoT setting, a
prefix is attached in one of the methods showcased in Table 5. Text generation
is done without sampling, with a maximum length of 2k tokens. For QuoteSUM,
the train and test sets are identical to the ones used by Schuster et al.
(2023). For MSMARCO, we utilize 2750 samples from the created dataset for
training, and a 1000 for the test set. We specify the training parameters in
Table 3.

##  Appendix B Models

The weights for the models we have used have all been retrieved from
HuggingFace, aside from GPT-4, for which we use Azure OpenAI:

  * •

Mistral 7B: mistralai/Mistral-7B-Instruct-v0.2, distributed under the Apache
License 2.0.

  * •

Flan-T5: google/flan-t5-xxl, distributed under the Apache License 2.0.

  * •

GPT-4: API access using Azure OpenAI, model `gpt-4-32k-0613`.

##  Appendix C Additional Results

Parameter | Decoder-Only | Encoder-Decoder
---|---|---
Max Target Len. | - | 4096
Max Seq. Len. | 8192 | 32768
Lora R | 16 | 16
Lora $\alpha$ | 32 | 32
Lora Dropout | 0.05 | 0.05
Lora Bias | None | None
Lora Modules | gate, down, up, $q$, $v$, $k$, $o$ | $q,v$
LR | 5e-5 | 2e-4
LR Scheduler | Linear | Linear
Weight Decay | 0 | 0.01
Precision | bfloat16 | bfloat16
Batch Size | 1 | 1
Epochs | 3 | 3
Warmup Ratio | 0 | 0.1
Table 3: The training parameters used for Lora finetuning for the decoder-only
and encoder-decoder models. Algorithm 1 Generating span attirubtions from a
passage-labeled dataset.

1:Input: Answer $a$, passages $\\{p_{i}\\}$

2:Extract named entities $\text{NER}(a)$ using spaCy Honnibal and Montani
(2017).

3:Find common substrings $CS(a,p_{i})$ between every passage $p_{i}$ and
answer $a$.

4:Each substring $sc\in CS(a,p_{i})$: contains $sc_{t}$ (string text) and
passage index $sc_{p}$.

5:Sort substrings by length across all passages into $SC(a)$ in a descending
order.

6:$a_{\text{marked}}$: A dictionary that maps the the letters of the answer
$a$ to either "marked" or "unmarked".

7:$\text{IsMarked}(sc,a,a_{\text{marked}})$: Returns true if the string $sc$
has already been marked or marked partially inside answer $a$, and false
otherwise.

8:$\text{Mark}(sc,a,a_{\text{marked}})$: Marks the string in the answer as
"marked" under passage $sc_{p}$.

9:for $sc\in SC(a)$ do

10: if $!\text{IsMarked}(sc,a,a_{\text{marked}})$ then

11: $\text{Mark}(sc,a,a_{\text{marked}})$

12: end if

13:end for

14:Format marked spans stored in $a_{\text{marked}}$ with the answer $a$ to
form $a_{\text{span}}$.

15:return $a_{\text{span}}$

Citation Level | Cited Answer
---|---
Passage |  | " Johnny Panic and the Bible of Dreams " is a song by the British
---
band Tears for Fears. They also have international hit singles
such as: "Mothers Talk", "Shout", "Everybody Wants to Rule the World",
"Head over Heels", and "I Believe" [1][5].
Sentence |  | " Johnny Panic and the Bible of Dreams " is a song by the British
---
band Tears for Fears [1]. They also have international hit singles
such as: "Mothers Talk", "Shout", "Everybody Wants to Rule the World",
"Head over Heels", and "I Believe" [5].
Span |  | " [ 1 Johnny Panic and the Bible of Dreams ] " [ 1 is a song by the British
---
band Tears for Fears ] . They also have [ 5 international hit singles ]
such as: [ 5 "Mothers Talk", "Shout", "Everybody Wants to Rule the World",
"Head over Heels", and "I Believe". ]
Table 4: Citation Level examples, with each row highlighting the approach for citing the source passages. CoT Method | Answer Prefix
---|---
Passage |  | Lets analyze the input passages.
---
The only relevant passages to the question are passages 1, 5.
Thus, the final answer is:
Sentence |  | Lets analyze the relevant information from the input passages.
---
From passage [1], we know that: " Johnny Panic and the Bible of Dreams "
is a song by the British band Tears for Fears .
From passage [5], we know that: They also have international hit singles
such as: "Mothers Talk", "Shout", "Everybody Wants to Rule the World", "Head
over Heels", and "I Believe".
Thus, the final answer is:
Span |  | Lets analyze the relevant spans of information from the input passages.
---
From passage [1], the relevant spans are the following:
Johnny Panic and the Bible of Dreams
is a song by the British band Tears for Fears
From passage [5], the relevant spans are the following:
international hit singles
"Mothers Talk", "Shout", "Everybody Wants to Rule the World", "Head
over Heels", and "I Believe".
Thus, the final answer is:
Table 5: Examples of Chain-of-Thought prompt prefix for the various Levels.
The model is instructed/finetuned to produce the snippets per CoT method.
Answer the question by summarizing the given sources while explicitly copying
---
spans from the sources. When copying a span, use brackets and the respective
source number to indicate that this span was copied. Use explicit copying
as much as possible and for all factual statements, while preserving fluency.
Make sure to use all relevant sources and properly quote them. Here
are some examples:
Question: how much power does a wind turbine produce?
[1] Compact wind acceleration turbine: It is generally thought that since…
[2] Sustainable architecture: roof ledge. Small-scale rooftop wind turbines
have …
[3] Turby wind turbine: can because horizontal axis (HAWT) types cannot change
…
Quoted summary: One source states the [ 1 amount of power produced by a wind
turbine is proportional to the cube of the wind speed ] . Other sources state
[ 2 Turbines for residential scale use ] [ 2 produce electricity at a rate of
900 watts to
10,000 watts ] , and [ 3 is specified to generate power in winds of
between 4 m/s (9 mph, 7.8kts) and 14 m/s (31 mph, 27.2kts) ] .
Question: a component is what?
[1] Modular programming: in Dart, Go or Java) is sometimes used instead of …
[2] Physical body: the system at a point in time changes from identifying the
object to…
Quoted summary: A [ 1 component is a piece of a whole system ] . Also, [ 2 A
component is an object completely within the boundary of a containing object.
]
Question: what is the title of tears for fears song?
[1] Johnny Panic and the Bible of Dreams (song): …
[2] Mark Crew: Will Ruin Your Life. The album was produced at his studio. …
[3] Raoul and the Kings of Spain: Raoul and the Kings of Spain is the fifth
studio …
[4] Everybody Loves a Happy Ending: Everybody Loves a Happy Ending …
[5] Songs from the Big Chair: Songs from the Big Chair is the second studio
album …
Quoted summary: Song titles by [ 1 the British band Tears for Fears ] include
the
[ 5 1985 ] [ 5 international hit singles "Mothers Talk", "Shout", "Everybody
Wants to
Rule the World", "Head over Heels", and "I Believe" ] , [ 4 1989’s "The Seeds
of Love" ]
, the [ 1 1990 single "Advice for the Young at Heart" ] with [ 1 "Johnny Panic
and the Bible
of Dreams" ] on [ 1 the B-side ] , and " [ 2 Lemon To A Knife Fight ] ",
" [ 2 Cheetah Tongue ] ", and " [ 2 Turn ] " from [ 2 2018 ]
Table 6: Citation Level prompt example, containing 2 fewshot examples, with the relevant passages for each question provided. Level | CoT | Model | Answer Quality | Citation Quality
---|---|---|---|---
BERT | HEM | RL | AVG. | ALCE F1 | CSCA | DOC F1 | SEM-F1 | AVG.
Span |  | GPT-4 | 91.71 | 67.58 | 66.08 | 75.12 | 77.06 | 81.15 | 85.33 | 67.57 | 77.78
None | Flan | 92.20 | 65.41 | 72.12 | 76.58 | 85.92 | 94.26 | 86.71 | 71.94 | 84.71
| Mist. | 91.61 | 60.93 | 70.03 | 74.19 | 86.95 | 94.97 | 84.69 | 65.27 | 82.97
| GPT-4 | 92.06 | 66.95 | 68.55 | 75.85 | 78.83 | 89.88 | 86.54 | 72.23 | 81.87
Span | Flan | 92.06 | 67.97 | 71.99 | 77.34 | 81.47 | 95.36 | 86.51 | 72.86 | 84.05
| Mist. | 92.04 | 65.98 | 70.67 | 76.23 | 74.60 | 96.17 | 85.87 | 72.19 | 82.21
| GPT-4 | 91.87 | 66.88 | 68.76 | 75.84 | 84.76 | 69.69 | 86.38 | 70.52 | 77.84
Sent. | Flan | 92.15 | 71.46 | 71.64 | 78.42 | 86.71 | 90.02 | 87.79 | 74.41 | 84.73
| Mist. | 92.01 | 63.32 | 70.22 | 75.19 | 73.56 | 90.86 | 85.40 | 69.69 | 79.88
| GPT-4 | 91.91 | 66.89 | 67.67 | 75.49 | 79.21 | 86.94 | 86.49 | 69.87 | 80.63
Pass. | Flan | 92.17 | 70.84 | 72.47 | 78.49 | 83.36 | 94.14 | 88.01 | 77.07 | 85.64
| Mist. | 91.88 | 63.90 | 70.91 | 75.56 | 80.14 | 94.12 | 85.03 | 68.63 | 81.98
Sent. |  | GPT-4 | 91.47 | 62.48 | 62.18 | 72.04 | 82.79 | - | 85.32 | 59.20 | 75.77
None | Flan | 91.84 | 60.43 | 69.22 | 73.83 | 79.37 | - | 81.18 | 53.66 | 71.40
| Mist. | 91.56 | 58.51 | 67.63 | 72.57 | 82.15 | - | 82.10 | 54.00 | 72.75
| GPT-4 | 91.92 | 62.06 | 64.76 | 72.91 | 82.39 | - | 86.71 | 61.74 | 76.94
Span | Flan | 91.43 | 57.43 | 66.60 | 71.82 | 71.33 | - | 83.19 | 54.46 | 69.66
| Mist. | 91.93 | 65.04 | 69.69 | 75.55 | 73.14 | - | 84.48 | 59.29 | 72.30
| GPT-4 | 91.96 | 64.59 | 67.36 | 74.64 | 86.42 | - | 86.51 | 64.44 | 79.12
Sent. | Flan | 92.37 | 69.36 | 72.87 | 78.20 | 85.38 | - | 87.03 | 66.77 | 79.72
| Mist. | 91.94 | 64.45 | 70.40 | 75.60 | 75.48 | - | 84.81 | 59.97 | 73.42
| GPT-4 | 91.69 | 61.19 | 64.42 | 72.43 | 84.25 | - | 85.92 | 61.79 | 77.32
Pass. | Flan | 91.64 | 62.97 | 68.62 | 74.41 | 72.57 | - | 85.05 | 60.43 | 72.68
| Mist. | 91.99 | 62.90 | 69.37 | 74.76 | 75.39 | - | 85.18 | 60.07 | 73.55
Pass. |  | GPT-4 | 91.55 | 61.13 | 62.76 | 71.81 | 75.38 | - | 85.19 | 44.31 | 68.29
None | Flan | 91.69 | 62.10 | 68.26 | 74.02 | 69.33 | - | 83.31 | 40.41 | 64.35
| Mist. | 91.74 | 63.64 | 67.91 | 74.43 | 75.30 | - | 83.93 | 40.08 | 66.44
| GPT-4 | 91.90 | 61.66 | 64.70 | 72.75 | 71.60 | - | 86.17 | 43.41 | 67.06
Span | Flan | 91.24 | 54.87 | 65.56 | 70.56 | 65.51 | - | 81.98 | 38.76 | 62.09
| Mist. | 91.99 | 64.25 | 69.14 | 75.12 | 70.65 | - | 83.52 | 39.78 | 64.65
| GPT-4 | 91.84 | 62.43 | 65.32 | 73.20 | 70.63 | - | 86.41 | 42.59 | 66.54
Sent. | Flan | 92.43 | 65.17 | 72.61 | 76.73 | 73.00 | - | 84.47 | 42.05 | 66.51
| Mist. | 92.14 | 62.54 | 70.59 | 75.09 | 69.92 | - | 84.53 | 40.61 | 65.02
| GPT-4 | 91.69 | 60.40 | 63.88 | 71.99 | 73.28 | - | 85.64 | 44.12 | 67.68
Pass. | Flan | 91.32 | 61.82 | 66.14 | 73.09 | 58.18 | - | 86.16 | 39.20 | 61.18
| Mist. | 92.15 | 64.54 | 70.83 | 75.84 | 72.92 | - | 84.64 | 40.76 | 66.11
Table 7: Results for the GPT-4 with fewshot and the smaller fintuned models over the QuoteSUM dataset. Each row contains both the values of the metrics per citation level and CoT method, and the average score on per metric type on their right. Mist. denotes Misrtal-7B, and Flan denotes FlanT5-XXL. Level | CoT | Model | Answer Quality | Citation Quality
---|---|---|---|---
BERT | HEM | RL | AVG. | ALCE F1 | CSCA | DOC F1 | SEM-F1 | AVG.
Span | None | Flan | 91.91 | 62.17 | 57.08 | 70.39 | 96.97 | 91.88 | 64.70 | 43.80 | 74.34
Mist. | 91.24 | 64.16 | 53.28 | 69.56 | 80.61 | 94.51 | 59.86 | 38.47 | 68.36
Span | Flan | 91.80 | 62.83 | 56.50 | 70.38 | 95.69 | 93.58 | 64.90 | 44.72 | 74.72
Mist. | 91.88 | 63.77 | 56.63 | 70.76 | 93.41 | 96.47 | 65.30 | 43.84 | 74.75
Sent. | Flan | 92.25 | 62.72 | 58.80 | 71.26 | 95.67 | 86.94 | 65.50 | 44.75 | 73.22
Mist. | 91.83 | 64.31 | 56.85 | 70.99 | 93.38 | 91.62 | 63.40 | 41.18 | 72.40
Pass. | Flan | 91.90 | 56.24 | 58.77 | 68.97 | 95.93 | 90.09 | 64.90 | 44.55 | 73.87
Mist. | 91.90 | 63.69 | 56.73 | 70.77 | 93.73 | 95.36 | 65.50 | 43.05 | 74.41
Sent. | None | Flan | 92.83 | 62.77 | 58.78 | 71.46 | 96.53 | - | 62.00 | 44.49 | 67.68
Mist. | 92.06 | 60.26 | 54.90 | 69.08 | 91.67 | - | 63.62 | 42.99 | 66.09
Span | Flan | 92.37 | 62.20 | 56.42 | 70.33 | 95.92 | - | 65.30 | 44.74 | 68.65
Mist. | 92.54 | 61.66 | 57.24 | 70.48 | 94.75 | - | 64.50 | 44.09 | 67.78
Sent. | Flan | 92.81 | 63.31 | 58.82 | 71.65 | 96.55 | - | 65.15 | 45.29 | 69.00
Mist. | 92.59 | 64.11 | 58.04 | 71.58 | 96.95 | - | 65.45 | 45.01 | 69.14
Pass. | Flan | 92.92 | 65.69 | 59.45 | 72.69 | 96.73 | - | 64.70 | 45.59 | 69.01
Mist. | 92.67 | 63.88 | 58.20 | 71.58 | 96.05 | - | 66.40 | 45.14 | 69.19
Pass. | None | Flan | 92.96 | 64.38 | 59.93 | 72.42 | 95.20 | - | 64.30 | 45.94 | 68.48
Mist. | 92.29 | 62.06 | 56.32 | 70.22 | 92.54 | - | 62.80 | 42.75 | 66.03
Span | Flan | 92.17 | 55.69 | 57.29 | 68.38 | 96.10 | - | 65.00 | 44.59 | 68.56
Mist. | 92.67 | 63.72 | 58.04 | 71.48 | 95.40 | - | 66.40 | 44.77 | 68.86
Sent. | Flan | 92.68 | 64.91 | 58.02 | 71.87 | 96.70 | - | 65.80 | 45.44 | 69.31
Mist. | 92.49 | 64.12 | 57.16 | 71.26 | 96.30 | - | 64.80 | 44.20 | 68.43
Pass. | Flan | 92.97 | 65.20 | 59.49 | 72.55 | 96.60 | - | 67.40 | 46.10 | 70.03
Mist. | 92.75 | 64.03 | 58.14 | 71.64 | 96.80 | - | 65.60 | 44.78 | 69.06
Table 8: Results for the smaller finetuned models over the MSMARCO dataset.
Each row contains both the values of the metrics per citation level and CoT
method, and the average score on per metric type on their right. Mist. denotes
Misrtal-7B, and Flan denotes FlanT5-XXL.


# Reconnection-Driven Energy Cascade in Magnetohydrodynamic Turbulence

Chuanfei Dong1,2${}^{~{},{\dagger}}$, Liang Wang2,1, Yi-Min Huang2,1,
Luca Comisso3, Timothy A. Sandstrom4, Amitava Bhattacharjee1,2
1Princeton Plasma Physics Laboratory, Princeton University, Princeton, NJ
08540, USA
2Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
3Department of Astronomy and Columbia Astrophysics Laboratory,
Columbia University, New York, NY 10027, USA
4NASA Ames Research Center, Mountain View, CA 94043, USA
Corresponding author. Email<EMAIL_ADDRESS>authors contributed
equally to this work.

> Magnetohydrodynamic turbulence regulates the transfer of energy from large
> to small scales in many astrophysical systems, including the solar
> atmosphere. We perform three-dimensional magnetohydrodynamic simulations
> with unprecedentedly large magnetic Reynolds number to reveal how rapid
> reconnection of magnetic field lines changes the classical paradigm of the
> turbulent energy cascade. By breaking elongated current sheets into chains
> of small magnetic flux ropes (or plasmoids), magnetic reconnection leads to
> a new range of turbulent energy cascade, where the rate of energy transfer
> is controlled by the growth rate of the plasmoids. As a consequence, the
> turbulent energy spectra steepen and attain a spectral index of -2.2 that is
> accompanied by changes in the anisotropy of turbulence eddies. The
> omnipresence of plasmoids and their consequences on, e.g., solar coronal
> heating, can be further explored with current and future spacecraft and
> telescopes.

## INTRODUCTION

Understanding the transfer of energy in magnetohydrodynamic (MHD) turbulence
is crucial for tackling many outstanding astrophysical problems, such as solar
and stellar coronal heating, star formation, cosmic-ray transport, and the
interstellar medium evolution. For more than half a century, it has been
widely accepted that the energy cascade in turbulent plasmas such as the Sun’s
atmosphere is controlled by magnetohydrodynamic wave interactions [1, 2, 3].
However, one essential feature of magnetohydrodynamic turbulence is the
ubiquitous presence of sheets of intense electric current (known as current
sheets), which are preferential locations for rapid breaking and reconnection
of magnetic field lines [4, 5], a fundamental physical process in magnetized
plasmas whereby stored magnetic energy is converted into heat and kinetic
energy of charged particles. It is yet an open question whether magnetic
reconnection can significantly change the transfer of energy from large to
small scales in a wide variety of astrophysical systems.

One of the most common approaches to investigating the effect of a physical
process on the turbulent energy cascade is to study the associated energy
spectra. It has been widely observed by different spacecraft and telescopes
that turbulent energy spectra can break at either visco-resistive or kinetic
scales [6, 7, 8]. Recent analytic studies suggest that magnetic reconnection
may also break the turbulent energy spectra and thus create a new range of
energy transfer when the growth time scale of the magnetic flux ropes (or
plasmoids), $1/\gamma_{p}$, becomes much shorter than the nonlinear eddy
turnover time, $\tau_{nl}$ [9, 10, 11, 12, 13, 14]. Plasmoids develop in
intense and elongated current sheets undergoing reconnection, and in MHD
turbulence, the spatial scales and aspect ratios of those current sheets are
controlled by the magnetic Reynolds number $R_{m}$, which quantifies the
relative magnitudes of plasma convection and resistive diffusion. At
sufficiently large $R_{m}$ ($>10^{5}$), as in the solar atmosphere, magnetic
reconnection is expected to be a ubiquitous phenomenon, as the thinning of
current sheets can lead to a copious formation of plasmoids via the tearing
instability [15, 16, 17, 18, 19], potentially changing the energy transfer
across scales.

Due to the complex, nonlinear nature of the turbulent energy cascade, direct
numerical simulations (DNS) are likely the best means to investigate the role
of magnetic reconnection in the energy transfer across scales and how it
changes the turbulent energy spectra. To date, no evidence for a new range of
the turbulent energy cascade due to magnetic reconnection has been provided by
DNS in realistic three dimensions (3D) [20, 21]. Such DNS are extremely
challenging, mainly due to the high grid resolution required to capture the
fine structure of the omnipresent current sheets in a turbulent plasma at
large $R_{m}$. In addition, MHD turbulence and magnetic reconnection are known
to behave differently in 2D and 3D [22, 20], therefore 3D DNS with large
$R_{m}$ are essential to fundamentally address this question.

## RESULTS

Here, we present the world’s largest 3D MHD turbulence simulation (at a cost
of $\sim 200$ million CPU hours) that self-consistently produces myriad fine
current sheets. An elongated ($1\times 1\times 2$) periodic box with $\sim
10000\times 10000\times 5000$ grid cells was adopted to resolve the thin
current sheets that develop at large magnetic Reynolds number, $R_{m}=10^{6}$
(and the effective magnetic Reynolds number $R_{m,{\rm eff}}=2\times 10^{5}$
based on the energy injection scales). We initialized the simulations with
uncorrelated, equipartitioned velocity and magnetic field fluctuations
superimposed by a strong mean magnetic field in the elongated direction (see
Materials and Methods for the detailed model setup). Such field configurations
are common in a variety of astrophysical systems such as solar/stellar coronae
and the interstellar medium. Compared with earlier 3D MHD turbulence
simulations, the most prominent development here is that the current study
reaches an unprecedented high-$R_{m}$ regime, such that the ubiquitous
reconnecting current sheets in a turbulent bath become unstable to tearing
instability.

The turbulent structures from our large-scale MHD simulation are visualized in
Fig. 1. In Fig. 1a, volume rendering of the current density $|\mathbf{J}|$
depicts MHD turbulence in the entire simulation domain at the fully developed
stage. Specifically, strong current sheets are ribbon-like (with large aspect
ratios) and are aligned with the mean magnetic field $B_{z0}$ due to the
parallel coherence of perturbations. Meanwhile, they are also current sheets
in the perpendicular plane aligned with in-plane perturbed magnetic fields.
These ribbon-shaped current sheets are subject to the tearing instability
during their dynamical evolution [15, 16, 17, 18, 19] (see Fig. 1d). The
current density $|\mathbf{J}|$ in different x-y slices can be seen in movie
S1. One of these current sheets undergoing magnetic reconnection is
highlighted in a zoomed-in subdomain (Figs. 1b and 1c), within which “ripples”
induced by the formation of plasmoids/magnetic flux ropes are present on the
current density isosurfaces. We cut a 2D slice across the embedded 3D magnetic
flux ropes for $J_{z}$ (Fig. 1c; also see Fig. 5 in Materials and Methods for
a selected flux-rope bundle at different viewing angles). The current sheet in
the subdomain exhibits features similar to those observed in previous
simulations of a single 3D reconnecting current layer [23, 24, 20]. This
subdomain ($0.025\%$ of the complete volume) is representative of the types of
coherent structures that develop throughout the entire domain. Reconnection-
produced magnetic flux ropes in 3D exhibit a more complex morphology than
their 2D counterpart. This is illustrated in Fig. 1d, in which the out-of-
plane current density, $J_{z}$, in one arbitrary $x$-$y$ slice, points out
that the reconnecting current sheet structures in 3D are fractured (or
disrupted) and differ significantly from the island-like morphology in 2D [23,
24, 25, 26, 27], despite the similar size of the flux ropes.

Our simulation identifies a novel mechanism for energy transfer through the
breakup of reconnecting current sheets into smaller fragments, which occurs
generically in the large-$R_{m}$ regime studied here (Fig. 1d). The role of
reconnection in controlling the energy transfer is supported by its imprint on
the energy spectrum of the turbulent cascade in Fig. 2. At large scales,
$k_{\bot}\lesssim k_{*}$, both the magnetic energy spectrum $E_{B}(k_{\bot})$
and the kinetic energy spectrum $E_{U}(k_{\bot})$ follow a power law with a
slope of $-3/2$ (see dashed fitted lines), in agreement with expectations for
the inertial range of a strong turbulent cascade mediated by Alfvén waves when
accounting for a reduction of nonlinearity due to dynamic alignment [28]. At
scales $k_{*}\lesssim k_{\bot}\lesssim k_{\eta}$, current sheets in this
tearing-mediated regime (i.e., reconnection-driven regime, as highlighted by
the shaded region in Fig. 2) exhibit numerous reconnecting structures, or
magnetic flux ropes, as shown in the inset. As a result, flux rope formation
through reconnection dominates over MHD wave interactions in controlling
turbulent energy transfer. Consequently, the energy spectrum in this range
(termed as sub-inertial range hereafter) becomes steepened and is
characterized by a spectral index of $-11/5$ when $R_{m}$ is sufficiently
large, which is broadly consistent with the theoretical predictions [9, 10,
11, 12]. It is noteworthy that the turbulent energy spectrum breaks at
$k_{*}\approx 1000$ in Fig. 2, which also agrees with the theoretical
prediction of tearing-disruption scale $k_{*}=R_{m,{\rm eff}}^{4/7}\approx
1000$ [9, 10, 11, 12, 29], where $R_{m,{\rm eff}}=2\times 10^{5}$ is the
effective magnetic Reynolds number due to the energy injection scales as the
energy spectrum peaks at $k_{\bot}\approx 30$. For comparison, we ran another
3D simulation with a relatively low magnetic Reynolds number, $R_{m}=8\times
10^{4}$, while keeping the initial condition identical. In this lower $R_{m}$
simulation, we observed that flux ropes are essentially absent (see Fig. S1 in
Supplementary Materials), which corroborates the lack of any imprint from
reconnection-driven cascades in the corresponding energy spectrum (the blue
curves in Fig. 2). We, therefore, conclude that a sufficiently large $R_{m}$
($>10^{5}$) is one of the prerequisites for evolving current sheets to become
tearing-unstable in a turbulent bath, consequently leading to the copious
formation of flux ropes before reaching the scale given by $k_{\eta}$.

We further investigate changes in the shape (i.e., anisotropy) of statistical
eddies in the sub-inertial range, one key element in making reconnection-
driven energy cascade dominate over energy cascade through MHD wave
interactions. In the top panels of Fig. 3, we report the dependence of the
half width $\xi$ (panel a) and aspect ratio $\xi/\lambda$ (panel b) of the
statistical eddies on the half thickness $\lambda$ based on second-order
structure functions (see Fig. 6 in Materials and Methods). While the scaling
$\xi/\lambda\propto\lambda^{-0.2}$ is observed in the traditional inertial
range, in the sub-inertial range the scaling changes to approximately
$\xi/\lambda\propto\lambda^{-0.1}$. The scaling is in contrast to the previous
theoretical prediction $\xi/\lambda\propto\lambda^{4/5}$ [11]. Similarly, the
bottom panels of Fig. 3 depict the relation between the half length
$\ell_{\parallel}$ and the half thickness $\lambda$. A scaling relation
$\ell_{\parallel}\propto\lambda^{1/2}$ is observed in the inertial range as
expected from the theory. However, the scaling in the sub-inertial range,
$\ell_{\parallel}\propto\lambda^{2/3}$ (and
$\ell_{\parallel}/\lambda\propto\lambda^{-1/3}$), again deviates from the
theoretical estimates [11]. It is noteworthy that the theoretical predictions
$\xi/\lambda\propto\lambda^{4/5}$ and
$\ell_{\parallel}/\lambda\propto\lambda^{1/5}$ in the sub-inertial range
require the assumption $\theta_{t}\sim\delta/\zeta$, where $\theta_{t}$ is the
alignment angle, and $\delta$ and $\zeta$ denote the inner layer width and the
wavelength of the fastest growing tearing mode, respectively [11]. The
difference between the theoretical predictions and the numerical experiments
indicates that the assumption $\theta_{t}\sim\delta/\zeta$ needs further
examination.

In order to address the underlying cause leading to this new range of the
turbulent energy cascade, we also investigate the difference in the energy
transfer between the inertial and sub-inertial ranges. For this purpose, we
calculate the cylindrical shell-to-shell magnetic energy transfer function,
$T_{bb}\left(Q,\
K\right)=-\int\mathbf{B}_{K}\cdot\left(\mathbf{u}\cdot\nabla\right)\mathbf{B}_{Q}{\rm
d}^{3}x$, where $\mathbf{B}_{K}$ contains all Fourier modes
$\tilde{\mathbf{B}}\left(\mathbf{k}\right)$ in the $K$-th perpendicular
wavenumber shell, $k_{K}<k_{\perp}<k_{K+1}$, and $K$ is an integral shell
number; $\mathbf{B}_{Q}\left(\mathbf{x}\right)$ is defined in the same way for
the $Q$-th shell (see Materials and Methods for a detailed description). The
red pixels in Fig. 4a represent positive magnetic energy transfer from the
$Q$-th shell to the $K$-th shell, and the blue pixels the opposite. The
diagonal pixels are white (i.e., empty) since there is no self-magnetic-
energy-transfer by construction. The $T_{bb}$ distribution on the
$k_{Q}$-$k_{K}$ phase plane in Fig. 4a is characterized by main features such
as local, forward transfer. Figs. 4b and 4c represent vertical cuts of
$T_{bb}$ normalized to

$N_{{\rm
inertial}}\sim\frac{E_{Q}}{\tau_{A}}(k_{Q}L)^{-1/4}\sim\frac{E_{Q}^{3/2}k_{Q}}{\left(\rho
V\right)^{1/2}}(k_{Q}L)^{-1/4}$ (1)

and

$N_{{\rm tearing}}\sim\gamma_{p}E_{Q}\sim\frac{\eta^{1/2}}{\left(\rho
V\right)^{1/4}}E_{Q}^{5/4}k_{Q}^{3/2},$ (2)

where $\tau_{A}$ is the Alfvénic time scale and $\gamma_{p}$ the linear growth
rate of the tearing instability (see Materials and Methods for details). The
factor $(k_{Q}L)^{-1/4}$ in Eq. (1) takes into account the reduction of
nonlinearity due to dynamic alignment [28]. Within the inertial range, energy
transfer is expected to be self-similar, i.e., the energy transfer from the
$Q$-th shell to the $K$-th shell operates in the same way as the transfer from
the ($Q+1$)-th shell to the ($K+1$)-th shell. This is consistent with the
observation in Fig. 4b that, when normalized to $N_{{\rm inertial}}$, the
three cuts in the inertial range (solid lines in Fig. 4a) strongly overlap. In
contrast, the three cuts in the tearing-mediated range (dotted curves in Fig.
4a) do not exhibit self-similarity in Fig. 4b, due to the underlying mechanism
being different than the classical, inertial-range energy cascade. However,
when $T_{bb}$ is normalized by $N_{{\rm tearing}}$, i.e., $T_{bb}/N_{{\rm
tearing}}$, the cuts in the tearing-mediated range (dashed curve) become self-
similar (Fig. 4c), thus confirming the existence of a sub-inertial range
within which the energy transfer is controlled by the tearing instability in
reconnecting current sheets.

## DISCUSSION

The present study suggests that the energy transfer in, e.g., the solar
atmosphere at small scales can be fundamentally different from the classic
paradigm of the turbulent energy cascade controlled by MHD wave interactions.
Our calculated magnetic energy spectrum captures a new range of reconnection-
driven energy cascade with a spectral index of $-2.2$. It is interesting to
point out that the line-of-sight magnetic field observations of a coronal hole
using magnetograms acquired by the Near InfraRed Imaging Spectrapolarimeter
(NIRIS) operating at the Goode Solar Telescope of the Big Bear Solar
Observatory (BBSO) showed that the turbulent magnetic energy spectra also
exhibit a new range with a spectral index of -2.2 [7] (also see Fig. S2 in
Supplementary Materials), but it requires further observations to rule out the
coincidence and then to identify the underlying mechanism that leads to the
spectral break. On the other hand, similar observations acquired by the
Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamic Observatory
(SDO) with a relatively low spatial resolution did not present a steepening in
the corresponding energy spectra, suggesting that a high spatial resolution
might be required to reveal this new range of energy transfer in the turbulent
solar atmosphere. To address the discrepancy between the observations of NIRIS
and HMI, we illustrated, with the same dataset in Fig. S2, that the spatial
resolution of observation can significantly influence the scientific findings.

In addition, spacecraft such as Solar Orbiter may also be able to reveal this
newly identified sub-inertial range in magnetic energy spectra through remote
sensing observations. Recently, Extreme Ultraviolet Imager (EUI) onboard Solar
Orbiter observed transient small-scale brightenings prevalent in the corona of
the quiet Sun termed ‘campfires’ [30]. It has been proposed that the majority
of campfire events observed by EUI are driven by magnetic reconnection, which
may play an important role in the coronal heating of the quiet Sun [31, 32].
Our 3D simulation results in the large-$R_{m}$ regime suggest that magnetic
reconnection is a ubiquitous process in the turbulent solar corona where the
$R_{m}$ is even larger and thus the current sheets can thin down to much
smaller scales and form the fractal structures within which copious formation
of plasmoids occurs. Hence, there should be many more reconnection sites than
observed, which can be revealed by (future) high-resolution extreme
ultraviolet images.

In addition to the solar coronal heating [31, 30], the omnipresence of
plasmoids associated with magnetic reconnection, identified in this study,
also has broad implications and consequences on particle acceleration in the
solar/stellar coronae, accretion disks, and jets from compact objects [33,
34], and on filament formation in the _Herschel_ maps of the Orion A giant
molecular cloud [35]. Moreover, the ubiquitous presence of dense plasmoids in
the interstellar medium can also play an important role in pulsar
scintillation [36]. Indeed, all the aforementioned astrophysical systems are
associated with large magnetic Reynolds numbers ($R_{m}>10^{5}$), and thus
could be characterized by the tearing-mediated turbulence identified in the
present work.

Figure 1: Reconnecting current sheets and magnetic flux ropes in MHD
turbulence. (a) Volume rendering of the current density $|\mathbf{J}|$ in the
entire domain at a stage when turbulence is fully developed. Myriad of current
sheets is evident in the plane perpendicular to the mean magnetic field
$B_{z0}$. Panels (b) and (c) depict one reconnecting current sheet and the
embedded flux-ropes in a small subdomain (within the boundaries
$[-0.45,-0.4]\times[0.45,0.5]\times[0.95,1.05]$). Panel (b) shows the volume
rendering of $|\mathbf{J}|$, while panel (c) displays magnetic field lines
(colored by $|\mathbf{B}|$) associated with the featured current sheet
(including magnetic flux ropes) and a $x$-$y$ slice view of the current
density component, $J_{z}$, along the mean magnetic field. (d) Out-of-plane
current density $J_{z}$ in a $x$-$y$ slice (at $z=-1$) of the 3D turbulence
simulation (left) compared with the corresponding result from a 2D simulation
(right). Copious formation of magnetic flux ropes/plasmoids occurs in both 3D
and 2D simulations despite the different morphology. Zoomed-in subdomains are
used to illustrate the increased morphological complexity that characterizes
the 3D simulation.

Figure 2: Steepening of energy spectra in reconnection-driven energy cascade.
Field-perpendicular magnetic and kinetic energy spectra, $E_{B}(k_{\bot})$ and
$E_{U}(k_{\bot})$, showing a standard inertial range with a slope
$E_{B,U}(k_{\bot})\propto k_{\bot}^{-3/2}$ and a reconnection-driven (or
tearing-mediated) sub-inertial range with a slope $E_{B,U}(k_{\bot})\propto
k_{\bot}^{-11/5}$ for large magnetic Reynolds number $R_{m}$ (orange curves).
The sub-inertial range is absent in an equivalent 3D simulation with lower
$R_{m}$ (blue curves). The shaded area emphasizes the reconnection-driven sub-
inertial range, with wavenumbers corresponding to the typical transverse
scales of the flux ropes. The left edge of the shaded range is identified at
$k_{*}\approx 10^{3}\approx R_{m,{\rm eff}}^{7/4}$, where $R_{m,{\rm
eff}}\approx 2\times 10^{5}$ is computed with the injection scale at $k\approx
30$. The right edge, marked by $k_{\eta}$, is the dissipation scale defined
such that $\eta\int_{0}^{k_{\eta}}E_{B}(k_{\bot})k_{\bot}^{2}{\rm d}k_{\bot}$
accounts for approximately half of the resistive dissipation power
$\eta\langle J^{2}\rangle$. The small volume rendering inset illustrates a
typical reconnecting current layer, that is ubiquitous in this range, with the
transverse size of an embedded flux rope annotated in red.

Figure 3: Scale-dependent anisotropy of statistical eddies. (a) The
perpendicular half width $\xi$ of the statistical eddies as a function of the
perpendicular half thickness $\lambda$. This function is obtained by equating
the second-order $\mathbf{B}$-trace structure function measured along $\xi=0$
and $\lambda=0$ axes in Fig. 6b. (b) Aspect ratio $\xi/\lambda$ of the
statistical eddies as a function of $\lambda$. (c) and (d) are the scaling
relation between the parallel half length $\ell_{\parallel}$ and $\lambda$,
obtained in a similar way but in the local coordinate plane $\xi=0$. The
fitted dashed and dot-dashed lines show two scaling laws in the traditional
inertial range and the tearing-mediated sub-inertial range of the turbulent
energy cascade.

Figure 4: Self-similarity of normalized cylindrical shell-to-shell energy
transfers. (a) Magnetic-to-magnetic energy transfer functions
$T_{bb}\left(Q,K\right)=-\int\mathbf{B}_{K}\cdot\left(\mathbf{u}\cdot\nabla\right)\mathbf{B}_{Q}{\rm
d}^{3}{x}$ in the $k_{Q}$-$k_{K}$ plane. (b)-(c) Normalized energy transfer
function values along different vertical cuts in specific giver $k_{Q}$
shells. The cut locations are marked in (a), with solid-line cuts chosen in
the inertial range and dotted-line cuts in the reconnection-driven (or
tearing-mediated) sub-inertial range. The normalization used in (b) and (c),
$N_{{\rm inertial}}=E_{Q}^{3/2}k_{Q}(k_{Q}L)^{-1/4}/(\rho V)^{1/2}$ and
$N_{{\rm tearing}}=\eta^{1/2}E_{Q}^{5/4}k_{Q}^{3/2}/(\rho V)^{1/4}$ (see
Materials and Methods), are appropriate for the inertial range and the
tearing-mediated sub-inertial range, respectively. The horizontal coordinates
in (b) and (c) are the shell number differences (integers), $K-Q$.

## MATERIALS AND METHODS

### Model description

Following procedures previously described in [25], the governing equations of
our numerical model are the dimensionless visco-resistive MHD equations:

$\partial_{t}\rho+\nabla\cdot(\rho{\bf{u}})=0\,,$ (3)
$\partial_{t}(\rho{\bf{u}})+\nabla\cdot(\rho{\bf{uu}})=-\nabla\left(p+B^{2}/2\right)+\nabla\cdot({\bf{BB}})+\nu{\nabla^{2}}(\rho{\bf{u}})\,,$
(4)
$\partial_{t}p+\nabla\cdot(p{\bf{u}})=(\gamma-1)\left({-p\nabla\cdot{\bf{u}}+\eta{\bf{J}}^{2}}\right)\,,$
(5)
$\partial_{t}{\bf{B}}=\nabla\times({\bf{u}}\times{\bf{B}}-\eta{\bf{J}})\,,$
(6)

where $\rho$, ${\bf{u}}$ and $p$ are the mass density, velocity, and pressure
of the plasma, respectively; ${\bf{B}}$ is the magnetic field and
${\bf{J}}=\nabla\times{\bf{B}}$ denotes the electric current density. The
kinematic viscosity and the magnetic diffusivity are denoted as $\nu$ and
$\eta$, respectively, while $\gamma$ is the adiabatic index.

### Model setup

We solve Eqs. (3)-(6) using the BATS-R-US MHD code [37] by adopting a fifth-
order scheme [38] in a domain $\\{(x,y,z):-L_{0}/2\leq x,y\leq
L_{0}/2,-L_{0}\leq z\leq L_{0}\\}$, where $L_{0}$ is set to unity. Periodic
boundary conditions are employed in all three directions. Lengths are
normalized to the box size $L_{0}$, velocities to the characteristic Alfvén
speed $V_{A}$, and time to $L_{0}/V_{A}$. We initialize the simulations by
placing uncorrelated, equipartitioned velocity and magnetic field fluctuations
in Fourier harmonics as follows,

$\displaystyle\psi$ $\displaystyle=$
$\displaystyle\sum\limits_{l,m,n}(a_{mn}/2\pi)\sin(2\pi mx/L_{x}+2\pi
ny/L_{y}+2\pi lz/L_{z}+\phi_{lmn})$ (7) $\displaystyle B_{x}$ $\displaystyle=$
$\displaystyle\frac{\partial\psi}{\partial
y}=\sum\limits_{l,m,n}a_{mn}(n/L_{y})\cos(2\pi mx/L_{x}+2\pi ny/L_{y}+2\pi
lz/L_{z}+\phi_{lmn})$ (8) $\displaystyle B_{y}$ $\displaystyle=$
$\displaystyle-\frac{\partial\psi}{\partial
x}=-\sum\limits_{l,m,n}a_{mn}(m/L_{x})\cos(2\pi mx/L_{x}+2\pi ny/L_{y}+2\pi
lz/L_{z}+\phi_{lmn})$ (9) $\displaystyle B_{z}$ $\displaystyle=$
$\displaystyle B_{z0}$ (10)

where we set

$a_{mn}=\frac{\sqrt{2}B_{0}}{N^{1/2}(m^{2}/L_{x}^{2}+n^{2}/L_{y}^{2})^{1/2}}$
(11)

and $L_{x}=L_{y}=L_{0}$, and $L_{z}=2L_{0}$. $\phi_{lmn}$ denotes random
phases for each mode. The summation is over the range $0\leq m\leq m_{max}$,
$-n_{max}\leq n\leq n_{max}$. However, for $m=0$, we limit to the range to
$0<n<n_{max}$. The range of $l$ is always $-l_{max}\leq l\leq l_{max}$. The
total number of modes are

$N=(2l_{max}+1)[(2m_{max}+1)(2n_{max}+1)-1]/2.$ (12)

We only sum over half of the Fourier space because the $(m,n)$ mode and the
$(-m,-n)$ mode are not independent. Here, we choose $m_{max}=n_{max}=10$ and
$l_{max}=5$. The above setup leads to

$\langle B_{\perp}^{2}\rangle=B_{0}^{2},$ (13)

where $\langle...\rangle$ represents the spatial average.

For the velocity $\mathbf{u}$, we use the similar expressions:

$\displaystyle u_{x}$ $\displaystyle=$
$\displaystyle\sum\limits_{l,m,n}b_{mn}(n/L_{y})\cos(2\pi mx/L_{x}+2\pi
ny/L_{y}+2\pi lz/L_{z}+\phi^{\prime}_{lmn})$ (14) $\displaystyle u_{y}$
$\displaystyle=$ $\displaystyle-\sum\limits_{l,m,n}b_{mn}(m/L_{x})\cos(2\pi
mx/L_{x}+2\pi ny/L_{y}+2\pi lz/L_{z}+\phi^{\prime}_{lmn})$ (15) $\displaystyle
u_{z}$ $\displaystyle=$ $\displaystyle 0$ (16)

where we set

$b_{mn}=\frac{\sqrt{2}u_{0}}{N^{1/2}(m^{2}/L_{x}^{2}+n^{2}/L_{y}^{2})^{1/2}}.$
(17)

Similarly, we have

$\langle u_{\perp}^{2}\rangle=u_{0}^{2}.$ (18)

The plasma density and pressure are initially set to constant values
$\rho=1.0$ and $p=1.6$, respectively. The constants $B_{0}$ and $u_{0}$
determine the strength of the initial velocity and magnetic field
fluctuations. In this work, we set $B_{0}=u_{0}=1$, which gives the initial
turbulent energy $E=\frac{1}{2}\langle
u_{\perp}^{2}+B_{\perp}^{2}\rangle=\frac{1}{2}(u_{0}^{2}+B_{0}^{2})=1$. The
mean magnetic field $B_{z0}=4$ leads to a plasma beta (ratio of the plasma
pressure to the magnetic pressure) $\beta\approx 0.2$. In the main 3D
simulation presented here, the magnetic Reynolds number is
$R_{m}=u_{0}L_{0}/\eta=10^{6}$, where $\eta=10^{-6}$ denotes the magnetic
diffusivity. The effective magnetic Reynolds number is $R_{m,{\rm
eff}}=2\times 10^{5}$ due to the energy injection scales as the energy
spectrum peaks at $k_{\bot}\approx 30$ (see Fig. 2). We set the viscosity
$\nu=10^{-6}$, such that the magnetic Prandtl number $P_{m}=\nu/\eta=1$. This
large $R_{m}$ value allows the copious formation of flux ropes within the
current sheets, as verified by visual inspection of the simulation data in a
tiny subdomain (see example in Fig. 5). For the 2D simulations, we adopt
similar initial setups but ignore the third dimension in $z$. For both 3D and
2D studies, we also ran one additional simulation with $R_{m}=8\times 10^{4}$
and a period box of $\sim 2000\times 2000\times 1000$ and $\sim 2000\times
2000$ grid cells, respectively (see Fig. S1 in Supplementary Materials).
Compared with the main 3D simulation (at a cost of about 200 million CPU hours
with approximately 0.5 trillion grid cells), the computational cost of the
rest simulations is negligible.

Figure 5: 3D magnetic flux ropes with elongated current sheets. Flux ropes
passing through a plasmoid located near the center of the current sheet
displayed in the intermediate 2D slice. The chosen views are along the $+x$
(panel a) and $-x$ (panel b) directions, respectively. Magnetic field lines
are colored accordingly to the magnetic field magnitude $|\mathbf{B}|$, while
the different $x$-$y$ slices are colored by the value of the electric current
density component, $J_{z}$, in the direction of the mean magnetic field.

### Second-order structure function

The anisotropy of statistical eddies shown in Fig. 3 is calculated based on a
second-order $\mathbf{B}$-trace structure function. To calculate this
structure function, we randomly sample a large number of data-point pairs. For
each pair of points, 1 and 2, the displacement
$\delta\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$ is projected onto a local,
scale-dependent coordinate system $\left(\ell_{\parallel},\ \xi,\
\lambda\right)$: $\hat{\ell}_{\parallel}$ is along the local mean field
$\mathbf{B}_{{\rm local}}=\left(\mathbf{B}_{1}+\mathbf{B}_{2}\right)/2$,
$\hat{\xi}$ is along the local perpendicular fluctuation field
$\delta\mathbf{B}_{\perp,{\rm local}}=\mathbf{B}_{{\rm
local}}\times\left[\left(\mathbf{B}_{1}-\mathbf{B}_{2}\right)\times\mathbf{B}_{{\rm
local}}\right]/|\mathbf{B}_{{\rm local}}|^{2}$, and $\hat{\lambda}$ completes
the right-handed coordinate system. We then accumulate the contributions from
all sampling pairs to compute the $\mathbf{B}$-trace structure function
$S_{2}=\left\langle\left|\mathbf{B}_{1}-\mathbf{B}_{2}\right|^{2}\right\rangle$,
an ensemble-averaged function of $\left(\ell_{\parallel},\ \xi,\
\lambda\right)$. As an example, here we show the cross-section of $S_{2}$ in
the $\ell_{\parallel}=0$ plane. The consequent perpendicular structure
function $S_{2\perp}\left(\xi,\ \lambda\right)$ is illustrated in Fig. 6a,
which exhibits clear anisotropy between the perpendicular half width, $\xi$,
of the statistical eddies and the perpendicular half thickness, $\lambda$. It
is immediately observed that Fig. 6a also depicts the perpendicular cross-
section of statistical eddies at different scales.

To further quantify the anisotropy, we investigate the aspect ratio of the
perpendicular eddies, $\xi/\lambda$. To this end, we measure $S_{2\perp}$
along the two perpendicular axes, $S_{2\perp}\left(\xi;\ \lambda=0\right)$ and
$S_{2\perp}\left(\lambda;\ \xi=0\right)$, as shown in Fig. 6b, and then find a
mapping between $\xi$ and $\lambda$ by equating $S_{2\perp}$ along the two
curves. This gives the result illustrated in Fig. 3.

Following a similar procedure, we also obtain the relation between the
parallel half length, $\ell_{\parallel}$, and the perpendicular half
thickness, $\lambda$, in Fig. 3 based on the structure function $S_{2}$ in the
$\xi=0$ plane.

Figure 6: Anisotropy of the statistical eddies perpendicular to the local mean
magnetic field. (a) Cross sections of the statistical eddies at different
values from the second-order structure function
$S_{2\perp}=\left<\left|\mathbf{B}_{1}-\mathbf{B}_{2}\right|^{2}\right>$. (b)
$S_{2\perp}$ as a function of the perpendicular half width $\xi$ or the
perpendicular half thickness $\lambda$ of the statistical eddies.

### Shell-to-shell energy transfers in the inertial and tearing-mediated sub-
inertial ranges

To further investigate the role of tearing instability in the turbulent energy
cascade, we calculated the cylindrical shell-to-shell magnetic energy transfer
function [39]

$T_{bb}\left(Q,\
K\right)=-\int\mathbf{B}_{K}\cdot\left(\mathbf{u}\cdot\nabla\right)\mathbf{B}_{Q}{\rm
d}^{3}x.$ (19)

Here,

$\mathbf{B}_{K}\left(\mathbf{x}\right)=\int_{k_{K}<k_{\perp}<k_{K+1}}\tilde{\mathbf{B}}\left(\mathbf{k}\right)e^{i\mathbf{k}\cdot\mathbf{x}}{\rm
d}^{3}k$ (20)

contains all Fourier modes $\tilde{\mathbf{B}}\left(\mathbf{k}\right)$ in the
$K$-th perpendicular wavenumber shell, $k_{K}<k_{\perp}<k_{K+1}$, and $K$ is
an integral shell number. $\mathbf{B}_{Q}\left(\mathbf{x}\right)$ is defined
in the same way for the $Q$-th shell. $T_{bb}\left(Q,K\right)$ then gives the
transfer rate of magnetic energy from the $Q$-th shell to the $K$-th shell
[40, 41], or, correspondingly, from the spatial scale $\sim 1/k_{Q}$ to the
scale $\sim 1/k_{K}$. We adopted 100 logarithmic bins along the $K$ and $Q$
directions between $k_{\perp}=2\pi$ and $k_{\perp}=10^{4}\pi$. The resulting
$T_{bb}$ values from our simulation on the $k_{Q}$-$k_{K}$ plane is shown in
Fig. 4a. The dominance of red pixels above the diagonal where $k_{K}>k_{Q}$
confirms the forward transfer from larger to smaller scales. Fig. 4a also
shows that transfers occur primarily close to the diagonal, indicating that
the energy transfer is mostly local.

### Normalization of the cylindrical shell-to-shell energy transfer rate

In the inertial range, we assume that Boldyrev’s turbulence theory with
dynamic alignment holds [28]. In this picture, turbulence eddies are
anisotropic in all three directions, with dimensions $\xi$ and $\lambda$
perpendicular to the local field, and $\ell_{\parallel}$ along the local
field. These scales are related as $\xi\sim L(\lambda/L)^{3/4}$ and
$\ell_{\parallel}\sim L(\lambda/L)^{1/2}$, where $L$ is the large-scale
length. The magnetic fluctuation $B_{\lambda}\sim B_{0}(\lambda/L)^{1/4}$ at
scale $\lambda$, where $B_{0}$ is the large-scale magnetic field. The energy
cascade of eddies at scale $\lambda$ occurs on the time scale
$\tau_{\lambda}\sim\ell_{\parallel}/V_{A0}\sim\xi/V_{A\lambda}$. Let $E_{Q}$
be the magnetic energy in the $Q$-th shell (i.e., the perpendicular wavenumber
$k_{\bot}$ satisfies $k_{Q}<k_{\perp}<k_{Q+1}$) and $V$ be the volume of the
domain. The magnetic energy density $E_{Q}/V$ on the $\lambda_{Q}\sim 1/k_{Q}$
scale is proportional to $B_{Q}^{2}/2$; therefore, the Alfvén speed
$V_{AQ}=B_{Q}/\sqrt{\rho}\sim E_{Q}^{1/2}/(V\rho)^{1/2}$ and the cascade time
$\tau_{Q}\sim\xi_{Q}/V_{AQ}\sim(k_{Q}L)^{1/4}(V\rho)^{1/2}/E_{Q}^{1/2}k_{Q}$.
This leads to the normalization $N_{{\rm
inertial}}\sim{E_{Q}}/{\tau_{Q}}\sim(k_{Q}L)^{-1/4}{E_{Q}^{3/2}k_{Q}}/{\left(\rho
V\right)^{1/2}}$. We expect $T_{bb}\left(Q,K\right)/N_{{\rm inertial}}$ in
Fig. 4b to be approximately independent of $Q$ in the inertial range.

In the tearing-mediated range, the energy cascade time is governed by
$1/\gamma_{p}$, where $\gamma_{p}$ is the linear growth rate of the tearing
instability. For a current sheet with a half thickness $\lambda$ and an
upstream Alfvén speed $V_{A\lambda}$, the linear growth rate [42]
$\gamma_{p}\sim(V_{A\lambda}\lambda/\eta)^{-1/2}V_{A\lambda}/\lambda$. For the
$Q$-th shell, $\lambda\sim 1/k_{Q}$, $V_{A\lambda}\sim
E_{Q}^{1/2}/(V\rho)^{1/2}$. Hence, we propose another normalization for the
reconnection-driven (or tearing-mediated) sub-inertial range, $N_{{\rm
tearing}}\sim\gamma_{p}E_{Q}\sim{\eta^{1/2}E_{Q}^{5/4}k_{Q}^{3/2}/\left(\rho
V\right)^{1/4}}$. We expect $T_{bb}\left(Q,K\right)/N_{{\rm tearing}}$ in Fig.
4c to be nearly independent of $Q$ in the tearing-mediated sub-inertial range.

In our simulation, the plasma density $\rho\simeq 1$ and the volume $V=2$. We
adopt the value $\rho V=2$ in the normalization for producing Fig. 4.

Supplementary Materials for
Reconnection-Driven Energy Cascade in Magnetohydrodynamic Turbulence

This PDF file includes:
Fig. S1. Current sheets in low-$\boldsymbol{R_{m}}$ 3D versus 2D MHD
turbulence.
Fig. S2. Cell-averaged simulation results and comparison to spectral
observations.
Section S1. Comparative simulations with a low magnetic Reynolds number.
Section S2. Exploring the connections to solar observations.
Movie S1. Animation of current density $|J|$.

## Section S1. Comparative simulations with a low magnetic Reynolds number.

For both 3D and 2D studies, we ran additional simulations with relatively low
magnetic Reynolds number, $R_{m}=8\times 10^{4}$ (or $R_{m,{\rm
eff}}=1.6\times 10^{4}$), at resolutions of $\sim 2000\times 2000\times 1000$
and $2000\times 2000$ grid cells, respectively. We set the magnetic
diffusivity $\eta=1.25\times 10^{-5}$ and viscosity $\nu=1.25\times 10^{-5}$,
such that the magnetic Prandtl number $P_{m}=\nu/\eta=1$. These two additional
simulations are used to comparatively demonstrate that only for large $R_{m}$
should magnetohydrodynamic turbulence be able to show the copious formation of
magnetic flux ropes (or plasmoids). Fig. S1 depicts the electric current
density component, $J_{z}$, in a 2D cut (at $z=-1$) of the low-$R_{m}$ 3D
simulation and in the low-$R_{m}$ 2D simulation. Different from Fig. 1d in the
main narrative, these low-$R_{m}$ simulations do not produce a significant
amount of magnetic flux ropes. This is because, for the low $R_{m}$, the
evolving current sheets cannot reach sufficiently fine scales to trigger the
tearing instability [15, 16, 17, 18, 19] before being dissipated at scales
$\sim 1/k_{\eta}$.

Figure 7: Current sheets in low-$\boldsymbol{R_{m}}$ 3D versus 2D MHD
turbulence. This figure (with low-$R_{m}$) is analogous to Fig. 1d in the main
narrative from large-$R_{m}$ simulations. Out-of-plane current density $J_{z}$
in a $x$-$y$ slice (at $z=-1$) of the 3D low-$R_{m}$ turbulence simulation
(left) compared with the corresponding result from the 2D low-$R_{m}$
simulation (right). The zoomed-in subdomains coincide with those in Fig. 1d.
Compared to Fig. 1d, it is clear that these low-$R_{m}$ simulations do not
produce the copious formation of plasmoids.

## Section S2. Exploring the connections to solar observations.

We mimic different instrument spatial resolutions by averaging the original 3D
simulation dataset at different levels (Fig. S2a). In the zoomed-in subdomains
of Fig. S2b, we can only observe plasmoids in the first two panels while they
are no longer observable in the last panel with 25$\times$25 cell average.
Consistent with Fig. 2 in the main narrative, the original dataset in Fig. S2c
can capture the new range of energy transfer in the magnetic energy spectrum
with a spectral index of $-2.2$ as has been observed by NIRIS (Fig. S2d) [7].
While the spectrum associated with the 5$\times$5 cell averaged dataset can
partially exhibit this new range of energy transfer, the spectra due to
25$\times$25 and 125$\times$125 cell averaging are no longer able to reproduce
observations, consistent with the discrepancy between NIRIS and HMI
observations at different spatial resolution. Here we want to reiterate that
it requires further observations to rule out the coincidence of the spectral
break from NIRIS observations and then to identify the underlying mechanism
that leads to the spectral break.

Figure 8: Cell-averaged simulation results and comparison to spectral
observations. (a) 2D slices of $|\mathbf{J}|^{2}$ (proportional to Joule
heating, $\eta|\mathbf{J}|^{2}$, and therefore thermal emission) from the
original and cell-averaged 3D datasets with $R_{m}=10^{6}$ at $z=-1$. (b)
Zoomed-in subdomains that are highlighted in the white boxes of Panel (a). (c)
Magnetic energy spectra from the original and cell-averaged 3D datasets with
$R_{m}=10^{6}$. (d) The magnetic energy spectrum derived from the high-
resolution NIRIS magnetograms of the Big Bear Solar Observatory (BBSO)
exhibits a new range of energy transfer with a spectral index of $k^{-11/5}$.
In contrast, the spectra calculated from the 25$\times$25 and 125$\times$125
cell-averaged 3D datasets do not show this new spectral range.

Movie S1. Animation of current density $|J|$. The electric current density
$|J|$ in different x-y slices through the entire 3D simulation domain for the
high-$R_{m}$ case.

## Acknowledgments

We acknowledge fruitful discussions with Muni Zhou, Silvio Cerri, Nan Liu,
Manasvi Lingam, Yuxi Chen, Greg Hammett, Steven Cowley, and Adam Burrows.
Resources supporting this work were provided by the NASA High-End Computing
(HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames
Research Center. We also would like to acknowledge high-performance computing
support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s CISL,
sponsored by NSF, and from National Energy Research Scientific Computing
Center, a DOE Office of Science user facility. Funding: This work was
partially supported by NASA grants 80NSSC19K0621 and 80NSSC21K1326, DOE grants
DE-SC0021205 and DE-SC0021254, and DOE under contract number DE-
AC02-09CH11466. Author contributions: All of the authors made significant
contributions to this work. C.D. and Y.M.H. designed the simulation for this
study; C.D. carried out the simulations of MHD turbulence; C.D., L.W., Y.M.H.,
L.C., and T.A.S. analyzed the simulation results; and C.D., L.W., Y.M.H.,
L.C., and A.B. wrote the paper; and all authors read the paper. Competing
interests: The authors declare no competing interests. Data and materials
availability: All data needed to evaluate the conclusions in the paper are
present in the paper and/or the Supplementary Materials.

## References

  * [1] P. S. Iroshnikov, Turbulence of a Conducting Fluid in a Strong Magnetic Field.  Astronomicheskii Zhurnal 40, 742 (1963).
  * [2] R. H. Kraichnan, Inertial-range spectrum of hydromagnetic turbulence.  Physics of Fluids 8, 1385 (1965).
  * [3] P. Goldreich, S. Sridhar, Toward a theory of interstellar turbulence. 2: Strong alfvenic turbulence.  The Astrophysical Journal 438, 763 (1995).
  * [4] W. H. Matthaeus, S. L. Lamkin, Turbulent magnetic reconnection.  Physics of Fluids 29, 2513 (1986).
  * [5] D. Biskamp, H. Welter, Dynamics of decaying two-dimensional magnetohydrodynamic turbulence.  Physics of Fluids B: Plasma Physics 1, 1964–1979 (1989).
  * [6] K. H. Kiyani, K. T. Osman, S. C. Chapman, Dissipation and heating in solar wind turbulence: from the macro to the micro and back again.  Philosophical Transactions of the Royal Society of London Series A 373, 20140155-20140155 (2015).
  * [7] V. I. Abramenko, V. B. Yurchyshyn, Analysis of quiet-sun turbulence on the basis of SDO/HMI and goode solar telescope data.  Monthly Notices of the Royal Astronomical Society 497, 5405-5412 (2020).
  * [8] R. Kieokaew, B. Lavraud, Y. Yang, W. H. Matthaeus, D. Ruffolo, J. E. Stawarz, S. Aizawa, C. Foullon, V. Génot, R. F. Pinto, N. Fargette, P. Louarn, A. Rouillard, A. Fedorov, E. Penou, C. J. Owen, T. S. Horbury, H. O’Brien, V. Evans, V. Angelini, Solar Orbiter observations of the Kelvin-Helmholtz waves in the solar wind.  Astronomy & Astrophysics 656, A12 (2021).
  * [9] N. F. Loureiro, S. Boldyrev, Role of magnetic reconnection in magnetohydrodynamic turbulence.  Physical Review Letters 118 (2017).
  * [10] A. Mallet, A. A. Schekochihin, B. D. G. Chandran, Disruption of sheet-like structures in alfvénic turbulence by magnetic reconnection.  Monthly Notices of the Royal Astronomical Society 468, 4862–4871 (2017).
  * [11] S. Boldyrev, N. F. Loureiro, Magnetohydrodynamic turbulence mediated by reconnection.  The Astrophysical Journal 844, 125 (2017).
  * [12] L. Comisso, Y.-M. Huang, M. Lingam, E. Hirvijoki, A. Bhattacharjee, Magnetohydrodynamic turbulence in the plasmoid-mediated regime.  The Astrophysical Journal 854, 103 (2018).
  * [13] N. F. Loureiro, S. Boldyrev, Nonlinear Reconnection in Magnetized Turbulence.  The Astrophysical Journal 890, 55 (2020).
  * [14] The first analytic study of the role of magnetic reconnection in MHD turbulence is attributable to Ref.[43]. This work, however, was based on an incorrect model of MHD turbulence [44].
  * [15] N. F. Loureiro, A. A. Schekochihin, S. C. Cowley, Instability of current sheets and formation of plasmoid chains.  Physics of Plasmas 14, 100703-100703 (2007).
  * [16] F. Pucci, M. Velli, Reconnection of Quasi-singular Current Sheets: The “Ideal” Tearing Mode.  The Astrophysical Journal Letters 780, L19 (2014).
  * [17] L. Comisso, M. Lingam, Y. M. Huang, A. Bhattacharjee, General theory of the plasmoid instability.  Physics of Plasmas 23, 100702 (2016).
  * [18] D. A. Uzdensky, N. F. Loureiro, Magnetic Reconnection Onset via Disruption of a Forming Current Sheet by the Tearing Instability.  Physical Review Letters 116, 105003 (2016).
  * [19] Y.-M. Huang, L. Comisso, A. Bhattacharjee, Plasmoid Instability in Evolving Current Sheets and Onset of Fast Reconnection.  The Astrophysical Journal 849, 75 (2017).
  * [20] A. Lazarian, G. L. Eyink, A. Jafari, G. Kowal, H. Li, S. Xu, E. T. Vishniac, 3d turbulent reconnection: Theory, tests, and astrophysical implications.  Physics of Plasmas 27, 012305 (2020).
  * [21] A. A. Schekochihin, MHD Turbulence: A Biased Review.  arXiv e-prints p. arXiv:2010.00699 (2020).
  * [22] D. Biskamp, Magnetohydrodynamic Turbulence (Cambridge University Press, 2003).
  * [23] W. Daughton, V. Roytershteyn, H. Karimabadi, L. Yin, B. J. Albright, B. Bergen, K. J. Bowers, Role of electron physics in the development of turbulent magnetic reconnection in collisionless plasmas.  Nature Physics 7, 539-542 (2011).
  * [24] Y.-M. Huang, A. Bhattacharjee, Turbulent Magnetohydrodynamic Reconnection Mediated by the Plasmoid Instability.  The Astrophysical Journal 818, 20 (2016).
  * [25] C. Dong, L. Wang, Y.-M. Huang, L. Comisso, A. Bhattacharjee, Role of the plasmoid instability in magnetohydrodynamic turbulence.  Physical Review Letters 121 (2018).
  * [26] T. Li, E. Priest, R. Guo, Three-dimensional magnetic reconnection in astrophysical plasmas.  Proceedings of the Royal Society A 477, 20200949 (2021).
  * [27] O. Pezzi, F. Pecora, J. Le Roux, N. Engelbrecht, A. Greco, S. Servidio, H. Malova, O. Khabarova, O. Malandraki, R. Bruno, et al., Current sheets, plasmoids and flux ropes in the heliosphere.  Space Science Reviews 217, 1–70 (2021).
  * [28] S. Boldyrev, Spectrum of magnetohydrodynamic turbulence.  Physical Review Letters 96 (2006).
  * [29] Although Ref. [43] predicted the same critical scale, it was derived based on an incorrect model of MHD turbulence[44].
  * [30] D. Berghmans, F. Auchère, D. M. Long, E. Soubrié, M. Mierla, A. N. Zhukov, U. Schühle, P. Antolin, L. Harra, S. Parenti, O. Podladchikova, R. Aznar Cuadrado, É. Buchlin, L. Dolla, C. Verbeeck, S. Gissot, L. Teriaca, M. Haberreiter, A. C. Katsiyannis, L. Rodriguez, E. Kraaikamp, P. J. Smith, K. Stegen, P. Rochus, J. P. Halain, L. Jacques, W. T. Thompson, B. Inhester, Extreme-UV quiet Sun brightenings observed by the Solar Orbiter/EUI.  Astronomy & Astrophysics 656, L4 (2021).
  * [31] Y. Chen, D. Przybylski, H. Peter, H. Tian, F. Auchère, D. Berghmans, Transient small-scale brightenings in the quiet solar corona: A model for campfires observed with Solar Orbiter.  Astronomy & Astrophysics 656, L7 (2021).
  * [32] D. Clery, ‘Campfires’ may drive heating of solar atmosphere.  Science 372, 557-558 (2021).
  * [33] J. F. Drake, M. Swisdak, H. Che, M. A. Shay, Electron acceleration from contracting magnetic islands during reconnection.  Nature 443, 553-556 (2006).
  * [34] L. Comisso, L. Sironi, Ion and electron acceleration in fully kinetic plasma turbulence.  The Astrophysical Journal Letters 936, L27 (2022).
  * [35] S. Kong, V. Ossenkopf-Okada, H. G. Arce, J. Bally, Á. Sánchez-Monge, P. McGehee, S. Suri, R. S. Klessen, J. M. Carpenter, D. C. Lis, F. Nakamura, P. Schilke, R. J. Smith, S. Mairs, A. Goodman, M. J. Maureira, The CARMA-NRO Orion Survey: Filament Formation via Collision-induced Magnetic Reconnection—the Stick in Orion A.  The Astrophysical Journal 906, 80 (2021).
  * [36] R. Narayan, The Physics of Pulsar Scintillation.  Philosophical Transactions of the Royal Society of London Series A 341, 151-165 (1992).
  * [37] G. Tóth, B. van der Holst, I. V. Sokolov, D. L. De Zeeuw, T. I. Gombosi, F. Fang, W. B. Manchester, X. Meng, D. Najib, K. G. Powell, Q. F. Stout, A. Glocer, Y.-J. Ma, M. Opher, Adaptive numerical algorithms in space weather modeling.  Journal of Computational Physics 231, 870-903 (2012).
  * [38] Y. Chen, G. Tóth, T. I. Gombosi, A fifth-order finite difference scheme for hyperbolic equations on block-adaptive curvilinear grids.  Journal of Computational Physics 305, 604-621 (2016).
  * [39] A. Alexakis, B. Bigot, H. Politano, S. Galtier, Anisotropic fluxes and nonlocal interactions in magnetohydrodynamic turbulence.  Phys. Rev. E 76, 056313 (2007).
  * [40] A. Alexakis, P. D. Mininni, A. Pouquet, Shell-to-shell energy transfer in magnetohydrodynamics. I. Steady state turbulence.  Phys. Rev. E 72, 1–10 (2005).
  * [41] O. Debliquy, M. K. Verma, D. Carati, Energy fluxes and shell-to-shell transfers in three-dimensional decaying magnetohydrodynamic turbulence.  Phys. Plasmas 12, 1–10 (2005).
  * [42] D. Biskamp, Magnetic Reconnection in Plasmas, vol. 3 (Cambridge University Press, 2000).
  * [43] V. Carbone, P. Veltri, A. Mangeney, Coherent structure formation and magnetic field line reconnection in magnetohydrodynamic turbulence.  Physics of Fluids A: Fluid Dynamics 2, 1487–1496 (1990).
  * [44] S. Boldyrev, N. F. Loureiro, Tearing Instability in Alfvén and Kinetic-Alfvén Turbulence.  Journal of Geophysical Research (Space Physics) 125, e28185 (2020).