Laure Berti-Équille

Ugo Comignani
email: ugo.comignani@lis-lab.fr

Explaining Automated Data Cleaning with CLeanEX

In this paper, we study the explainability of automated data cleaning pipelines and propose CLeanEX, a solution that can generate explanations for the pipelines automatically selected by an automated cleaning system, given it can provide its corresponding cleaning pipeline search space. We propose meaningful explanatory features that are used to describe the pipelines and generate predicate-based explanation rules. We compute quality indicators for these explanations and propose a multi-objective optimization algorithm to select the optimal set of explanations for user-defined objectives. Preliminary experiments show the need for multi-objective optimization for the generation of high-quality explanations that can be either intrinsic to the single selected cleaning pipeline or relative to the other pipelines that were not selected by the automated cleaning system. We also show that CLeanEX is a promising step towards generating automatically insightful explanations, while catering to the needs of the user alike.

Introduction

When it comes to real-world data, inaccurate, noisy, uncertain, or incomplete data are the norm rather than the exception. In many applications of machine learning (ML), the cost of an error can be high. However, explaining the results of an ML pipeline is as crucial as reducing the impact of dirty data input and estimating the uncertainty in the predictions of a model. Data scientists have to spend considerable time integrating data from multiple sources, manually curating and preparing the data with their expertise of the domain. They use various libraries and tools to correct erroneous values, impute missing ones, eliminate duplicate records or disambiguate conflicting data, to prevent unreliable data being delivered downstream to a machinelearning model. Although data cleaning in ML pipelines is still considered to be "intractable" or "AI-hard" -as it is difficult to fully automate and often requires human expertiseseveral solutions of automated data cleaning have been proposed recently [START_REF] Shang | [END_REF]Krishnan et al., 2017; . While important, however, we argue that a critical missing piece of this line of work be the explainability of automated cleaning methods, i.e., why/how a particular data value is detected as an anomaly, how it is fixed, and why it is fixed the way it is fixed.

Current automated cleaning methods hardly ever explain their choice in building automatically data cleaning pipelines. However, the lack of convincing explanations may severely limit the scope of an end-to-end ML pipeline (including data cleaning and data preparation) and the broader adoption of ML-based Artificial Intelligence for critical decision-making. In this paper, therefore, we study the explainability of automated data cleaning pipelines and propose CLeanEX, a solution that can generate explanations for a pipeline generated by an automated data cleaning system that can provide its pipeline search space. CLeanEX extends the automated cleaning system Learn2Clean1 we proposed in [Berti-Equille, 2019a;Berti-Équille, 2019b], based on Q-Learning, a model-free reinforcement learning technique adapted for automating data cleaning and data preparation.

The paper proceeds as follows. In Section 2, we discuss related work. Section 3 presents our definitions and describe how to generate explanations from a sequence of cleaning tasks. Experiments are reported in Section 4. Section 5 concludes the paper.

As illustrated in Fig. 1, data cleaning and preparation require a sophisticated sequence of tasks T i for the detection and elimination of a variety of intricate data quality problems (e.g., duplicates, inconsistent, missing, and outlying values). Generally, the users may not know which preprocessing methods can be applied to optimize the final results downstream. The selection of the optimal sequence of tasks would require for him/her executing all possible methods for each potentially relevant preprocessing task, as well as all the possible combinations of the methods with different orderings and configurations (see all the alternative pipelines represented in grey in the figure). Successfully automating this daunting manual process would definitively improve data scientists' every day life. Related to our work, we review recent advances for scalable and automated data cleaning and explainable systems.

ML-based data and pipeline curation. A recent line of work is to use machine learning to improve the efficiency and reliability of data cleaning and data repairing [START_REF] Yakout | [END_REF]Rekatsinas et al., 2017;Krishnan et al., 2017;[START_REF] Shang | [END_REF]. [START_REF] Yakout | [END_REF] train ML models and evaluate the likelihood of recommended replacement values to fix erroneous and missing ones. [Rekatsinas et al., 2017] propose a framework for data repairing based on probabilistic inference to handle integrity constraints or external data sources seamlessly, with quantitative and statistical data repairing methods. [Krishnan et al., 2017] propose a method for data cleaning optimization for user-defined ML-based analytics. Their approach selects an ensemble of methods (statistical and logic rules) and boosting to combine error detection and repair. AutoML approaches can optimize the hyperparameters of a considered ML model, but they support only a limited number of preprocessing steps with by-default methods. Recently, Alpine Meadow [START_REF] Shang | [END_REF] combines an AutoML approach and a cost model to select candidate logical ML pipeline plans (as in DB query optimization). Multi-armed bandits are used to select promising logical ML pipeline plans, and Bayesian Optimization is used to fine-tune the hyper-parameters of the selected models in the search space. However, these approaches do not provide explanations justifying how the data preparation pipeline search space is built and pruned. We argue that more efforts should be devoted to proposing a principled and efficient explainable automated data cleaning approach to help the user in understanding the sequence of data preparation tasks selected by the automated cleaning system.

Explainable systems for data quality and data profiling. The explanation of predictive models has been an overgrowing field in the last years [START_REF] Guidotti | [END_REF][START_REF] Pedreschi | [END_REF]. Despite the numerous works on black box models explanation, few of them address the problem of explaining the choices performed during data cleaning tasks [START_REF] Bertossi | [END_REF]. [START_REF] Rammelaere | [END_REF] proposes to explain repairs performed by the user over small data sets by extracting the most relevant conditional functional dependencies (CFDs) according to the repairs. However, while this approach can be adapted to explain the output of an automated data cleaning task, the CFDs can be hard to understand for the user and they have a limited expressiveness, which might be detrimental to explain the result of non rule-based repairing tasks. The T-REx framework [Deutch et al., 2020] provides explanations taking the form of a ranking of the set of Denial Constraints (DCs) according to their influence during the repairing process. This approach has to provide the same set of DCs that is used as input of the data repairing task. In the context of Entity Resolution (ER) [START_REF] Ebaid | [END_REF], local explanations have been proposed to explain why a tuple pair was predicted to be a duplicate. In this context, other approaches include identifying a ranked list of features that contribute to the prediction as a duplicate (such as LIME [Ribeiro et al., 2016]) or as a rule that holds in the vicinity of the tuple-pair with high precision (such as Anchor [START_REF] Ribeiro | [END_REF]). While the current approaches explain the effect of a single repairing task, we choose to address the problem of explaining the full sequence of cleaning tasks (the cleaning pipeline). As this sequence has been selected and generated by an automated data cleaning system, we argue that the explanations should also offer the perspective related to the other candidate pipelines that were not selected by the system using various indicators to characterize the quality of an explanation.

Our Definitions

In this work, we represent the pipeline exploration space of the automated cleaning agent A as a tree. Each cleaning task updates the input data set and the updated version is represented as a node in a tree structure: the cleaning tree. As illustrated in Fig. 1, the cleaning tree is composed of branches representing alternative cleaning pipelines among which the cleaning agent A may select the optimal one with respect to an ML goal and a quality performance metric. Each pipeline has features that can be used and exposed to explain the agent decision independently of the ML quality performance metric optimization. Explanation rules can be generated using such features. More formally, we define the cleaning tree as follows. Definition 1 (Cleaning tree). Let D be a data set, T = {t 1 , . . . , t n } be the set of data preprocessing tasks that can be executed by the automated cleaning agent A using an ML model M that maximizes the model quality performance metric q. A cleaning tree for A, M , and q over D is an oriented tree C A,M,q = (V, v r , E) such that: • V is the set nodes formed by triples v p = (D p , q(M, D p ), s p ) with: -D p , a data set resulting from the application of a cleaning task t p to the parent node data set or, for the root node, the input data set D; q(M, D p ), the value of the quality metric of the ML model M measured over the pre-processed data set D p ; s p a vector associating, for each explanatory feature f ∈ F, its value s p [f ] for the current node; • v r the root node of the cleaning tree C A,M,q such that v r ∈ V and v r = (D, q(M, D), s r ); • E is the set of edges in C A,M,q , with each edge e from a node v = (D, q(M, D), s) to one of its children v = (D , q(M, D ), s ) being labelled with the task used to clean and update D as D . We define the bijection φ : E → V × V which maps each edge e ∈ E to the pair of nodes (v, v ) connected by e. Definition 2 (Cleaning pipeline). Given a cleaning tree C A,M,q = (V, v 0 , E) and a node v n ∈ V. A cleaning pipeline p vn over C A,M,q is a path σ E = e 1 , . . . , e n of edges in E such that there exists a sequence of nodes from

V, σ V = v 0 , v 1 , . . . , v n-1 , v n such that each edges e i ∈ σ E connects the node v i-1 to the node v i ∈ σ V : ∀e i ∈ σ E , φ(e i ) = (v i-1 , v i ).
After the automated cleaning agent A selects the node with the best quality metric value, the path from the root to this selected node corresponds to the best pipeline p opt(M,q) that is the optimal sequence of tasks evaluated by A for the ML model M and quality metric q.

To ground our explanations, we use explanatory features: the values of the explanatory features changing from one node to the next in a cleaning pipeline will serve as key elements for generating explanations of an automatically selected cleaning strategy. For each pipeline, we compute the following explanatory features, denoted s p [f ] in Definition 1:

• cost: The normalized cost of the pipeline;

• data quality improvement: The percentage of data quality problems solved by the pipeline (e.g., remove 100% of missing values by imputation). Note that quantifiable data quality indicators have to be defined beforehand to characterize the main dimensions of data quality (i.e., consistency, accuracy, etc. Please see [Comignani et al., 2020;Berti-Équille, 2007] for more detail to specify and compute data quality indicators); • distortion: The distortion as defined in [START_REF] Dasu | [END_REF] as the Mahalanobis distance between the original and cleaned version of the data set; • satisfaction: The satisfaction of ML model requirements by the pipeline defined as a Boolean: e.g., for regression, satisfaction equals 1 if linearity, multivariate normality, no or little multicollinearity, no auto-correlation, and homoscedasticity constraints are satisfied by the cleaned data set; • corr ratio: The fraction of the number of pipelines sharing the same tasks over the sum of their respective ranks and the total number of explored pipelines; and • non corr ratio: The fraction of the number of pipelines that do not share the same task over the sum of their respective ranks and the total number of explored pipelines. Other features could be added but the ones we consider define a representative and relevant signature that can be used for generating meaningful explanations of a cleaning pipeline independently from the ML model and the model quality metric used by the automated cleaning system. Example 1. We illustrate the notations with an example using the House Price data set from Kaggle2 with 81 variables and 1.46k observations to predict the SalePrice attribute with regression and MSE as the quality metric. We

Explaining via Feature Changes and Comparisons

Next, we define two types of explanation predicates: (1) Change predicates that characterize the evolution of the explanatory features, and (2) Comparison predicates that characterize the relative changes by comparing with other pipelines in the cleaning tree. Change Predicates.

Given a cleaning tree C A,M,q = (V, v r , E), a set of features F, two distinct nodes v, v ∈ V and s, s their respective vectors of features values. We define the following set of change predicates:

• succ(v, v ) iff there exists an oriented path from v to v ; • increase(f, v, v ) iff: s[f ] < s [f ] ∧ succ(v, v ) • decrease(f, v, v ) iff: s[f ] > s [f ] ∧ succ(v, v ) • stable([f 0 , . . . , f n ], v, v ) iff: ∀f i ∈ [f 0 , . . . , f n ], s[f i ] = s [f i ] ∧ succ(v, v ) • equiv(v, v ) iff: ∀f ∈ F, s[f ] = s [f ].
Comparison Predicates. We define the following set of comparison predicates to express that a value of feature f for node v is respectively: (i) more, less than, or as a% of the feature f values of the other nodes in the cleaning tree with a, a decimal constant in [0,1]; (ii) the best value, the worst value, or (iii) different from another node v feature value within a certain distance d in [0,1]: Figure 2: Example of a selected cleaning tree (in green) in the pipeline space of our illustrative example (a), with its corresponding explanatory features (b), and explanation tree (c) composed of the sets of rule IDs described in Table 2.

• more([f 0 , . . . , f n ], v, a) iff: ∀f i ∈ [f 0 , . . . , f n ], a = |{v ∈V|s[fi]>s [fi]}| |V|-1 • less([f 0 , . . . , f n ], v, a) iff: ∀f i ∈ [f 0 , . . . , f n ], a = |{v ∈V|s[fi]<s [fi]}| |V|-1 • as([f 0 , . . . , f n ], v, a) iff: ∀f i ∈ [f 0 , . . . , f n ], a = |{v ∈V|s[fi]=s [fi]}| |V|-1 • most([f 0 , . . . , f n ], v) iff: ∀f i ∈ [f 0 , . . . , f n ], ∀v ∈ V s.t. v = v , s[f i ] > s [f i ] • least([f 0 , . . . , f n ], v) iff: ∀f i ∈ [f 0 , . . . , f n ], ∀v ∈ V s.t. v = v , s[f i ] < s [f i ] • δ(f, v, v , a) iff: |s [f ] -s[f ]| = d
Additionally to the axioms of the predicate logic, we define the following set of axioms for our explanation language:

• ∀p 1 , p 2 ∈ {more, as, less} s.t.

p 1 = p 2 : p 1 (f, v, 1) ≡ p 1 (f, v, 1) ∧ p 2 (f, v, 0) • more(f, v, a) ∧ less(f, v, a ) ≡ more(f, v, a) ∧ as(f, v, 1 -(a + a )) ≡ as(f, v, 1 -(a + a )) ∧ less(f, v, a ) • most([f 0 , . . . , f n ], v) ≡ i∈[0,n] more(f i , v, 1) • least([f 0 , . . . , f n ], v) ≡ i∈[0,n] less(f i , v, 1) • equiv(v 1 , v 2 ) ⇔ f ∈F stable(f, v 1 , v 2 )
Now, we can define formally the explanation rules that will be automatically generated.

Definition 3 (Explanation Rule). Given a cleaning tree C A,M,q = (V, v r , E), a feature f in F, v and v two nodes in V, three decimals constants, a, a , a ∈ [0, 1], and the explanation predicates: p 0 ∈ {increase; stable; decrease}; p 1 , p 2 ∈ {more; as; less} such that p 1 = p 2 ; p 3 ∈ {most; least}. We define the following types of explanation rules: Type P explanation rules are path explanations between two nodes v and v . Type B explanation rules are behavior change explanations with a constant representing the δ change value of some feature f between two nodes v and v . Type C explanations are comparative explanations with respect to the change of other nodes. Type E explanation rules are related to extreme value description. Type S explanation rules are stability explanation rules between v and v .

P v,v : ∃v 1 , . . . , v n ∈ V, succ(v, v 1 ) ∧ . . . ∧ succ(v n , v ) B f,v,v : p 0 (f, v, v ) ∧ δ(f, v, v , a) C f,v : p 1 (f, v, a ) ∧ p 2 (f, v, a ) E f,v : p 3 (f, v) S f,v,v : stable(f, v, v ) S v,v : equiv(v, v )
These explanation rules can be either focused on one feature f for the S f,v,v explanation rules, or all the features covered by C A,M,q for the S v,v explanation rules. Table 2 gives the explanation rules automatically generated for the cleaning pipeline KN N → ZSB. Finally, all the generated rules can be organized into an explanation tree mirroring the cleaning tree defined earlier, as follows.

Definition 4 (Explanation Tree). Let C A,M,q = (V, v r , E) be a cleaning tree and R be a set of explanation rules such that C A,M,q |= R. An explanation tree for R in C A,M,q is an oriented tree X C = (V, E) such that:

• V is the set nodes such that each node v Ri is formed by the subset of rules R i ⊆ R such that each r ∈ R i is a rule of type C f,vi or E f,vi ; • E is the set of edges in X C , with each edge e from a node v Ri to one of its children v Rj being labelled with the set of rules R i,j such that each r ∈ R i,j is a rule of type B f,vi,vj , S f,vi,vj or S vi,vj .

The branch of the explanation tree corresponding to the selected pipeline of the example is provided in Fig. 2(c).

Quality Indicators of Explanations

Given an explanation tree X C , its set of predicates, denoted P C including the set of predicates P ml of the form most([f 0 , . . . , f n ], v) or least([f 0 , . . . , f n ], v), f eatures(p), the function returning the set of features covered by a predicate p, and delta(r) the function returning the value a of the delta predicate δ(f, v, v , a) in a rule r. We select the best explanations to provide to the user using the following relevant quality indicators defined as follows. Definition 5 (Polarity). The polarity of X C is defined as:

polarity(X C ) = ∀p∈P ml |f eatures(p)| ∀p ∈P C |f eatures(p )| .
Definition 6 (Distancing). Given the set X of explanation trees over C A,M,q such that X C ∈ X. The distancing of X C is defined as:

distancing(X C ) = ∀ delta predicate r in X C delta(r) ∀ delta predicate r in X delta(r)
.

Definition 7 (Surprise). Given a set X of explanation trees over C A,M,q such that X C ∈ X. The surprise of X C is defined as:

surprise(X C ) = ∀X ∈X ∀ rule r of type B in X delta(r) ∀predicate p in X f eatures(p)
.

Definition 8 (Diversity). Given S the set of all predicates symbols excepting δ in the rules of X C (defined in Def. 3 for P,B,C,E and S types of explanations), and P s the set of all predicates with predicate symbol s in the rules of X C .The diversity of X C is defined as:

diversity(X C ) = -

Multi-Objective Optimization

Next, we use these quality indicators as optimization objectives to select the optimal explanations of the automated choice of a cleaning pipeline. When dealing with more than one dimension to optimize, there may be many incomparable sets of explanation. Therefore, we adapt the notion of explanation plan and Pareto plan as follows. Definition 9 (Explanation Plan). Explanation plan π i , associated to x i , a branch of an explanation tree X C , is a tuple π i = polarity(x i ), distancing(x i ), surprise(x i ), diversity(x i ) . Definition 10 (Explanation Sub-Plan). Explanation plan π i is the sub-plan of another plan π j if their associated explanation sub-trees satisfy x i ⊆ x j .

Algorithm 1 Approximated Pareto-optimal Explanations Input: The user-defined optimization objective O, the precision value α, the maximum number of rules of an explanation k, the explanation tree X , and the pipeline to explain x Output: The set of explanations Π α for x 1 Π α ←-π x 2 for all explanation branches B of the tree X \ x do Definition 11 (Dominance). Plan π 1 dominates π 2 if π 1 has better or equivalent values than π 2 in every quality indicator. The term better is equivalent to greater for maximization objectives (e.g., diversity or polarity), and lower form in minimization ones (e.g., distancing or surprise) depending on the user's preferences. Furthermore, plan π 1 strictly dominates π 2 if π 1 dominates π 2 and the values of indicators for π 1 and π 2 are not equal. Definition 12 (Pareto Plan). Plan π i is Pareto if no other plan strictly dominates π i . The set of all Pareto plans is denoted as Π.

Now, we define our problem as a multi-objective optimization (MOO) problem as follows: for an explanation tree X C , the problem is to find all the explanation branches x, such that each branch satisfies:

• diversity(x) is maximized; • surprise(x) is optimized; • polarity(x) is optimized; • distancing(x) is optimized;
Note that while we always maximize diversity, we may either minimize or maximize the distancing, surprise, and polarity based on the user's needs. The main challenge in designing an algorithm for finding optimal explanations of automated data cleaning is the multi-objective nature of the problem. A multi-objective problem can be easily solved if it is possible to combine all quality indicators into one or if the optimization of one indicator leads an optimized value of other indicators. However, in our problem, the indicators may be conflicting, i.e., optimizing one does not necessarily lead to an optimized value for others and they cannot be combined into one single indicator. Therefore, we propose Algorithm 1 based on approximation to select the Pareto-optimal explanation rules. The algorithm makes less enumerations with a theoretical guarantee on the quality of results. Its pruning mechanism uses the precision value α. In the special case of α = 1, the algorithm operates exhaustively. If α > 1, the algorithm prunes more and hence is faster. In the latter case, a new plan is only compared with all plans that generate the same result. But a new plan are only inserted into the buffer if no other The algorithm 1 exploits a dynamic programming approach. The algorithm begins by constructing a plan for each single branch of an explanation tree (lines 2 to 3). It keeps all non α-dominated plans in a buffer. Then, it examines all the branches of the explanation tree with depth 1 up to k (lines 6 to 10). After each iteration, it removes α-dominated plans from the buffer. Finally, it returns the buffer content as result.

Preliminary Experiments

Based on our observations from the automated cleaning pipeline search spaces, we generated two types of cleaning tree structures with 10 and 100 nodes respectively and a maximal depth varying from 1 to 5 levels from root to leaves. The number of nodes per level follows a Beta distribution (alpha=2, beta=2.35). The values of the six explanatory features were generated randomly five times. We ran five trials of picking one cleaning pipeline to explain among all the candidate ones for each of the 250 configurations (2 × 5 × 5 × 5) and we averaged the results. We perform all experiments on a laptop Dell XPS machine with an Intel Core i7-7500U quad-core, 2.8 GHz, 16 GB RAM, powered by Windows 10 64-bit with Python. CLeanEX, data sets, and more detailed experiments are available at https://github.com/ucomignani/cleanex. Fig. 3 illustrates different Type C explanations for explaining the same cleaning pipeline. Each set of explanations has been selected as it optimizes specific objectives in terms of quality indicators. We can observe that in general, no correlation exists between the optimized value of the objectives and the explanation set can be very different. Thus, each objective should be optimized independently and according to the user's preferences. Fig. 4 shows the result of our algorithm (with α = 1 exhaustive search) for the top-10 ranked explanations depending on the optimization objectives. The ranking depends on the explanation size and the maximal depth of the tree: 2-3 rules for the explanations maximizing all objectives for depth from 2-3 (Max all). In Fig. 5, we can also examine the effect of the tree structure (in terms of the number of nodes and maximal depth) on execution time of CLeanEX to find the optimal explanations. We observe that the execution times for trees containing 10 and 100 nodes have similar values as their respective depth increases. The reason is that the number of subsets of candidate rules for a given branch increases quickly when the depth of the tree also increases. Inversely, reducing the depth of the branches increases the number of branches to explore, for which the depth generally decreases, reducing the time complexity for the cases with few nodes. The execution time of our algorithm is spent on the computation of the quality indicators, the generation of the rules, and the exploration of a large number of candidates rules sets. However, changing α will reduced this last dimension while providing an approximation guarantee.

Conclusions

In this paper, we propose CLeanEX, a framework to generate explanations for automated data cleaning. The key advantages of our approach is that our explanation system is:

(1) model-agnostic: it can be applied to any ML model, any set of data cleaning tasks, and any automated cleaning agent that can provide its cleaning pipeline search space; (2) logicbased: explanations are comprehensible to humans with diverse expertise and can be extensible to handle causal reasoning; (3) both local and global: it can explain both some part or the whole cleaning pipeline; (4) model quality-independent: it provides a reliable set of explanations independently from the ML model quality performance metric that is used by the automated cleaning agent for the selection of the optimal cleaning pipeline; and (5) user-defined: optimization is based on the user's objectives.

Figure 1 :

 1 Figure 1: CLeanEx Overview

  with n s , the number of nodes features over predicates with symbol s: n s = ∀p∈Ps |f eatures(p)|, and N , the total number of nodes features: N = ∀s∈S ∀p∈Ps |f eatures(p)|.

3π

  B ←-construct plan(B) 4 if π B is not α-dominated wrt O by any other plan in Π α 5 then Π α .add(π B ) 6 for d ∈ [1, k] do 7 for all explanation branches b of d rules in X \ x do 8 π b ←-construct plan(b) 9 if b is not α-dominated wrt O by any other plan in Π α 10 then Π α .add(π b ) 11 return Π α

C58:Figure 3 .Figure 5 .

 35 Figure 3. Optimal sets of explanations with their quality indicators depending on the multiple (left) single (right) objectives

Table 1 :

 1 Learn2Clean results on House Price data set for OLS regression over SalePrice attribute.

	Rank#	Sequence of Actions for Data Preparation	MSE	Time (ms)
	1	KN N → ZSB	2.20E-25	45.1
	2	M M → KN N → ZSB	1.44E-12	64.2
	3	IQR → M ICE → AD	0.231	154.9
	4	IQR → M ICE	0.236	94.1
	5	IQR → LC → M ICE	0.372	141.5
	ran Learn2Clean [Berti-Equille, 2019a] as the automated
	cleaning agent and reported the top-5 cleaning strategies pre-
	sented in Table 1 generated by Learn2Clean for OLS regres-
	sion. The best strategy has been selected based on minimal
	MSE. The MSE value, execution time (in milliseconds), and
	tasks of the explored cleaning pipelines are reported in the
	table. Data cleaning strategies include various tasks such
	as K-nearest neighbors (KNN) or MICE imputation, outlier
	detection/removal using InterQuartile Range (IQR), Z-score-
	based (ZSB), or LOF methods, approximate duplicate re-
	moval (AD), MinMax normalization (MM), and linearly cor-
	related (LC) feature exclusion (see [Berti-Equille, 2019a] for
	more details about the preprocessing methods). These strate-
	gies can be represented as the cleaning tree of Fig. 2(a) for
	which each node v

i has explanatory feature signature s i presented in the table of Fig.

2(b

). In the table, we observe cost= 0 for the root and cost= 1 for s 9 as the maximal cost for the leaf node of the pipeline IQR → M ICE → AD. The pipeline KN N → ZSB selected by Learn2Clean is represented in green. corr ratio of KN N in s 1 is 2 5+1+2 = 0.25 as it appears in 2 pipelines out of 5 with rank# 1 and 2 respectively (similarly for ZSB). no corr ratio of KN N is 3 5+3+4+5 = 0.176 as it does not appear in 3 pipelines ranked #3, 4, and 5 respectively.

Explanatory Features (a) Cleaning Tree root

	v4.2 (b) P : su cc(root, v1) ∧ su cc(v1,v4) ,0.625,s7 ) 0 0.36 0.16 0.20 0.66 0.52 0.60 0.51 0.82 1 0.91 c o s t s 1 s 2 s 3 s 4 s 5 s 6 s 7 s 8 s 9 s 10 S : equ iv(root, v1) C : δ(cost, root, v4,0.43) ∧ decrease(cost, root, v4) C : δ(dq_imp, root, v4,0.4) ∧ increase(dq_imp, root, v4) 0.42 0.3 0 0 0.15 0.42 0.3 0 0.25 0.176 0.11 0.72 0 0.167 0.21 0.89 0.72 0 0.176 0.25 0.82 0.47 0.2 0.25 0.176 0.05 0.66 0 0.25 0.176 0.97 0.42 0 0.176 0.25 0.88 0.1 0.8 0.10 0.26 0.63 0.72 0.2 0.25 0.176 0.91 0.85 0 0.125 0.24 0.34 0.97 1 0.176 0.25 d is to rt io n d q _ im p n o n _ c o rr _ ra ti o c o rr _ ra ti o s a ti s fa c ti o n C : δ(distortion, root, v4,0.17) ∧ increase(distortion, root, v4) C : δ(satisfaction, root, v4,1) ∧ increase(satisfaction, root, v4) C : δ(corr_ratio, root, v4,0.55) ∧ decrease(corr_ratio, root, v4) C : δ(non_corr_ratio, root, v4,0.05) ∧ increase(non_corr_ratio, root, v4) C : δ(cost, v1,v4,0.43) ∧ decrease(cost, v1,v4) C : δ(dq_imp, v1,v4,0.4) ∧ increase(dq_imp, v1,v4) C : δ(distortion, v1,v4,0.17) ∧ increase(distortion, v1,v4) C : δ(satisfaction, v1,v4,1) ∧ increase(satisfaction, v1,v4) C : δ(corr_ratio, v1,v4,0.55) ∧ decrease(corr_ratio, v1,v4) C : δ(non_corr_ratio, v1,v4,0.05) ∧ increase(non_corr_ratio, v1,v4)

(c

) Explanations for the selected cleaning pipeline KNN-> ZSB

	feature	R root = {1; 12; 13}
	signature #	R root,v1 = {3; 5; 6; 7}
		R v1 = {14; 15; 16; 17} R v1,v4 = {4; 8; 9; 10; 11}
		R v4 = {2; 18; 19; 20; 21}

Table 2 :

 2 Explanation rules generated for pipeline#1 KN N → ZSB

	Rule#	Type	Rule
	0	P	succ(root, n1) ∧ succ(n1, n4)
	1	E	least([cost, corr ratio, non corr ratio], root)
	2	E	most([cost, dq imp, distortion, satisfaction], n4)
	3	S	stable([dq imp, distortion, satisfaction], root, n1)
	4	S	stable([corr ratio, non corr ratio], n1, n4)
	5	B	increase(cost, root, v1) ∧ δ(cost, root, v1, 0.36)
	6	B	increase(corr ratio, root, v1) ∧ δ(corr ratio, root, v1, 0.25)
	7	B	increase(non corr ratio, root, v1)
			∧δ(non corr ratio, root, v1, 0.026)
	8	B	increase(cost,v1, v4) ∧ δ(cost, v1, v4, 0.3)
	9	B	increase(dq imp, v1, v4) ∧ δ(dq imp, v1, v4, 0.4)
	10	B	increase(distortion, v1, v4) ∧ δ(distortion, v1, v4, 0.17)
	11	B	increase(satisfaction, v1, v4) ∧ δ(satisfaction, v1, v4, 0.2)
	12	C	as([dq imp, distortion, satisfaction], root, 0.5)
	13		less([dq imp, distortion, satisfaction], root, 0.5)
	14	C	less(cost, v1, 0.5)
	15	C	more(cost, v1, 0.5)
	16	C	as([dq imp, distortion, satisfaction, corr ratio,
			non corr ratio], v1, 0.5)
	17	C	less([dq imp, distortion, satisfaction, corr ratio,
			non corr ratio], v1, 0.5)
	18	C	as(corr ratio, v4, 0.5)
	19	C	more(corr ratio, v4, 0.5)
	20	C	as(non corr ratio, v4, 0.5)
	21	C	more(non corr ratio, v4, 0.5)

https://github.com/LaureBerti/Learn2Clean

https://www.kaggle.com/datasets


P G Hughes

V Ábrahám

A Akimkin

L Banzatti

F Chen

Cruz-Sáenz

Jacob Aaron White

Á Kóspál

A G Hughes

P Ábrahám

V Akimkin

A Banzatti

L Chen

F Cruz-Sáenz De Miera

A Dutrey

M Flock

S Guilloteau

A S Hales

T Henning

K Kadam

D Semenov

A Sicilia-Aguilar

R Teague

ALMA and VLA Observations of EX

Keywords:

HAL is

INTRODUCTION

Corresponding author: J.A. White jwhite@nrao.edu Giant molecular clouds are the nurseries in which stars are born. The earliest phases of mass accumulation take place in the densest regions of the cloud cores and the rotation of these in-falling cores can lead to the formation of young stellar objects (YSOs) embedded in accretion disks.

A growing number of YSOs have been observed to have extreme outbursts (e.g., Audard et al. 2014), increasing the brightness by up to several magnitudes and the total luminosity by a factor of 10 -100. The prototypical outbursting YSO, FU Orionis, entered an outburst in 1936 (Herbig 1966) and is still slowly fading (Kenyon et al. 2000). This outburst has motivated the study of other similar systems (named FUors) that have also experienced outbursts of various timescales and intensities. In addition to the FUor class of young stars, there are systems that exhibit shorter duration and weaker outbursts called EXors (named after EX Lupi, discussed in detail below). Theoretical considerations (e.g., Zhu et al. 2009;D'Angelo & Spruit 2012;Vorobyov & Basu 2015) suggest that EXors/FUors, and Sun-like stars in general, can build up a significant fraction of their mass during periods of enhanced accretion, referred to as episodic accretion. The outbursts of both EXors and FUors are thought to be due to this episodic accretion of material from their circumstellar disks onto the protostars. While it is generally accepted that the observed outbursts are due to accretion, the exact triggering and transport mechanism(s) that deliver material from the disk on to the star over short time periods typical for FUors/EXors is still unknown (Audard et al. 2014).

Understanding the outburst mechanisms of protostars is important not only for building a complete picture of stellar evolution, but also for the potential implications for the planet formation process around low-mass stars. If these types of outbursts are a common byproduct of the star formation process then they will undoubtedly also impact disk evolution and therefore the conditions in which planets form. Outbursts have been shown to potentially change the chemistry and mineralogy of the surrounding circumstellar disk ( Ábrahám et al. 2009;Rab et al. 2017;Molyarova et al. 2018); could spur the growth of small solids through, e.g., evaporation and recondensation from a rapid evolution of the ice line (Cieza et al. 2016); and offer a potential solution to the luminosity problem (Hartmann & Kenyon 1996).

EX Lupi is a young 1 -3 Myr M0 star (0.6 M , 0.7 L quiescent luminosity) located 157.7 ± 0.9 pc away in the Lupus 3 cloud (Lombardi et al. 2008;Sipos et al. 2009;Frasca et al. 2017;Bailer-Jones et al. 2018). There have been multiple observed outbursts in the EX Lupi system, the most recent of which was in 2008. During this last outburst, EX Lupi was monitored in the optical and infrared. The outburst caused the annealing of amorphous silicate particles to crystalline grains, which were transported to the outer comet-forming zone ( Ábrahám et al. 2009, 2019). The direct effects of the outburst on the outer disk, however, are still yet to be fully characterized (e.g., Sicilia-Aguilar et al. 2012). The outer disk was only recently spatially resolved by Hales et al. (2018) using the Atacama Large Millimeter Array (ALMA) at 232 GHz. They found no indications of substructure in the disk and used radiative transfer modeling to constrain the total dust mass to 1.0 × 10 -4 M . Kóspál, et al. (2020) used spectropolarimetric monitoring of EX Lupi post-outburst with the Canada-France Hawaii Telescope (CFHT) to place constraints on the stellar magnetic fields. Their preliminary analysis suggests the surface magnetic field is 3 kG, stronger than what has been observed for some embedded FUor/EXor objects (e.g., 1 kG observed in FU Orionis Donati et al. 2005) and more in line with what is seen in Classical T Tauri stars with field strengths of up to a few kG (Johns-Krull 2007;[START_REF] Bouvier | Protostars & Planets[END_REF]Johns-Krull et al. 2009). The strength and stability of the stellar magnetic field in EX Lupi can influence the stability of the accretion columns (Sicilia-Aguilar et al. 2015) and has significant implications on possible outburst mechanisms.

In this paper, we present long wavelength observations of EX Lupi. In Sec. 2, we describe the data and calibration procedures. In Sec. 3, we outline the model fitting approach used and in Sec. 4 we discuss the results.

OBSERVATIONS

The analysis presented here uses a combination of new radio observations with the NSF's Karl G. Jansky Very Large Array (VLA), new ALMA 100 GHz observations, and ALMA 232 GHz data from literature (Hales et al. 2018). All of the observations are summarized in Table 1 and the new ones are detailed below.

VLA Observations

The VLA observations PI White) were centered on EX Lupi using J2000 coordinates RA = 16 h 03 m 7.09 s and δ = -40 • 18 05. 10. The Scheduling Blocks (SB) for each observation were executed on different days but all used the B configuration with 26 antennas and baselines ranging from 0.21 to 11.1 km. Quasar J1607-3331 was used for bandpass and gain calibration. Quasar 3C286 was used as a flux calibration source. Data were reduced using the CASA 5.4.1 pipeline, which included bandpass, flux, and phase calibrations [START_REF] Mcmullin | Astronomical Data Analysis Software and Systems XVI[END_REF]). The absolute flux calibration of the VLA at these wavelengths is typically ∼ 5%.1 All of the SBs used the same sources for calibration.

Table 1. Summary of the observations. The ALMA 232 GHz data is from Hales et al. (2018) with no additional calibration procedures. The peak fluxes were measured from the CLEANed images and the total flux was calculated in the CASA task uvmodelfit. The stated uncertainties do not include the absolute flux calibration uncertainty, which is ≤ 10% at these frequencies for ALMA and ≤ 5% for VLA. * The 45 GHz VLA observations resulted in a non-detection, therefore we include a 3σ upper level limit on the peak flux. The 15 GHz data from the VLA were acquired in Semester 19A on 2019 March 09 and had a total onsource time of 948 s. The observations used a Ku Band tuning setup with 3 × 2 GHz basebands and rest frequency centers of 13 GHz, 15 GHz, and 17 GHz. This gives an effective frequency of 15 GHz (1.99 cm). The data were imaged with a natural weighting and cleaned using CASA's CLEAN algorithm down to a threshold of 1 2 the RMS noise. The 15 GHz observations achieve a sensitivity of 7 µJy beam -1 . The size of the resulting synthesized beam is 1. 88 × 0. 42 (∼ 180 au) at a position angle of 10.5 • . Due to the large beam, the emission is only marginally resolved along the minor axis of the beam. The peak flux is 0.044 mJy beam -1 as measured in the CLEANed image. We used the CASA task uvmodelfit and a disk model to obtain a total flux of 0.050 ± 0.008 mJy.

The 45 GHz VLA data were acquired on 2019 May 13 and 2019 July 10 and had a total combined on-source time of 2345 s. Both 45 GHz SBs used a Q Band tuning setup with 4 × 2.048 GHz basebands and rest frequency centers of 41 GHz, 43 GHz, 47 GHz, and 49 GHz. This gives an effective frequency of 45 GHz (6.7 mm) for the Q band. The 45 GHz observations were concatenated in CASA and the data were imaged with a natural weighting and cleaned using the CLEAN algorithm down to a threshold of 1 2 the RMS noise. Together, these data achieve a sensitivity of 43 µJy beam -1 . The resulting synthesized beam is 0. 36 × 0. 16 (∼ 40 au) at a position angle of 10.2 • .

The low declination of EX Lupi leads to an elongated synthesized beam with the VLA (EX Lupi is located at δ = -40 • and the VLA is at a latitude of 34 • N). The atmospheric fluctuations at low altitude can also be much greater, leading to a poorer phase calibration and thus impacting the quality of the reconstructed images. These factors, coupled with less time on source than initially requested for the Q Band, led to a non-detection at 45 GHz. We note that this does not necessarily mean the disk is not observable at these frequencies.

ALMA Observations

EX Lupi was observed with ALMA (2017.1.00224.S, PI Kóspál) on 2018 January 27, March 17, and March 19 with a phase center of RA= 16 h 03 m 5.5 s δ = -40 • 18 25. 4. The total on-source integration time was 11640 s. The nominal antenna configuration was C43-4, with baselines between 14 m and 1398 m. Two spectral setups were used with spectral windows for different molecular lines (which will be presented in a later paper). In both setups, a 1.875 GHz wide window was centered at 100.2 GHz (2.99 mm) to measure the continuum emission of EX Lupi. The data were manually calibrated using the CASA 5.4.1 pipeline. The procedure included offline water vapor radiometer calibration, system temperature correction, and bandpass, phase and amplitude calibrations. Quasars J1427-4206 and J1517-2422 were used for pointing, bandpass, and flux calibration, and Quasar J1610-3958 was used for phase calibration. The sampled visibilities were Fourier transformed, creating dirty images of the source, which yields the skymodel of the target convolved with the point source function (PSF) of the beam. Dirty images for each spectral window were used to determine the frequencies without line emission and as input for the uvcontsub routine in CASA to obtain the continuum emission.

The 100 GHz observations were concatenated in CASA and the data were imaged with a natural weighting and CLEANed using the CLEAN algorithm down to a threshold of 1 2 the RMS noise. Together, these data achieve a sensitivity of 13 µJy beam -1 . The resulting synthesized beam is 1. 09 × 0. 81 (∼ 150 au) at a position angle of 87.8 • . The disk is marginally resolved. The peak flux is 1.98 Jy beam -1 as measured in the CLEANed image. We used the CASA task uvmodelfit and a disk model to obtain a total flux of 2.72 ± 0.013 µJy.

Literature Data

In addition to the new ALMA and VLA data presented here, we also use the ALMA 232 GHz continuum data from Hales et al. (2018) in our analysis. The observations were made on 2016 July 25, have a CLEANed synthesized beam of 0. 32 × 0. 26, and sensitivity of σ rms = 38 µJy beam -1 . The disk was well resolved at this frequency. No further calibration or processing was performed outside the procedure listed in Hales et al. (2018).

MODEL FITTING

To constrain the parameters of the dust in EX Lupi's circumstellar disk, we followed a radiative transfer (RT) model fitting approach similar to Hales et al. (2018). We use the RT code RADMC-3D 0.41 (Dullemond et al. 2012) with the Python interface radmc3dPy 2 to set the code input parameters for a given disk model. We keep the following parameters fixed in the fitting procedure: inclination i = 32 • , position angle PA= 65 • , and inner disk radius r = 0.05 au (Hales et al. 2018); flaring parameter ψ = 0.09 (Sipos et al. 2009); stellar temperature T = 3859 K and stellar radius R = 1.6 R (Frasca et al. 2017).

We adopt a disk model similar to that of a typical T Tauri protoplanetary disk (Andrews et al. 2009) with a disk density given by:

ρ = Σ(r, φ) H p (2π) exp - z 2 2H 2 p , (1)
where Σ is the surface density profile, H p is the pressure scale height, and z is the height above the disk midplane. The disk's surface density profile follow a power-law profile with an exponential outer tapering:

Σ(r) = Σ 0 r R c -γ exp - r R c 2-γ , ( 2
)
where Σ 0 is the surface density at the inner radius of 0.05 au, R c is the characteristic radius of the disk, and γ is the power-law exponent of the radial surface density profile. The pressure scale height is defined as:

H p = h c r 100 au 1+ψ , (3)
2 http://www.ast.cam.ac.uk/ juhasz/radmc3dPyDoc/index.html where ψ is the disk flaring parameter and h c is the ratio of the pressure scale height over radius at 100 au (see Hales et al. 2018).

In order to converge on the best fit disk parameters, we use a Metropolis-Hastings Markov Chain Monte Carlo (MCMC) model fitting approach. The free parameters considered in the modeling are: the total disk mass with a gas-to-dust ratio of 100:1, M disk , characteristic radius, R c , power law exponent of the surface density profile, γ, and scale height ratio at 100 au, h c . We perform the MCMC modelling in the image plane (see, e.g., Booth et al. 2016;White et al. 2017). After a trial model is selected, we use RADMC-3D and the mctherm command to calculate the dust temperature and then use the image command to generate a ray-traced continuum image projected to the inclination and position angle of the disk. The image is then attenuated by the primary beam and convolved with the dirty beam for a given observational setup. A χ 2 for each trial model is calculated as

χ 2 = (Data -M odel) 2 σ 2 , ( 4
)
where σ is the observed σ rms for a given observation multiplied by the synthetic beam size in pixels (see Booth et al. 2016). To fit multiple wavelength's data simultaneously, the χ 2 at each wavelength needs to be calculated and weighted. We adopt an equal weighting for each wavelength, and all of the χ 2 values are averaged together. A given trial model is then accepted if a random number drawn from a uniform distribution [0,1] is less than α, where

α = min(e 1 2 (χ 2 i -χ 2 i+1 ) , 1). ( 5
)
Fitting the 232 GHz, 100 GHz, and 15 GHz data simultaneously requires a dust opacity file extended to larger grain sizes. We use the OpacityTool 3 program (Toon & Ackerman 1981;Woitke et al. 2016) to get more realistic dust absorption and scattering parameters. This program calculates dust opacities by using a volume mixture of 60% amorphous silicates (e.g., Dorschner et al. 1995), 15% amorphous carbon (e.g. Zubko et al. 1996), and a 25% porosity. Bruggeman mixing is used to calculate an effective refractory index and a distribution of hollow spheres with a maximum hollow ratio of 0.8 (Min et al. 2005) is included to avoid Mie theory scattering artifacts. We further assume the disk is populated by 0.1 -21000 µm grains following a power-law size distribution of s -3.5 .

3

The OpacityTool Software was obtained from https://dianaproject.wp.st-andrews.ac.uk/data-resultsdownloads/fortran-package/

We were unable to adequately reproduce any wavelength's data with this approach. We tried varying the weight of each wavelength's χ 2 and could only begin obtaining reasonable results if the weighting for the 15 GHz data was set to an arbitrarily small value. This indicates that the approach is not well suited for fitting all three datasets simultaneously, or that there are different emission mechanisms at longer wavelengths. Therefore, we decided to exclude the 15 GHz data and discuss other approaches to fitting it in Sec.4.1.

To constrain the properties of the 232 GHz and 100 GHz observations, we use all the same approach and assumptions outlined above but change the particle size population to 0.1 -10000 µm grains. We adopted an equal weighting for the two data sets and ran an MCMC fit with 100×1000 link chains minus 100 each for burn-in. This approach of only fitting the ALMA data was able to well reproduce the observations. The most probable values of the free parameters and 95% credible region (CR) are summarized in Table 2. The resulting best fit model and residuals are shown in Fig. 1 and the posterior distribution functions (PDF) are shown in Fig. 2.

DISCUSSION

Model Fitting Results

Our RT models of the 232 and 100 GHz ALMA data were able to well reproduce the observations. When fitting both of the datasets simultaneously, we find most probable values of: M disk = 0.01 M , R c = 45 au, γ = 0.25, and h c = 2.5 au (the corresponding 95% confidence intervals are listed in Table 2). The total flux of the most probable models at 232 and 100 GHz are 18.2 and 2.45 mJy, respectively. Most of the results are consistent with the values reported in Hales et al. (2018), within the uncertainties, which only fit the ALMA 232 GHz data. We find though that the most probable value of the characteristic radius is about a factor of 2 larger. This difference in radius could be due to the larger beam size in the 100 GHz observations, along with equal weighting, which is forcing the models to be larger. The spectral index between 232 and 100 GHz is α 1.3-3.0mm = 2.20 ± 0.11 (the uncertainty includes a 10% absolute flux calibration uncertainty at each frequency). The spectral index is consistent with 2.19 ± 0.47 as reported in Hales et al. (2018) and calculated within the ALMA Band 6 spectral windows. Extrapolating the peak flux of the 100 GHz data to 45 GHz, along with the new beam size and α = 2.2, gives a peak flux lower than the achieved σ rms of the 45 GHz observations. This calculated flux shows that the nondetection at 45 GHz is not useful for constraining the thermal disk properties.

If the 15 GHz observations are tracing the thermal emission from large grains, then it is possible these large grains are located at a different area of the disk than the smaller grains probed by ALMA. To test this, we tried fitting the VLA 15 GHz data alone with the same RT approach as for the ALMA datasets. This approach allows for a different disk geometry of the large grains, but still requires the total flux to be well fit to the data. We note though that due to the large beam size at 15 GHz, the disk is not resolved along the major axis and would be only marginally resolved along the minor axis. The results are summarized in Table 2, the best fit model and residuals are shown in Fig. 3, and the PDF is shown in Fig. 4.

While the 15 GHz only model does indeed well reproduce the data, the resulting disk parameters seem improbable. As segregation by grain size is possible due to radial drift and settling, a different disk geometry alone is not immediately disqualifying. The scale height, however, seems un-physically large at 89 au. Larger grains would be expected to settle in the midplane of the disk, meaning the scale height would likely be the same or smaller than is observed with the mm grains. The best fit total disk mass is also 3.5× larger than when fitting to the mm data alone. Such a disk may have emission at 45 GHz, depending on the spectral index, and we report a non-detection at 45 GHz. We conclude that even though a thermal emission model can technically fit the data well at 15 GHz, the resulting disk parameters to do so are highly improbable. Alternate sources of 15 GHz emission are explored in Sec. 4.3.

Millimeter Observations & Grain Growth

EX Lupi lies between the Lupus 3 and Lupus 4 star forming regions (Cambrésy 1999). ALMA surveys of protoplanetary disks in the Lupus star forming regions have found total disk masses to be ∼ 10 -3 M (Ansdell et al. 2016). Our model fitting finds a total disk mass nearly an order of magnitude higher, assuming a gasto-dust ratio of 100:1. The mass discrepancy could be due to the assumptions made in the mass calculation in Ansdell et al. (2016) such as the optical depth and that here we include a full RT calculation. The actual gasto-dust ratio in EX Lupi could also be much lower than assumed. If the difference in the masses is indeed real, it could be explained by EX Lupi being younger, more heavily accreting, or from an inherent difference in the disks of EXor/FUor-type system from that of typical protoplanetary disks. EX Lupi's disk mass is smaller than that of FUors, as expected, but the characteristic radius is similar. Other ALMA studies have found that EXor/FUor disks tend to be more compact than that Table 2. Summary of the most probable model parameters and the corresponding 95% confidence intervals from the posterior distributions in brackets. The first row is the ALMA 232 GHz and 100 GHz data fit simultaneously, the second row is the 15 GHz VLA data fit alone, the third row is fitting a point source model on top of the RT calculated 15 GHz disk emission. The disk mass is the total disk mass assuming a gas-to-dust ratio of 100:1. The scale height is measured at 100 au (see Hales et al. 2018) of typical protoplanetary disks (e.g., Cieza et al. 2018, Kóspál et al. in prep.). EX Lupi experienced an outburst in 2008 and has since returned to a quiescent state. Therefore, all of the ALMA and VLA data presented here are indicative of a post-outburst circumstellar environment. Ábrahám et al. (2009) found that the outburst increased the crystallinity of the disk grains and transported them from the inner regions to the outer regions. Their work shows that even a short-lived outburst, such as an EXor-type outburst, can have observable effects on the circumstellar disk.

The water ice line (or snow line) is the radial location in a disk where water reaches its condensation temperature and freezes out on to grains in the disk. The exact temperature at which this occurs depends on other disk properties, such as gas pressure, but is typically between 150 to 175 K in the midplane of a disk (Lecar et al. 2006). EX Lupi's recent outburst increased the bolometric luminosity by about 36× ( Ábrahám et al. 2009). This enhanced luminosity can heat the disk and push the ice line(s) out to further radial distances. This was observed in ALMA observations of V883 Ori, where the ice line moved from a presumed 1 -5 au pre-outburst to a current location of 42 au as measured on the surface of the disk (Cieza et al. 2016). As an outburst fades and the temperature drops, the ice line will move back inwards to the pre-outburst location. From Ábrahám et al. (2009, Fig. 4 in supplementary materials), we estimate that during EX Lupi's outburst the ice line was located at 20 -30 au on the surface of the disk and 1 -2 au in the midplane of the disk. Using Spitzer observations of the H 2 O spectra, Banzatti et al. (2012) find the ice line moved from 1.3 au during outburst to 0.6 au after outburst.

Our RT modeling shows that the ice line is now located at 3 -4 au on the surface of the disk and 0.2 -0.3 au in the disk's midplane, indicating a significant shift inwards post-outburst. To check if the disk could have reasonably cooled off enough to make this change in the position of the ice line, we explore the disk relaxation timescales (see Lin & Youdin (2015) and Section 2.3 of Flock et al. 2017). The density and temperature we take from our best-fit RT model, the gas-to-dust ratio is set to 100:1, and we choose a representative cooling wavelength opacity at 18 µm of κ = 52.7 cm 2 /g. The relaxation time curve is shown in Fig. 5.

The black line represents the relaxation timescale in years for the best fit disk model as a function of radial position in the disk, with the outburst and post-outburst ice lines indicated in gray. This curve shows that EX Lupi's recent outburst, of duration of ∼ 1 year, could only significantly heat the midplane of the disk out to radii of 1-2 au, consistent with data from during the outburst ( Ábrahám et al. 2009;Banzatti et al. 2012). While the current data cannot be used to infer the location of the ice line prior to the outburst, the outburst and postoutburst data together show evidence that the location of the ice line can indeed shift during EXor-type events. The plot also shows that the ∼ 10 yr time between the end of the outburst and the observations is more than enough to relax the disk up to 10 au to a new thermal state.

An ice line can be a source of turbulence in a disk, around which grain growth can occur (e.g., [START_REF] Bouvier | Protostars & Planets[END_REF]Ros & Johansen 2013;Zhang et al. 2015). There- fore as the ice line moves inwards after an outburst it provides a potential mechanism to spur grain growth throughout the disk. The associated timescales for water to deposit onto grains near the ice line will be much shorter than the relaxation timescale outlined above, meaning grain growth from water freeze-out or deposition can be commensurate with the moving position of the ice line (e.g., Brown & Charnley 1990). Thermal grains are inefficient emitters at wavelengths longer than their grain size. Therefore, using the detection of the disk in the ALMA 100 GHz data, and calculated spectral index of α = 2.20 ± 0.01 indicating the disk is optically thin, we infer that grains of at least 3 mm may be present in the disk post-outburst. Our RT models include grains of sizes up to 1 cm. While this is not a confirmation of 1 cm grains it does indicate the data is at least consistent with the presence of up to 1 cm grains. Significant grain growth is possible with EX Lupi's evolving ice line, although we note that there is no pre-outburst ALMA and VLA data to compare it to. If EXor-like outbursts, which are short lived and can repeat many times during pre-main-sequence stellar evolution (e.g., Vorobyov & Basu 2015;[START_REF] White | Science with a Next Generation Very Large Array[END_REF], are indeed common for most stars then appreciable grain growth can occur throughout the disk both early and often. Quick and abundant grain growth to mm-cm sized particles is an important step in the planet formation process and can drive the growth of larger planets and planetesimals (e.g., Morbidelli & Raymond 2016;Johansen & Lambrechts 2017;Hughes & Boley 2017). Significant early grain growth and the presence of planetesimals is also a possible explanation for the gaps seen in HL Tau (Brogan et al. 2015) or surveys of nearby protoplanetary disks with the ALMA DSHARP survey (e.g., Andrews et al. 2018).

Long Wavelength Central Emission

In Sec. 3, we concluded that the 15 GHz VLA data is likely dominated by something other than thermal disk emission. A more likely scenario is that at 15 GHz we are seeing thermal disk emission plus a combination of nonthermal disk emission, stellar winds, jets, or stellar emission. To test this scenario, we can fit a point source emission model on top of the thermal disk emission. Given the resolution of the 15 GHz data (1. 88 × 0. 42), we are unable to differentiate emission from a point source and a moderately extended (i.e., a few au) region. Taking the best fit RT model from the ALMA data alone, we use RADMC-3D to calculate the emission at 15 GHz, resulting in a very faint disk with a total flux of only ∼ 0.001 mJy. Starting with this disk model, the "additional" emission required to reproduce the data can be assumed to be coming from approximately the center of the disk. This leaves only 3 free parameters: the flux of the central emission, the X-offset, and the Y-offset. We adopt an MCMC modelling approach similar to the one outlined in Sec. 3 but now just add the flux of the central emission on top of the already calculated RT model. The results are summarized in Table 2, the best fit model and residuals are shown in Fig. 6, and the PDF is shown in Fig. 7. We find a most probable flux of 0.055 mJy, an X-offset of 0. 025, and a Y-offset of 0. 26. Due to the nature of the large beam at 15 GHz (1. 88×0. 42), the location of the central emission flux is still broadly consistent with being peaked on the star itself.

Role of the Magnetic Field

Some of the most promising outburst theories for EXors rely on the role of the magnetic field. Armitage (2016) predicted that the outbursts could be explained by changes in the polarity and strength of the stellar magnetic fields at the kG-level. In a competing scenario, D'Angelo & Spruit (2012) proposed an instability which can lead to quasi-periodic oscillations in the inner disk and associated recurrent outbursts. This instability can occur when the accretion disk is truncated close to the co-rotation radius by the strong magnetic field of the star. Stable accretion columns linked to a very strong magnetic field in EX Lupi were noted in Sicilia-Aguilar et al. (2015).

In principle, if kG-level stellar magnetic fields were present in EX Lupi then there could be significant nonthermal emission that would yield brightness temperatures easily detectable at long wavelengths with VLA. The central emission flux of 0.055 mJy corresponds to a Rayleigh-Jeans brightness temperature of T B ≈ 4 × 10 7 K, assuming the size of the emitting region is uniformly spread out over the surface of the star (we note that depending on the emission mechanism the actual emitting region could range from a small localized area of the star to several stellar radii). This brightness temperature is much larger than the effective temperature of EX Lupi (T = 3859 K), which indicates that significant non-thermal stellar emission could be present as was seen in optical line emission (e.g., Sicilia-Aguilar et al. 2012).

One potential source of such large brightness temperatures is synchrotron emission from relativistic or nearly relativistic electrons being accelerated by EX Lupi's magnetic field. If synchrotron emission is present, then the flux should peak near the critical frequency (e.g., Hughes et al. 2019),

ν crit = γ 2 qB 2πm e , ( 6
)
where γ is the Lorentz factor and is assumed to be ∼ 1, q is the electron charge, B is the magnetic field strength, and m e is the electron mass (we note that γ here is not the same as the disk power-law exponent used in Eqn. 2). Adopting a magnetic field strength of 3 kG gives ν crit = 8.4 GHz. Since the flux should quickly drop off at frequencies larger than the critical frequency, it would be unlikely to observe such a large flux at 15 GHz. However, if the electrons are slightly more relativistic, with γ = 1.34, then ν crit = 15 GHz. Alternatively, Armitage (2016) predict that significant changes in magnetic field strength could be a driver for the episodic accretion in EXor-type stars. Therefore, the magnetic field strength could potentially be different than measured previously with CFHT (Kóspál, et al. 2020). If this is indeed the case, and ν crit is ∼ 15 GHz, then a magnetic field strength of ∼ 5.4 kG can be inferred.

The effective temperature (see, e.g., Pacholczyk 1970) of electrons emitting synchrotron radiation at a fre- quency ν, is given by:

T e ≈ 2πm e ν qB 1/2 m e c 2 3k ( 7
)
where k is the Boltzmann constant. This yields 2 -3 × 10 9 K, depending on which value for the magnetic field is used. The electron temperature is about a factor of 50 larger than the observed brightness temperature implying the emission is optically thin. At face value, synchrotron emission is a potential source of the 15 GHz emission. An important consideration though is that the size of the emitting region could be extended much further from the stellar surface. At further stellar separations, the magnetic field strength will become weaker causing ν crit to fall below 15 GHz and the expected synchrotron emission at 15 GHz to be much lower as well.

While EX Lupi is indeed still a pre-main-sequence star, its spectral classification is M0 and its mass indicates that when it reaches the main-sequence it will still be an M-type star. M-type stars are notorious for flares and magnetic activity, but were not observed at radio wavelengths until recent technological advancements (Berger et al. 2001). While the exact emission mechanisms are still debated, they are thought to be primarily due to electron cyclotron maser instabilities (ECMI) and/or gyrosynchrotron emission.

ECMI emission is similar to the auroral emission observed on most solar system planets (Turnpenney et al. 2018). If the 15 GHz emission in EX Lupi was due to ECMI, then the emission is expected to peak at the fundamental critical frequency of ν crit = 8.4 GHz and fall off rapidly at higher frequencies. ECMI requires a relatively stable magnetic field configuration, such as in planetary or brown dwarf magnetic fields. Therefore, even though EX Lupi's magnetic field strength could have changed between the CFHT and VLA observations such that ν crit ∼ 15 GHz, ECMI is a highly unlikely source of the 15 GHz emission due to the required stability.

If the 15 GHz emission in EX Lupi was due to gyrosynctrotron emission, then it is likely due to magnetic reconnection events releasing a large number of energetic particles (e.g., Williams et al. 2014;Hughes et al. 2019). The surface of a pre-main-sequence accreting Mtype star, such as EX Lupi, is undoubtedly a turbulent environment where these processes could dominate. The expected emission at 15 GHz depends on the size of the emitting region (which is typically much smaller than the stellar radius), the magnetic field strength, and the electron energy index δ. Assuming the gyrosynchrotron emission is optically thin, then the spectral index can be used as a tracer for δ. Following Hughes et al. (2019), the relation between all of these parameters is given by:

R em = 132 d pc GHz ν F µJy T ef f , (8)
where R em is the emitting region in units of R , F µJy is the observed flux in µJy, and T eff is given by: T eff = 2.2×10 9 (sin θ) -0.36-0.06δ 10 -0.31δ ν ν crit 0.50+0.085δ

(9) where θ is assumed to be ∼ 90 • (see Dulk 1985, for the derivation and further details on the relationship between expected flux and the size of the emitting region). Fig. 8 shows the size of the emitting region and magnetic field strength for various values of δ (which cannot be constrained with the available data). If we assume a typical lower level limit of δ ∼ 2 and a magnetic field strength range of 3 -6 kG then the corresponding size of the emitting region is ∼ 0.3 -0.5 R . This is a significant fraction of the stellar surface (R * = 1.6 R ) and much larger than is typically expected for mainsequence M-dwarf stars which have reconnection regions typically of order 0.01 R . If δ > 2, which is typically the case in M-dwarfs, then the size of the emitting region quickly becomes larger than the star. EX Lupi could, however, have a very small value of δ or have the gyrosynchrotron emission come from a large area where it is actively accreting disk material. Follow-up observations at higher/lower frequencies will enable spectral index constraints, and thus a value for δ, to better determine if gyrosynchrotron is a potential source of the observed 15 GHz emission.

Radio Jets and Non-thermal Disk Emission

Some EXor/FUor systems also have very bright radio jets (e.g., Z CMa and L1551 IRS 5, Poetzel et al. 1989;Rodríguez et al. 2003). While EX Lupi did show significant X-ray activity leading up to and during the most recent outburst, the emission is most likely stemming from accretion shocks instead of jets (see, e.g., Grosso et al. 2010). No radio jets have been previously reported in EX Lupi. If the 15 GHz emission was indeed due to a radio jet, then we would expect it to have an approximately flat spectral index leading it to possibly be detectable at 45 GHz. Since we did not detect anything at 45 GHz, and no radio jet was reported previously, we note the likelihood of the 15 GHz emission being jetdriven is low.

Centrally located disk winds are another possible source of long wavelength emission. Kóspál et al. (2011) found that the hydrogen emission within 1 au is likely due to disk winds in EX Lupi. In optical line spectra, there is further evidence of disk winds in both outburst and quiescence (Sicilia-Aguilar et al. 2012, 2015;Banzatti et al. 2019). Hales et al. (2018) observed large scale asymmetries in the outer CO gas disk which could be due to a molecular outflow. Similar to jets, disk winds would likely be peaked at frequencies < 15 GHz and have a flat to slightly positive spectral index. Followup observations at < 15 GHz are necessary in order to measure the spectral index and determine if winds or jets are a possible source of observed emission.

Free-free emission from an ionized disk, typically more prevalent in the latter stages of disk dissipation when photoevaporation becomes significant, can come from EUV and X-ray irradiation. EX Lupi had significant X-ray activity before and during the outburst, which could have significantly ionized its inner disk. Although unlikely, if this irradiation continued after the outburst ended then free-free emission could still be present in the EX Lupi system. Free-free emission can also arise from accretion shocks propagating through the disk when the accretion rate becomes much larger (e.g., Hartmann & Kenyon 1996). Sicilia-Aguilar et al. (2015) find evidence of accretion shocks via optical spectroscopy and Grosso et al. (2010) noted that the X-ray emission seen in EX Lupi is also likely due to accretion shocks. The X-ray luminosity of EX Lupi was 1.7 × 10 30 erg s -1 during the outburst (Grosso et al. 2010). Using the relation for X-ray luminosity to expected radio flux from free-free emission outlined in [START_REF] Pascucci | [END_REF]:

F 3.5 cm = 2.4 × 10 -29 51 d 2 L x [µJy], (10)
where F 3.5 cm is the 3.5 cm flux, d is the distance to the source in pc and L X is the X-ray luminosity in erg s -1 , we find an expected flux of 4.3 µJy. Free-free emission should have an approximately flat spectral index meaning the expected flux at 15 GHz (2 cm) should also be ∼ 4.3 µJy. It is therefore possible that free-free emission accounts for up to 10% of the observed 15 GHz flux. We note though that since the post-outburst X-ray luminosity should be lower than that observed by Grosso et al. (2010), the expected post-outburst radio flux should be lower as well. Considering all of the potential sources of the 15 GHz emission, we find that the most likely scenario is a combination of (gyro)synchrotron and free-free emission.

We note though that follow-up observations at lower frequencies are needed to confirm the spectral index. It is indeed possible that there is a combination of more emission mechanisms present at 15 GHz. In order to disentangle all the potential sources of emission, 1-15 GHz observations with both high angular resolution and sensitivity are necessary. The VLA is currently the leading facility in this regime and the observations presented are already pushing the limits of its capabilities. Therefore, proposed future facilities, such as the ngVLA [START_REF] White | Science with a Next Generation Very Large Array[END_REF], will be key to fully understanding the underlying emission mechanisms in EX Lupi and connecting them in the broader context of EXor/FUors in general.

SUMMARY

In this paper, we presented ALMA and VLA continuum observations of the EX Lupi disk in its postoutburst state. We fit radiative transfer models of the circumstellar dust and find the models are consistent with grain growth up to at least 3 mm, and possibly as high as 1 cm. The grain growth could have been spurred by the recent outburst, which ended in 2008. The best fit disk model has an ice line located at 0.2 -0.3 au in the disk's midplane, and the associated cooling timescales show there as been adequate time for the ice line to migrate inward from the observed position during the outburst. The most probable value for the total disk mass is 0.01 M , accompanied by a relatively compact characteristic dust radius of 45 au. The size is in agreement with other studies that find EXor/FUor disks to be more compact than that of typical protoplanetary disks.

We find a best fit flux of 0.055 mJy at 15 GHz, significantly more than can be explained by thermal disk emission alone. We explored several potential sources of the emission and conclude that it is likely primarily due to (gyro)synchrotron emission coming from strong stellar magnetic fields and/or non-negligible free-free emission from accretion shocks and disk heating through X-ray emission. ence Foundation operated under cooperative agreement by Associated Universities, Inc.

Facilities: ALMA, VLA Software: CASA 5.4.1 [START_REF] Mcmullin | Astronomical Data Analysis Software and Systems XVI[END_REF]

Figure 1 .

 1 Figure 1. Top: Dirty image, best fit RT model, and residuals for the ALMA 232 GHz observations. Bottom: Dirty image, best fit RT model, and residuals for the ALMA 100 GHz observations. The two datasets were fit simultaneously with equal weighting. The black ellipses in the bottom left of each figure represent the synthesized beam.

Figure 2 .

 2 Figure 2. Posterior distribution of the best fit model for fitting the ALMA 232 GHz and 100 GHz simultaneously. The most probable values for each parameter are denoted by the blue lines.

Figure 3 .

 3 Figure 3. Left: VLA 15 GHz dirty image. Middle: 15 GHz RT disk model Right: Data minus model residuals. The black ellipses in the bottom left of each figure represent the synthesized beam.

Figure 4 .

 4 Figure 4. Posterior distribution of the best fit model of the 15 GHz emission alone shown in Fig. 3. The most probable values are denoted by the blue lines.

Figure 5 .

 5 Figure 5. The relaxation time in years as a function of disk radii along the midplane for the best fit RT model of EX Lupi (see Eqns. 10-13 in Flock et al. 2017). The locations of the ice line along the disk midplane both during ( Ábrahám et al. 2009) and post-outburst are also denoted in gray.

Figure 6 .

 6 Figure 6. Left: VLA 15 GHz dirty image. Middle: 15 GHz RT disk model (as calculated from fitting to both ALMA datasets) added to the best fit central emission model. Right: Data minus model residuals. The black ellipses in the bottom right represents the synthesized beam.

Figure 7 .

 7 Figure 7. Posterior distribution of the 15 GHz central emission on top of the RT calculated disk emission from the ALMA data. The most probable values are denoted by the blue lines.

Figure 8 .

 8 Figure 8. The size of the emitting region and magnetic field strength for various electron energy indices, δ, assuming gyrosynchrotron emission at 15 GHz. The solid black line at δ = 2 represents a typical lower level estimate for δ.

  .

	Dataset(s)	Disk Mass	Characteristic Radius	γ	Scale Height
		(M )	(au)		(au)
	232 GHz + 100 GHz	0.0099	45	0.25	2.5
		[0.0067, 0.012]	[39, 68]	[-0.36, 0.77]	[1.9, 3.8]
	15 GHz	0.035	56	0.12	89
		[0.011, 0.053]	[13, 140]	[-1.8, 1.2]	[21, 98]
	Dataset	Flux	X-offset	Y-offset
		(mJy)	( )	( )
	15 GHz	0.055	0.025	0.26
		[0.045, 0.066]	[-0.09, 0.11]	[-0.010, 0.55]

For a note on the VLA flux calibration uncertainty, see science.nrao.edu/facilities/vla/docs/manuals/oss/performance/fdscale.

We thank the anonymous referee for feedback that improved this paper. This project received support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No 716155 (SAC-CRED). DS acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1833: "Building a Habitable Earth" (SE 1962/6-1). EV and VA ac-knowledge the support of the Large Scientific Project of the Russian Ministry of Science and Higher Education "Theoretical and experimental studies of the formation and evolution of extrasolar planetary systems and characteristics of exoplanets" (No. 13.1902.21.0039). On behalf of the SACCRED project we are thankful for the usage of MTA Cloud (https://cloud.mta.hu/) that helped us achieving the results published in this paper. This paper makes use of the ALMA data from projects ADS/JAO.ALMA#2017.1.00224.S and ADS/JAO.ALMA#2015.1.00200.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Sci-