Information-Theoretic Testing and Debugging of Fairness Defects in Deep Neural Networks

Verya Monjezi

University of Texas at El Paso

Gang Tan

Pennsylvania State University
Ashutosh Trivedi

University of Colorado Boulder

Saeid Tizpaz-Niari

University of Texas at El Paso

The deep feedforward neural networks (DNNs) are increasingly deployed in socioeconomic critical decision support software systems.
DNNs are exceptionally good at finding minimal, sufficient statistical patterns within their training data.
Consequently, DNNs may learn to encode decisions—amplifying existing biases or introducing new ones—that may disadvantage protected individuals/groups and may stand to violate legal protections.
While the existing search based software testing approaches have been effective in discovering fairness defects, they do not supplement these defects with debugging aids—such as severity and causal explanations—crucial to help developers triage and decide on the next course of action.
Can we measure the severity of fairness defects in DNNs?
Are these defects symptomatic of improper training or they merely reflect biases present in the training data?
To answer such questions, we present : an information-theoretic testing and debugging framework to discover and localize fairness defects in DNNs.

The key goal of  is to assist software developers in triaging fairness defects by ordering them by their severity.
Towards this goal, we quantify fairness in terms of protected information (in bits) used in decision making.
A quantitative view of fairness defects not only helps in ordering these defects, our empirical evaluation shows that it improves the search efficiency due to resulting smoothness of the search space.
Guided by the quantitative fairness, we present a causal debugging framework to localize inadequately trained layers and neurons responsible for fairness defects.
Our experiments over ten DNNs, developed for socially critical tasks,
show that  efficiently characterizes the amounts of discrimination,
effectively generates discriminatory instances (vis-a-vis the
state-of-the-art techniques), and localizes layers/neurons with significant biases.

§ INTRODUCTION

AI-assisted software solutions—increasingly implemented as deep neural networks [1] (DNNs)—have made substantial inroads into critical software infrastructure where they routinely assist in socio-economic and legal-critical decision making [2].
Instances of such AI-assisted software include software deciding on recidivism, software predicting benefit eligibility, and software deciding whether to audit a given taxpayer.
The DNN-based software development, driven by the principle of information bottleneck [3], involves a delicate balancing act between over-fitting and detecting useful, parsimonious patterns.
It is, therefore, not a surprise that such solutions often encode and amplify pre-existing biases in the training data. What's worse, improper training may even introduce biases not present in the training data or irrelevant to the decision making.
The resulting fairness defects may not only disadvantage protected groups [4, 5, 6, 7, 8], but may stand to violate statutory requirements [9].

This paper presents ,
an information-theoretic testing and debugging framework for fairness defects in deep neural networks.

Quantifying Fairness. Concentrated efforts from the software engineering and the machine learning communities have produced a number of successful fairness testing frameworks [10, 11, 12, 13].
These frameworks characterize various notions of fairness—such as group fairness [14] (decision outcome for various protected groups must be similar) and individual fairness [15] (individuals differing only on protected attributes must receive similar outcome)—and employ search-based testing to discover fairness defects.
While a binary classification of fairness is helpful in discovering defects, developers may require further insights into the nature of these defects to decide on the potential “bug fix”.
Are some defects more severe than others? Whether these defects stem from biases present in the training data, or they are artifacts of an inadequate training? Is it possible to find an alternative explanation of the training data that does not use protected information?

Individual discrimination is a  well-studied [16, 11, 17, 18] causal notion of fairness that defines a function being discriminant towards an individual (input) if there exists another individual (potentially counterfactual), differing only in the protected features, receives a more favorable outcome.
We present a quantitative generalization of this notion as the quantitative individual discrimination (QID). We define QID as the amount of protected information—characterized by entropy metrics such as Shannon entropy and min entropy—used in deriving an outcome.
Observe that a zero value for the QID measure implies the absence of the individual discrimination.
The QID measure allows us to order various discriminating inputs in terms of their severity, as in an application that is not supposed to base its decisions on protected information, inputs with higher dependence indicate a more severe violation.
Our first research question (RQ1) concerns the usefulness of QID measure in finding inputs with different severity.

Search-Based Testing.
Search-based software testing  provide scalable optimization algorithms to automate discovery of software bugs.
In the context of fairness defects, the search of such bugs
involves finding twin inputs exhibiting discriminatory instances.
The state-of-the-art algorithms for fairness testing [19, 17, 18] explore the input space governed by a binarized feedback, resulting in a discontinuous search domain. On the other hand, QID-based search algorithms can benefit a smooth (quantitative) feedback during the optimization, resulting in a more guided search.
Our next research question (RQ2) is to investigate whether this theoretical promise materializes in practice in terms of discovering richer discriminating instances than using classic notions of discrimination.

Causal Explanations.
While the discriminating instances (ordered by their severity) provide a clear evidence of fairness defects in the DNN, it is unclear whether these defects are inherent in the training data, or whether they are artifacts of the training process.
Inspired by the notion of “the average causal effects” [20] and Audee framework [21] for bug localization in deep learning models, we develop a layer and neuron localization framework for fairness defects. If the cause of the defects is found to be at the input layer, it is indicative of discrimination existing in the training data.
On the other hand, if we localize the cause of the defect to some internal layer, we wish to further prod the DNN to extract quantitative information about neurons and their counterfactual parameters that can mitigate the defect while maintaining the accuracy.

This debugging activity informed our next research question (RQ3):
is it possible to identify a subset of neurons and their causal effects on QID to guide a mitigation without affecting accuracy?

Experiments.  implements a search algorithm (Algorithm <ref>) to discover inputs that maximize QID and a causal debugging algorithm (Algorithm <ref>) to localize layers and neurons that causally affect the amounts of QID.
Using $10$ socio-critical DNNs from the literature of algorithmic fairness, we show that  finds inputs
that can use significant amounts of protected information in the decision making; outperforms three
state-of-the-art techniques [19, 17, 18] in generating discriminatory instances;
and localizes neurons that guides a simple mitigation strategy to reduce QID down to $15\%$ of reported initial QID with at most $5\%$ loss of accuracy.
The key contributions of this paper are:

* Quantitative Individual Discrimination. We introduce an information-theoretic characterization of discrimination, dubbed quantitative individual discrimination (QID), based on Shannon entropy and Min entropy.

* Search-based Testing. We present a search-based algorithm to discover circumstances under which the DNNs exhibit severe discrimination.

* Causal Debugging. We develop a causal fairness debugging based on the language of interventions to localize the root cause of the fairness defects.

* Experimental Evaluation. Extensive experiments over different datasets and DNN models that show feasibility, usefulness, and scalability (viz-a-viz state-of-the-art). Our framework can handle multiple protected attributes and can easily be adapted for regression tasks.

§ PRELIMINARIES

Fairness Terminology.

We consider decision support systems as binary classifiers where a prediction label is favorable if it gives a desirable outcome to an input (individual). These favorable predictions may include higher income
estimations for loan, low risk of re-offending in parole assessments,
and high risk of failing a class.
Each dataset consists of a number of attributes (such as income,
experiences, prior arrests, sex, and race) and a set
of instances that describe the value of attributes for each individual.
According to ethical and legal requirements, data-driven software should not discriminate
on the basis of an individual's protected attributes such as sex, race, age, disability,
colour, creed, national origin, religion, genetic information, marital status, and sexual orientation.

There are several well-established fairness definitions. Group fairness requires that the statistics of ML outcomes for different protected groups to be similar [14] using metrics
such as equal opportunity difference (EOD), which is the difference between
the true positive rates (TPR) of two protected groups. Fairness through unawareness (FTU) [15] requires removing protected attributes during training.
However, FTU may provide inadequate support since protected attributes can influence the prediction via a non-protected collider attribute (e.g., race and ZIP code).
Fairness through awareness (FTA) [15] is
an individual fairness notion that requires that two individuals that deemed similar (based on their non-protected attributes) are treated similarly. Our approach is geared toward individual fairness.

Individual Discrimination.
Causal discrimination, first studied in Themis [10],
measures the difference between two subgroups via counterfactual queries.
It samples individuals with the protected attributes set to $A$ and compares the outcome to a counterfactual scenario where the protected attributes is set to $B$.
Individual discrimination (ID) is a prevalent notion that adapts counterfactual queries to find an individual such that their counterfactual with a different protected attributes receives more favorable outcome.
This fairness notion is used by
the state-of-the-art fairness testing to generate fairness defects [22, 19, 17, 18] and closely related to situation testing notion [23].
While standard group fairness metrics (e.g., AOD/EOD) are already quantitative, the quantitative measures do not exist for individual fairness. We propose to adapt information theoretic tools to provide quantitative measures for individual fairness.

Information-Theoretic Concepts.

The notion of Rényi entropy [24], $H_\alpha(X)$ quantifies the uncertainty (randomness) of a system responding to inputs $X$.
In particular, Shannon entropy ($\alpha{=}1$) and min-entropy ($\alpha{=}\infty$) are two important subclasses of Rényi entropy. Shannon entropy ($H_1$) measures the expected amounts of uncertainty over finitely many events whereas min entropy ($H_{\infty}$) measures the uncertainty over single maximum likelihood event.

Consider a deterministic system (like pre-trained DNN) with a finite set of responses and assume that the input $X$ is distributed uniformly. Thus, the system induces an equivalence relation over the input set $X$ such that two inputs are equivalent if their system outputs are approximately close, i.e. $x {\sim} x'$ iff $\dnn(x)\approx_{\epsilon}\dnn(x')$. Let $X_o$ denote the equivalence class of $X$ with output $o$. Then, the remaining uncertainty after observing the output of DNN over $X$ can be written as:

\[
H_1(X|O)=\sum_{O=o}\frac{|X_o|}{|X|}.\log_{2}(|X_o|)~~~~~~~~(\text{Shannon entropy})
\]

where $|X|$ is the cardinality of $X$ and $|X_o|$ is the size of equivalence class of output $o$. Similarly, the min-entropy is given as

\[
\]

where $|O|$ is the number of equivalence classes over $X$ [25, 26, 27]. Given that the initial entropy is equal to $\log_2(|X|)$ for both entropies, the amount of information from $X$ used by the system to make decisions are

\begin{eqnarray*}
I_1(X;O) &=& \log_2(|X|) - H_1(X|O), \text { and }\\ I_{\infty}(X;O) &=& \log_{2}(|O|),
\end{eqnarray*}

under Shannon- and min- entropies with $I_{1} \leq I_{\infty}$.

Quantitative Notion and Fairness.

Our approach differs from these state-of-the-art techniques [28, 22, 19, 17, 29, 18]
in that it extends the individual discrimination notion with quantitative information flow
that enables us to measure the amount of discrimination in the ML-based software systems.

Given non-protected attributes and ML outcomes,
the Shannon entropy measures the expected amount of individual discrimination over all possible responses varying protected classes, whereas, min entropy measures the amount over a single response from the maximum likelihood class.

Consider a dataset with $16$ different
protected values—sex(2), race(2), and age(4)—distributed uniformly,
and suppose that we have $4$ individuals in the system. We perturb the protected attributes of these individuals to
generate $16$ counterfactuals and run them through the DNN to get their prediction scores.
Suppose that the outputs have all $16$ outputs in the same class for the first individual (absolutely fair) and
have $16$ classes of size one for the second individual (absolutely discriminatory).
For the third individual, let the outputs be in $4$ classes with $\set{4,4,4,4}$
elements in each class (e.g., there is one output class per each age group).
For the fourth individual, let us consider outputs to be in $5$ classes with $\set{8,4,2,1,1}$
elements in each class (e.g.,  race=1  the output is class $1$; else,
sex=1  the output is $2$; else
age={1,2}  the output is $3$; else
there is one output class for each age={3,4}).

We work with the following notions of discrimination.

* Individual Discrimination Notion. The individual discrimination used by the
state-of-the-art techniques can only
distinguish between the first individual and the rest, but they cannot distinguish
among individuals two to four. In fact, these techniques generate
tens of thousands of individual discriminatory instances in a short amount of time [19, 17, 18].
However, they fail to prioritize test cases for mitigation and
cannot characterize the amounts of discrimination (i.e., their severity).

* Shannon Entropy. Using the Shannon entropy, we have the initial
fairness to be $4.0$ bits, a maximum possible discrimination. The remaining fairness
of DNN are $4.0$, $0.0$, $2.0$ and $2.125$, for the first to fourth individuals, respectively.
The discrimination is the difference between the initial and remaining fairness,
which are $0.0$, $4.0$, $2.0$ and $1.875$ for the first to fourth individuals, respectively.
It is important to note that beyond the two extreme cases, Shannon entropy deems perturbations
to the third individual (rather than the fourth) create a higher amount of discrimination.

* Min Entropy. The initial fairness via min entropy is also $4$ bits.
The conditional min entropy is $\log \frac{16}{1}=4.0$, $\log \frac{16}{16}=0.0$,
$\log \frac{16}{4}=2.0$, and $\log \frac{16}{5}=1.7$, for the four individuals, respectively.
The amounts of discrimination thus are $0.0$, $4.0$, $\log 4=2.0$ and $\log 5=2.3$, respectively. Beyond the two extreme cases where both entropies agree,
the min entropy deems perturbations to the fourth individual create a higher amount of discrimination. This is intuitive since the discrimination on the fourth case is more subtle, complex, and significant. Therefore, the ML software developers might prioritize those cases characterized by the min entropy.

§ OVERVIEW

Workflow of . Given a DNN and relevant input dataset, quantifies QID discrimination via testing and applies causal debugging to localize and mitigate QID discrimination.

in Nutshell. Figure <ref> shows an overview of our framework .
It consists of two components: (1) an automatic test-generation mechanism based on
search algorithms and (2) a debugging approach that localizes the neurons with significant impacts on fairness
using a causal algorithm. First, searches through the space of input dataset
to find circumstances on the non-protected attributes under which the DNN-under-test shows a significant dependency on
the protected attributes in making decisions.
In doing so, it works in global and local phases. In the
global phase, the search explores the input space to increase the amount of discrimination in each step
of search. On the other hand, in the local phase it exploits the promising seeds from the global phase to generate as many discriminatory instances as possible.

The key elements of search is a threshold-based clustering algorithm used for computing both gradients and objective functions that provide a smooth feedback.
The search characterizes the quantitative individual discrimination (QID)
and returns a set of interesting inputs. Second, uses those inputs to localize neurons with
the largest causal effects on the amounts of discrimination. In doing so, it intervenes [30] over a set of
suspicious neurons.

For every neuron, our debugging approach forces the neurons to be active ($n{>}0$) and non-active ($n{=}0$) over the test cases as far as the functional accuracy of DNN remains in a valid range. Then, it computes the difference between the amounts of QID in these two cases to characterize the causal effects of the neuron on fairness.

reports top $k$ neurons that both have positive impacts (i.e., their activation
reduces the amounts of discrimination) and negative impacts (i.e, their activation
increases the amounts of discrimination).
A potential mitigation strategy is to intervene
to keep a small set of neurons activated (for the positive neurons) or deactivated (for the negative neurons).

Test Cases.
Consider the adult census income [31] dataset with a pre-trained model with $6$ layers [17] to overview in practice. We ran for $1$ hours and obtain $230,593$ test cases. It discovered $36$ clusters
from the initial of $14$ clusters, and the amounts of QID are $4.05$ and $2.64$ bits
for min entropy and Shannon entropy, respectively, out of a total of $5.3$ bits of information from the protected attributes. Considering to order the test cases, we have $6$ test case with maximum QID discrimination of $5.3$ bits. In addition, we have $29$ and $112$ test cases cases with $5.2$ and $5.1$ bits of QID discrimination. The reported numbers are averaged (and rounded) for $10$ runs.

Localization and Mitigation. uses the generated test cases to localize layers and neurons with a significant causal contribution
to the discrimination. For the census dataset, it identifies the second layer as the layer with largest sensitivity to protected attributes.
Among the neurons in this layer,
found that $15$th neuron has the largest negative influence on fairness (the discrimination decreased by $19.6\%$ when it is deactivated) and $19$the neuron has the largest positive influence on fairness (the discrimination decreased by $17.6\%$ when it is activated).
Following this localization, a simple mitigation strategy of activating or deactivating these neuron reduces the amounts of QID discrimination by
$20\%$ with $3\%$ accuracy loss.

Comparison to the State-of-the-art. We compare to the
state-of-the-art techniques in terms of generating individual discrimination (ID) (rather than the quantitative notion)
per each protected attribute. Our goal is to evaluate whether the clustering-based search is effective in generating discriminatory instances.
We run and baseline for $15$ minutes, and report average of results over $10$ runs.
The baseline includes Aequitas [19], ADF [17],
and NeuronFair [18]. Considering sex as the protected attribute in the census dataset, generated $79.0k$ instances whereas Aequitas, ADF, and NeuronFair generated $10.4k$, $18.2k$, and $21.6k$ discriminatory instances, respectively.
Overall, generate more ID instances in all cases with more success rates. However, is slower in finding the first ID
instance in order of a few seconds (in average), since our approach does not generate ID instances in global phase. When considering the time to the first 1,000 instances, has significantly outperformed the state-of-the-art. We conjecture that the improvements are due to smooth search space via quantitative feedback.

§ PROBLEM STATEMENT

We consider DNN-based classifiers with the set of input variables $A$ partitioned into protected set of variables $Z$ (such as race, sex, and age) and non-protected variables $X$ (such as profession, income, and education).
We further assume the output to consist of $t$ prediction classes.

A deep neural network (DNN) encodes a function $\Dd: X \times Z \to [0,1]^t$ where $X = X_1 \times X_2 \cdots \times X_n$ is the set of non-protected input variables, $Z = Z_1 \times Z_2 \cdots \times Z_r$ is the set of protected input variables, and
the output is $t$-dimensional probabilistic vector corresponding to $t$ prediction classes. The predicted label $\Dd_{\tt \ell}(x, z)$ of an input pair $(x, z)$ is the index of the maximum score, i.e. $\Dd_{\tt \ell}(x, z) = \max_{i} \Dd_{\tt \ell}(x, z)(i)$.
We assume that the set of protect input variables are finite domain, and we let $m$ be the cardinality of the set of protected variables $Z$.

A DNN $\Dd$ is parameterized by the input dimension $n{+}r$, the output dimension $t$, the depth of hidden layers $N$,  and the weights of its hidden layers $W_1, W_2, \ldots, W_N$. Our goal is to test and debug a pre-trained neural network with known parameters and weights.
Let $D_{i}$ be the output of layer $i$ that implements an affine mapping from the output of previous layer $D_{i-1}$ and its weights $W_{i-1}$ for $1 \leq i \leq N$ followed by

* a fixed non-linear activation unit (e.g., ReLU defined as $D_{i-1} \mapsto \max \set{W_{i-1}.D_{i-1}, 0}$) for $1 \leq i < N$, or

* a SoftMax function that maps scores to probabilities of each class for $i=N$.
Let $D_i^{j}$ be the output of neuron $j$ at layer $i$.

Individual Discrimination. We say a DNN $\Dd$ is biased based on causal discrimination notion [16, 11, 17, 18] if

\[
\exists z_1, z_2 \in Z, x \in X~s.t.~ \Dd_\ell(x, z_1) \neq \Dd_\ell(x, z_2),
\]

for $z_1 \neq z_2$ of protected inputs. Intuitively, the idea is to find an individual such that their counterfactual with different protected attributes such as race receives a different outcome.

Quantitative Individual Discrimination. In the setting of fairness testing,
it is often desirable to quantify the amounts of bias for individuals. We define the notion quantitative individual discrimination (QID) based on the equivalence classes induced from the output of DNN over protected attributes. Formally, $\qid(Z,X=x) = \seq{Z_1,\ldots,Z_k}$ that is the quotient space of $Z$ characterized by the DNN outputs under an individual with non-protected value $x$. Using this notion, we say a pair of protected values $z,z'$ are in the same equivalence class $i$ (i.e., $z,z' \in Z_i$) if and only if $\Dd(z,x)\approx\Dd(z',x)$.

Given that $Z$ is uniformly distributed
and $\Dd$ is a deterministic function,
we can quantify the QID notion for an individual $(z,x)$ according to the Shannon
and min entropy, respectively:

\[
Q_1(Z,x) {=}  \log_2(m){-}\sum_{i=1}^{k} \frac{|\qid_i(Z,x)|}{m}.\log_2(|\qid_i(Z,x)|)
\]

\[
Q_{\infty}(Z,x) = \log_2(m) - \log_2(\frac{m}{k}) = \log_2(k).
\]

where $m$ is the cardinality of $Z$,
$|\qid_i(Z,x)|$ is the size of equivalence class $i$, and
$k$ is the number of equivalence classes.

Debugging/Mitigating DNN for QID. After characterizing the amounts of discrimination via $\qid$, our next step is to localize a set of layers and neurons that causally effect the output of DNN to have $k$ equivalence classes.

Causal logic [30] provides a firm foundation to reason about the causal relationships between variables. We consider a structural causal model (SCM) with exogenous variables $U$ over the unobserved
input factors, endogenous variables $V$ over $(X,Z,D_i^j)$; and the set of functions $\Ff$ over the set $V$ using the DNN function $\Dd$ and exogenous variables $U$. Using the SCM, we aim to estimate the average causal effect (ACE) [30] of neuron $D_i^j$ on the QID.

A primary tool for performing such computation is called  logic [20]. We write ($i,j,y$) to indicate that the output of neuron $j$ at layer $i$ is intervened to stay $y$. In doing so, we remove the incoming edges to the neuron and force the output of neuron to take a pre-defined value $y$, but we are not required to control back-door variables due to the feed-forward structure of DNN. Then, the ACE of neuron $D_i^j$ on the quantitative individual discrimination with min entropy can be written as $E[Q_{\infty}|$($i,j,y$),$k,l]$,
which is the expected QID after intervening on the neuron given that the non-intervened DNN characterized $k$ classes with an accuracy of $l$. Our goal is to find neurons with the largest causal effects on the QIDs, requiring that such interventions are faithful to the functionality of DNN.

Given a deep neural network model $\Dd$ trained over a dataset $A$ with
protected ($Z \subset A$) and non-protected ($X \subset A$) attributes;
the search problem is to find a single non-protected value $x \in X$
such that the quantitative individual discrimination (QID), for a chosen measure $Q_1$ or $Q_{\infty}$,
is maximized over the $m$ protected values $\Sigma=\{z_1,\ldots,z_m\}$.
Given the inputs $(\Sigma,x)$ characterizing the maximum QID,
our debugging problem is to find a minimal subset of layers
$l \subset \{1,\ldots,N\}$ and neurons $D_{l}^{(J)}$ for $J \subseteq |W_{l}|$
such that the average causal effects of $D_l^{J}$ on the QID
are maximum.

§ APPROACH

Characterizing Quantitative Individual Discrimination.
Given a DNN $\Dd$ over a dataset $A$, our goal is to characterize the worst-case
QID over all possible individuals. Since min entropy characterizes the amounts of discrimination
from one prediction with $Q_{\infty}(Z,x)\geq Q_1(Z,x)$, it is a useful notion to prioritize
the test cases. Therefore, we focus on $Q_{\infty}$ and propose the following objective function:

$\max_{x \in X}~2^{Q_\infty(Z,x)} + (1 - exp(-0.1*\delta))$

where $2^{Q_\infty(Z,x)}=k$ and $\delta$ is the maximum distance between equivalence classes, normalized with the exponential function to remain between $0$ and $1$.
The term is used to break ties when two instances characterize the same number of classes, by preferring one with the highest distance.
Overall, the goal is to find a single value of non-protected attribute $x$ such that the neural network model $\Dd$ predicts many distinguishable
classes of outcomes when $x$ is paired with $m$ protected values.
However, finding those inputs requires an exhaustive search in the exponential set of subsets of input space,
and hence is clearly intractable. We propose a gradient-guided search algorithm that aims to search the space of input variables (attributes) to
maximize the number of equivalence classes and generate as many discrimination instances as possible.

Search Approach.
Our search strategy
consists of global and local phases as in some of the prior work [17, 29, 18].
The goal of the global phase is to find the maximum quantitative individual discrimination via gradient-guided clustering. The local phase uses the promising instances to generate a maximum number of discriminatory instances (ID).

Global Phase.
Given a current instance $x$, the global stage first uses $m$ different values from the space of protected attributes, while keeping the values of non-protected attributes the same.
Then, it receives $m$ prediction scores from the DNN and partitions them into $k$ classes. We adapt a constrained-based clustering with $\epsilon$ where two elements cannot be in the same cluster if their scores differ more than $\epsilon$. Now, the critical step is to perturb the current instance over a subset of non-protected attributes with a direction
that will likely increase the number of clusters induced from the perturbed instance in the next step of global search.

In doing so, we first compute the gradients of DNN loss function for a pair of instances (say $a,a'$) in the cluster with the maximum
elements. The intuition is that we are more likely to split the largest cluster into $2$ or more sub-clusters and increase
the number of partitions in the next step. For the pair of samples, we use the non-protected attributes that
have the same direction of gradients $d$ since it shows the high sensitivity of loss function with respect to
small changes on those common features of the pair.
If we were to use gradients of opposite directions, we will neutralize
the effects of gradients since we only perturb one instance over the non-protected attributes.
Finally, we perturb the current sample $x$ to generate $x'$ using the direction $d$ and step size $s_g$.

Local Phase.
Once we detect an instance with more than $2$ clusters, we enter a local phase where the goal
is to generate as many discriminatory instances (ID) as possible. In our quantitative approach,
we say that an unfavorable decision for an individual $x$ is discriminatory if there is a counterfactual
individual $x'$ that received a favorable outcome. Similar to the state-of-the-art [19, 17, 18], we use non-linear
optimizer that takes an initial instance $x$, a step function to generate the next instance around the neighborhood
of the current instance, and an objective function that quantifies the discrimination of the current instance.

Since our approach uses a continuous objective

based on the characteristics of clusters, it enables us to guide the local search to generate discriminatory instances.

Search Procedure.
Algorithm <ref> sketches our search algorithm to quantify the amounts of bias. We first use the clustering (KMeans algorithm)
to partition the data points into $p$ groups (line 1).
Next, we run the algorithm until the time-out $T$ reaches where in each iteration we seed a sample randomly from one of the partitions $p$ (line 2).
Then, we proceed into global and local phases of

In the global phase, we first use the seed instance $x$ and perturb the protected attributes $P$ to generate $m$
possible instances $X_m$ with different protected values and the same non-protected ones. Then, we run $X_m$
through the DNN model to get the probability scores $S_m$ (line 5). Then, we cluster the scores $S_m$ into
$k$ groups using the tolerance $\epsilon$ (line 6). Afterward, we compute the gradients over two random instances
(from the cluster with the largest size) and use a subset of non-protected features that have the largest number of agreements
on the gradient directions. Then, we use those features and their directions to
perturb the inputs $x$ and generate the next instance $x'$ (line 7-10).
Then, we enter the local phase if the number of clusters or distance between are increased (line 11).

In the local phase, we use the general-purpose optimizer, known as  [32],
which takes an initial seed $x$, an objective function, a step function, and the maximum number of local iterations $N_l$;
it returns the generated instances during the optimization (Line 12-19).
In the objective function shown with  (Line 12-17),
we generate $m$ instances with the same non-protected values but different protected ones (Line 13).
We generate prediction scores for those instances and cluster them with tolerance parameter $\epsilon$ (line 14).
Then, we compute the difference between the indices of two clusters with the smallest and largest scores (line 15).
Finally, we record the generated sample and return the difference as the evaluation of optimizer at the current sample
(line 16-17). The step function is shown with  (Line 18) where it guides the
optimizer to take one step in the input space. Our step function uses a random sample
from a different cluster compared to the current sample. Then, it computes the
normalized sum of gradients and perturbs it using the smallest gradients
to remain in the neighborhood of the current sample.

Debugging Approach. Since it is computationally difficult to intervene over all possible neurons in a DNN, we first adapt a layer-localization technique from the literature of DL framework debugging [33, 21] where we detect a layer with the largest sensitivity to the protected attributes. Let $D_i(z,x)$ be the output of layer $i$ over protected value $z$ and non-protected value $x$.
Let $\Delta_i(x): \Real^{|D_i|} \times \Real^{|D_i|} \to \Real$
be the distance between the outputs of DNN at layer $i$ as triggered by
$m$ different protected values and the same non-protected value $x$,
and let $\delta_i$ be the $\max_{x}\Delta_i(x)$.
The rate of changes in the sensitivity of layer $i$ (w.r.t protected attributes) is

\[
\rho_i = \frac{\delta_i - \max_j\delta_j}{\max_j\delta_j+\epsilon},~with~0\leq~j<~i
\]

where $\delta_0=0.0$ and $\epsilon=10^{-7}$ (to avoid  division-by-zero [33, 21]). Let $l = \argmax_i \rho_i$ be the layer index with the maximum rate of changes. Our next step is to localize neurons in the layer $l$ that have significant positive or negative effects on fairness. Let $V_l^j$ be the set of possible values for neuron $j$ at later $l$ (recorded during the layer localization).
We are interested in computing the average causal effects when the neuron $D_l^j$ is activated vs. deactivated, noting that such interventions might affect the functionality of DNN. Therefore, among a set of intervention values,
we choose one activated value $v_1\in V_l^j > 0$ and one deactivated value $v_2\in V_l^j\approx0$, considering the functional accuracy of DNN within $\epsilon$ of original accuracy $A$.
Therefore, we define average causal difference (ACD) for a neuron $D_l^j$ as:

\[
\eE[Q_{\infty} \mid \texttt{do}(l,j,v_1>0),k,A]{-}\eE[Q_{\infty} \mid \texttt{do}(l,j,v_2\approx0),k,A],
\]

where  notation is used to force the output of neuron $j$ at layer $l$ to a fix value $v$.
We then return the neuron indices with the largest positive (aggravating discrimination) and smallest negative (mitigating discrimination).
Let $\hat{i}$ and $\hat{j}$ be the layer and neuron with the largest positive $ACD$.
One simple mitigation strategy is thus to
deactivate the neuron $D_{l}^{\hat{j}}$, expecting to reduce QID by $ACD/k$ percentage.
Similarly, activating the neuron with the smallest negative $ACD$ is expected
to reduce QID by $ACD/k$ percentage.

Debugging Procedure. Algorithm <ref> shows the debugging aspect of .
Given a set of test cases from the search algorithm, we first use a notion of distance (e.g, $\Delta=L_1$) to compute
the difference between any pair of protected values $z,z'$ w.r.t the outputs of layer $l \in \{1,\ldots,N\}$ (line 1).
Then, we compute the rate of changes (line 2-3) and return a layer $l$ with the largest change (line 4). We compute various
statistics on the output of every neuron $i$ at layer $l$ such as the minimum, maximum, average, average$\pm$std. dev, average$\pm2*$std. dev, etc (line 5).
Among those values, we take the smallest and largest values such that the intervention on the neuron $i$ at layer $l$ has the minimal impacts on the accuracy of DNN (line 6).
Finally, we compute the average causal difference (line 7) and return the indices of layer, neurons with large negative influence, and neurons with large positive influence.

§ EXPERIMENTS

Datasets and DNN models.
We consider $10$ socially critical datasets from the literature of algorithmic fairness.
These datasets and their properties are described in Table <ref>.
For the DNN model, we used the same architecture as the literature [17, 18, 29] and
trained all datasets on a six-layers fully-connected neural network with $\seq{64, 32, 16, 8, 4, 2}$ neurons. We used the same hyperparameters for the all training with num_epochs, batch_size, and learning_rate are set to $1000$, $128$, and $0.01$, respectively.
The accuracy of trained models are reported in Table <ref>.

Technical Details. We implemented with TensorFlow v2.7.0 and scikit-learn v0.22.2.
We run all the experiments on an Ubuntu 20.04.4 LTS OS  sever with AMD Ryzen Threadripper PRO 3955WX 3.9GHz 16-cores X 32 CPU and two NVIDIA GeForce RTX 3090 GPUs. We choose the values $10$, $1000$, $0.025$, $1$, and $1$ for max_global, max_local, $\epsilon$, $s\_g$ , and $s\_l$ in Algorithm <ref>, respectively, and take the average of $10$ multiple runs for all experiments. In Algorithm <ref>, we used $L_1$-norm, $10^{-7}$, $0.05$, and $3$ for $\Delta$, $\epsilon_1$, $\epsilon_2$, and $k$, respectively.

Research Questions.
We seek to answer the following three questions using our experimental setup.

* Can characterize the amounts of information from protected attributes used for the inferences?

* Is the the proposed search algorithm effective and efficient (vis-a-vis the
state-of-the-art techniques) in generating individual discrimination instances?

* Can the proposed causal debugging guide us to localize and mitigate the amounts of discrimination?

Our open-source tool with all experimental subjects are
publicly accessible:

* <https://github.com/armanunix/Fairness-testing>

§.§ Characterizing QID via Search (RQ1)

An important goal is to characterize the amount of information from protected attributes used during the inference of DNN models. Table <ref> shows the result of
experiments to answer this research question. The left side of table shows the initial characteristics such as the number of protected values ($m$), the maximum possible amounts of discrimination ($Q_I$) based on ($\epsilon^{-1},m$), and the initial number of clusters found using samples from the dataset ($K_I$).
The right side of table shows the results after running our search for $1$ hour. The column #$I$ is the number of QID instances generated, and $K_{F}$ is the maximum number of clusters discovered by . The column T$_{K_{F}}$ is the time taken to find the maximum number of clusters from an input with initial clusters $K_{I}$ (in seconds). The columns $Q_{\infty}$ and $Q_{1}$
are the quantitative individual discrimination based on min entropy and Shanon entropy, respectively.
The columns #$I_{K_{F}^1}$, #$I_{K_{F}^2}$, and #$I_{K_{F}^3}$ show the number of test cases with the highest, second-highest, and third-highest QIDs, respectively, that order test cases with their QID severity.
Overall, the results show that can find $3.4\times$ more clusters (in average) from the initial characteristics within one minute of search. The DNN for Students dataset showed the largest increase in the number of clusters going from $1.9$ to $10.9$. found that Adult Income Census dataset has the largest amounts of QID where $4.05$ out of $5.3$ bits ($76.4\%$) from protected variables are used to make decisions. The German Credit dataset with $1.61$ out of $4.0$ bits ($40.0\%$) showed the least amounts of discrimination. For test-case prioritizing, the column #$I_{K_{F}^1}$ shows our approach to be useful in finding a small percentage of generated test cases with the worst-case discrimination. In $7$ out of $10$ experiments,
found less than $50$ test cases with severe discrimination out of hundreds of thousands inputs.

Answer RQ1:
The search algorithm is effective in characterizing the
amounts of discrimination via QID. Within $1$ hour, it increased the number of clusters by $3.4\times$ in average, and found instances that used up to
$76\%$ of protected information ($4.05$ out of $5.3$ bits) to infer DNN outcomes.
is useful to prioritize test cases with their severity where it generates less than $50$ test cases with the maximum QID among hundreds of thousands test cases.

§.§ Individual Discriminatory Instances (RQ2)

In this section, we compare the efficiency and effectiveness of our search algorithm to the state-of-the-art techniques in searching individual discrimination ($ID$) instances (as defined in Section <ref>).
Our baselines are
Aequitas [19], ADF [17], and NeuronFair [18].
We obtained the implementations of these tools from their GitHub repositories and configured them according to the prescribed setting to have the best performance. Following these techniques, we report the results for each protected attribute separately.
Table <ref> shows the results of baselines and in runs of 15 minutes. The results are averaged over $10$ repeated runs. The column #$ID$ is the total number of generated individual discriminatory instances. The column  is the success rate of local stage of searches. We exclude the global success rate since the goal of global phase in our search is to maximize QID whereas the local phase focuses to generate many $ID$ instances.
We calculate success rate as the number of $ID$ found over the total number of generated samples. The columns  and   are the amount of time (in seconds) taken to find the first $ID$ instance and to generate $1,000$ individual discriminatory instances (note: $N/A$ in column  means that the tool did not generate $1,000$ $IDs$ in the experiment timeout of $900$ seconds in the average of 10 runs).

The result shows that outperforms the-state-of-the-art in generating many ID instances. In particular, finds $27.1\times$, $16.0\times$, and $16.0\times$ more $ID$s in the best case
compare to Aequitas, ADF, and NeuronFair, respectively.
also generates $3.2\times$, $2.3\times$, and $2.6\times$ more $ID$s in the worst case compare to Aequitas, ADF, and NeuronFair, respectively.
The success rate of local search are $20.6\%$, $33.0\%$, $29.6\%$, and $78.2\%$ in average for Aequitas, ADF, NeuronFair, and , respectively. For the time taken to find the first $ID$, Aequitas achieves the best result with an average of $0.03$ (s) whereas it took $1.46$ (s) in average to find the first $ID$. In average, was found to take the lowest time to generate $1000$ $IDs$ with $57.2$ (S), while ADF took $179.1$ (s), NeuronFair took $135.4$ (s), and Aequitas took the longest time at $197.7$ (s).
Overall, our experiments indicate that is effective in generating $ID$ instances compared to the three state-of-the-art techniques, largely due to the smoothness of the feedback during the local search.

Answer RQ2:
Our experiments demonstrate that outperforms the state-of-the-art fairness testing techniques [19, 17, 18]. In the best case, our approach found $20\times$
more individual discrimination (ID) instances than these techniques with almost $3\times$ more success rates in average.
However, we found that is slower than those techniques in finding the first ID instance in the order of a few seconds.

§.§ Causal Debugging of DNNs for Fairness (RQ3)

We perform experiments over the DNN models to study whether
the proposed causal debugging approach is useful in identifying layers and neurons that significantly effect the amounts of discrimination as characterized by $QID$. Table <ref> shows the results of experiments (averaged over $10$ independent runs). The first two columns show the localized layer and its influence (i.e., $l$ and $\rho$ in Algorithm <ref>). The next six columns show top $3$ neurons with the positive influence on fairness (i.e., activating those neurons reduce the amounts of discrimination based $Q_{\infty}$). The last six columns show top $3$ neurons with the negative influence on fairness (i.e., activating those neurons increase the amounts of discrimination based $Q_{\infty}$).
The layer index $2$ is more frequently localized than other layers where the layers $3$, $4$, and $5$ are localized once.
Overall, the average causal difference (ACD) ranges from $4\%$ to $55\%$ for neurons with positive fairness effects
and from $0.6\%$ to $18.3\%$ for neurons with negative fairness effects.

Guided by localization, intervenes to activate neurons with positive fairness influence or de-activate
those with negative influence.
Table <ref> shows the results of this mitigation strategy. The columns $A$ and $K$ show the accuracy and the number of clusters (averaged over a set of random test cases) reported by before mitigation over the DNN model. The columns $A^{=0}$
and $K^{=0}$ are accuracy and the number of clusters reported after mitigating the DNN model by de-activating the neuron with the highest negative fairness impacts (as suggested by Neuron$^{-}_{1}$ in Table <ref>). Similarly, the columns $A^{>0}$
and $K^{>0}$ are accuracy and the number of clusters reported after mitigating the DNN model by activating the neuron with the highest positive fairness impacts (as suggested by Neuron$^{+}_{1}$ in Table <ref>).
The results indicate that the activation interventions
can reduce QID discrimination by at least $5\%$ with $3\%$ loss of accuracy and up to $64.3\%$ with $2\%$ loss of accuracy.
The de-activation, on the other hand, can improve the fairness by at least $6\%$ with $1\%$ loss of accuracy and up to $27\%$ with $2\%$ loss.

Answer RQ3:
The debugging approach implemented in identified neurons that
have at least $5\%$ and up to $55\%$ positive causal effects on the fairness and those which have at least $0.6\%$ and up to $18.3\%$
negative causal effects. A mitigation strategy followed by the localization can reduce the amounts of discrimination by at least $6\%$ and up to $64.3\%$ with less than 5% loss of accuracy.

§ DISCUSSION

Limitation. In this work, we consider all set of protected values
and perturb them to generate counterfactual. Various perturbations of protected attributes
may yield unrealistic counterfactuals and contribute towards false positives
(an over-approximation of discrimination). This limitation can be mitigated by supplying domain-specific constraints (Age$<$YY$\implies$NOT(married)): we already apply some common-sense constraints
(e.g., to ensure valid range of age). In addition, similar to any dynamic testing methods,
our approach might miss discriminatory inputs and is prone to false negatives. The probability of missing relevant inputs can be contained under a suitable statistical testing (e.g., Bayes factor). In addition,
our debugging approach is similar to pin-pointing suspicious code fragments and is based on causal reasoning of its effect in decision making rather than correlation. But, it is not to furnish explanations or interpretations of black-box DNN functions.

Threat to Validity. To address the internal validity and ensure
our finding does not lead to invalid conclusion, we follow established
guideline and take average of repeated experiments.
To ensure that our results are generalizable and address external validity,
we perform our experiments on $10$ DNN models taken from the literature of fairness testing.
However, it is an open problem whether these datasets and DNN models are sufficiently
representative for fairness testing.

§ RELATED WORK

Fairness Testing of ML systems.
Themis [10] presents a causal discrimination notion
where they measure the difference between the fairness metric of two subgroups
by counterfactual queries; i.e., they sample individuals with the protected attributes
set to A and compare the outcome to a counterfactual scenario where the protected attributes
are set to B. Symbolic generation (SC) [28, 22]
presents a black-box testing that approximates the ML models with decision trees and
leverage symbolic execution over the tree structure
to find individual discrimination (ID).
AEQUITAS [19] uses a two-step approach
that first uniformly at random
samples instances from the input dataset to find a discriminatory instance and then
locally perturb those instances to further generate biased test cases.
ExpGA [16] proposed a genetic algorithm (GA)
to generate ID instances in natural language processes. The proposed technique
used a prior knowledge graph to guide the perturbation of protected attributes in the NLP tasks.
While these techniques are black-box, they potentially suffer
from the lack of local guidance during the search.
ADF [17] utilized the gradient of the loss function as guidance in generating ID instances.
The global phase explores the input space to find diverse set of individual discrimination
whereas the local phase exploits each instance to generate many individual discriminatory (ID) instances
in their neighborhoods. EIDIG [29] follows similar
ideas to ADF, but uses different computations of gradients.
First, it uses the gradients of output (rather than
loss function) to reduce the computation cost at each iteration. Second, it uses momentum
of gradients in global phase to avoid local optima.
NeuronFair [18] extends ADF and EIDIG to support
unstructured data (e.g., image, text, speech, etc.) where the protected attributes
might not be well-defined. In addition, NeuronFair is guided by the DNN's
internal neuron states (e.g., the pattern of activation and deactivation) and their
activation difference. Beyond the capability of these techniques, quantifies the amounts of discrimination,
enables software developers to prioritize test cases, and searches multiple protected attributes at one time.

Beyond the scope of this paper, a body of prior work [42, 43, 44, 23, 45, 46] considered testing for group fairness. Fairway [43] mitigates biases
after finding suitable ML algorithm configurations. In doing so, they used
a multi-objective optimization (FLASH) [47].
Parfait-ML [45] searches the hyperparameter
space of classic ML algorithms via a gray-box evolutionary algorithm to characterize
the magnitude of biases from the hyperparameter configuration.

Debugging of Deep Neural Network.
Cradle [33] traced the execution graph of
a DNN model over two different deep-learning frameworks and used the differences
in the outcomes to localize what backend functions might cause
a bug. However, since Cradle did not use causal analysis,
it showed a high rate of false positive. Audee [48]
used a similar approach, but it leveraged causal-testing
methods. In particular, it designed strategies to intervene
in the DNN models and tracked how the intervention affected the observed
inconsistencies. We adapted the layer localization of Cradle
and Audee; but our causal localization is developed using
logic for a meta-property (fairness). Audee
used a simple perturbation of neuron values for
functional correctness (i.e., any inconsistency shows a bug) without
considering the accuracy or the severity of neuron contributions to a bug.

In-process Mitigation.
A set of work considers in-process algorithms to mitigate biases in ML
predictions [49, 50, 51]. Adversarial debiasing [49] and
Prejudice remover [50] improve fairness by adding constraints to model parameters or the loss function. Exponentiated gradient [51] uses a meta-learning algorithm to infer a family of classifiers that maximizes accuracy and fairness. Different than these approaches, we develop a mitigation approach
that is specialized to handle neural networks for individual fairness. This setting allows us to exploit the layer-based structure of NNs toward causal reasoning and mitigation. We believe that our approach can be extended with in-process mitigation techniques to maximize fairness in the DNN-based decision support systems.

Formal Methods. We believe that this paper can connect to the rich literature of formal verification and its application. Here, we provide two examples. FairSquare [52] certifies
a fair decision-making process in probabilistic programs using a novel verification
technique called the weighted-volume-computation algorithm. SFTREE [53] formulated the problem of inferring fair decision tree as a mixed integer linear programming and apply constraint solvers iteratively to find solutions.

Fairness in income, wealth, and taxation. We develop a fairness testing and debugging approach that is uniquely geared toward handling regression problems. Therefore, our approach can be useful to study and address biases in income and wealth distributions [54] among different race and gender. Furthermore, our approach can be useful to study fairness in taxation (e.g., vertical and horizontal equities [55, 56]). We left further study in these directions to future work.

§ CONCLUSION

DNN-based software solutions are increasingly being used in socio-critical applications where a bug in their design may lead to discriminatory behavior.
In this paper, we presented : an information-theoretic model to characterize the amounts of protected information used in DNN-based decision making.
Our experiments showed that the search and debugging algorithms, based on the quantitative landscape, are effective in discovering and localizing fairness defects.

[1]
I. J. Goodfellow, Y. Bengio, and A. Courville, Deep Learning.1em
plus 0.5em minus 0.4emCambridge, MA, USA: MIT Press, 2016,

[2]
S. Ranchordás and L. Scarcella, “Automated government for vulnerable
citizens: Intermediating rights,” SSRN Electronic Journal, 2021.

[3]
N. Tishby, F. C. Pereira, and W. Bialek, “The information bottleneck method,”
arXiv preprint physics/0004057, 2000.

[4]
S. M. Julia Angwin, Jeff Larson and L. Kirchne, “Machine bias,”
2021, online.

[5]
D. A. Elyounes, “" computer says no!": The impact of automation on the
discretionary power of public officers,” Vand. J. Ent. & Tech. L.,
vol. 23, p. 451, 2020.

[6]
N. Escher and N. Banovic, “Exposing error in poverty management technology:
A method for auditing government benefits screening tools,” Proc.
ACM Hum. Comput. Interact., vol. 4, no. CSCW, pp. 064:1–064:20, 2020.
[Online]. Available: <https://doi.org/10.1145/3392874>

[7]
J. Slemrod, “Group equity and implicit discrimination in tax systems,”
National Tax Journal, vol. 75, no. 1, pp. 201–224, 2022.

[8]
D. A. Brown, “The IRS is targeting the poorest americans,” August 2021,
[Online; posted 27-July-2021]. [Online]. Available:

[9]
C. Thomas and A. Pontón-Núñez, “Automating judicial discretion:
How algorithmic risk assessments in pretrial adjudications violate equal
protection rights on the basis of race,” Minnesota Journal of Law &
Inequality, vol. 40, no. 2, p. 371, 2022.

[10]
R. Angell, B. Johnson, Y. Brun, and A. Meliou, “Themis: Automatically testing
software for discrimination,” in Proceedings of the 2018 26th ACM
Joint meeting on european software engineering conference and symposium on
the foundations of software engineering, 2018, pp. 871–875.

[11]
S. Galhotra, Y. Brun, and A. Meliou, “Fairness testing: testing software for
discrimination,” in FSE, 2017, pp. 498–510.

[12]
A. Aggarwal, P. Lohia, S. Nagar, K. Dey, and D. Saha, “Black box fairness
testing of machine learning models,” in Proceedings of the 2019 27th
ACM Joint Meeting on European Software Engineering Conference and Symposium
on the Foundations of Software Engineering, 2019, pp. 625–635.

[13]
A. Sharma and H. Wehrheim, “Automatic fairness testing of machine learning
models,” in IFIP International Conference on Testing Software and
Systems.1em plus 0.5em minus 0.4emSpringer, 2020, pp.

[14]
M. Hardt, E. Price, and N. Srebro, “Equality of opportunity in supervised
learning,” in NIPS, 2016.

[15]
C. Dwork, M. Hardt, T. Pitassi, O. Reingold, and R. Zemel, “Fairness through
awareness,” in Proceedings of the 3rd innovations in theoretical
computer science conference, 2012, pp. 214–226.

[16]
M. Fan, W. Wei, W. Jin, Z. Yang, and T. Liu, “Explanation-guided fairness
testing through genetic algorithm,” in Proceedings of the 44th
International Conference on Software Engineering, ser. ICSE '22.1em
plus 0.5em minus 0.4emNew York, NY, USA: Association for Computing
Machinery, 2022, p. 871–882. [Online]. Available:

[17]
P. Zhang, J. Wang, J. Sun, G. Dong, X. Wang, X. Wang, J. S. Dong, and T. Dai,
“White-box fairness testing through adversarial sampling,” in
Proceedings of the ACM/IEEE 42nd International Conference on Software
Engineering, 2020, pp. 949–960.

[18]
H. Zheng, Z. Chen, T. Du, X. Zhang, Y. Cheng, S. Ti, J. Wang, Y. Yu, and
J. Chen, “Neuronfair: Interpretable white-box fairness testing through
biased neuron identification,” in 2022 IEEE/ACM 44th International
Conference on Software Engineering (ICSE), 2022, pp. 1519–1531.

[19]
S. Udeshi, P. Arora, and S. Chattopadhyay, “Automated directed fairness
testing,” in Proceedings of the 33rd ACM/IEEE International Conference
on Automated Software Engineering, 2018, pp. 98–108.

[20]
M. Glymour, J. Pearl, and N. P. Jewell, Causal inference in statistics: A
primer.1em plus 0.5em minus 0.4emJohn Wiley & Sons, 2016.

[21]
Q. Guo, X. Xie, Y. Li, X. Zhang, Y. Liu, X. Li, and C. Shen, “Audee: Automated
testing for deep learning frameworks,” in 2020 35th IEEE/ACM
International Conference on Automated Software Engineering (ASE).1em
plus 0.5em minus 0.4emIEEE, 2020, pp. 486–498.

[22]
A. Agarwal, P. Lohia, S. Nagar, K. Dey, and D. Saha, “Automated test
generation to detect individual discrimination in ai models,” arXiv
preprint arXiv:1809.03260, 2018.

[23]
J. Chakraborty, S. Majumder, and T. Menzies, “Bias in machine learning
software: Why? how? what to do?” in Proceedings of the 29th ACM Joint
Meeting on European Software Engineering Conference and Symposium on the
Foundations of Software Engineering, ser. ESEC/FSE 2021.1em plus
0.5em minus 0.4emNew York, NY, USA: Association for Computing
Machinery, 2021, p. 429–440. [Online]. Available:

[24]
A. Rényi et al., “On measures of entropy and information,” in
Proceedings of the fourth Berkeley symposium on mathematical statistics
and probability, vol. 1, no. 547-561.1em plus 0.5em minus 0.4emBerkeley, California, USA, 1961.

[25]
G. Smith, “On the foundations of quantitative information flow,” in
Foundations of Software Science and Computational Structures,
L. de Alfaro, Ed.1em plus 0.5em minus 0.4emBerlin, Heidelberg:
Springer Berlin Heidelberg, 2009, pp. 288–302.

[26]
M. Backes, B. Köpf, and A. Rybalchenko, “Automatic discovery and
quantification of information leaks,” in 2009 30th IEEE Symposium on
Security and Privacy.1em plus 0.5em minus 0.4emIEEE, 2009, pp.

[27]
Y. Noller and S. Tizpaz-Niari, “Qfuzz: Quantitative fuzzing for side
channels,” ser. ISSTA 2021.1em plus 0.5em minus 0.4emNew York,
NY, USA: Association for Computing Machinery, 2021, p. 257–269. [Online].
Available: <https://doi.org/10.1145/3460319.3464817>

[28]
A. Aggarwal, P. Lohia, S. Nagar, K. Dey, and D. Saha, “Black box fairness
testing of machine learning models,” in Proceedings of the 2019 27th
ACM Joint Meeting on European Software Engineering Conference and Symposium
on the Foundations of Software Engineering, ser. ESEC/FSE 2019, 2019, p.
625–635. [Online]. Available: <https://doi.org/10.1145/3338906.3338937>

[29]
L. Zhang, Y. Zhang, and M. Zhang, “Efficient white-box fairness testing
through gradient search,” in Proceedings of the 30th ACM SIGSOFT
International Symposium on Software Testing and Analysis, ser. ISSTA 2021,
2021, p. 103–114. [Online]. Available:

[30]
J. Pearl, Causality.1em plus 0.5em minus 0.4emCambridge
university press, 2009.

[31]
D. Dua and C. Graff, “UCI machine learning repository,” 2017. [Online].
Available: <https://archive.ics.uci.edu/ml/datasets/census+income>

[32]
D. C. Liu and J. Nocedal, “On the limited memory bfgs method for large scale
optimization,” Mathematical programming, vol. 45, no. 1, pp.
503–528, 1989.

[33]
H. V. Pham, T. Lutellier, W. Qi, and L. Tan, “Cradle: cross-backend validation
to detect and localize bugs in deep learning libraries,” in 2019
IEEE/ACM 41st International Conference on Software Engineering (ICSE).1em plus 0.5em minus 0.4emIEEE, 2019, pp. 1027–1038.

[34]
ProPublica, “Compas software ananlysis,”
<https://github.com/propublica/compas-analysis>, 2021, online.

[35]
“UCI:statlog (german credit data) data set,” 2000. [Online]. Available:

[36]
“UCI:default of credit card clients data set,” 2009. [Online]. Available:

[37]
“UCI:heart disease data set,” 2001. [Online]. Available:

[38]
“Bank marketing uci,” 2017. [Online]. Available:

[39]
“UCI: Diabetes patient records,” 1994. [Online]. Available:

[40]
“Student performance data set,” 2014. [Online]. Available:

[41]
“Medical expenditure panel survey,” 2014. [Online]. Available:

[42]
R. K. Bellamy, K. Dey, M. Hind, S. C. Hoffman, S. Houde, K. Kannan, P. Lohia,
J. Martino, S. Mehta, A. Mojsilović et al., “AI fairness 360:
An extensible toolkit for detecting and mitigating algorithmic bias,”
IBM Journal of Research and Development, vol. 63, no. 4/5, pp. 4–1,

[43]
J. Chakraborty, S. Majumder, Z. Yu, and T. Menzies, “Fairway: a way to build
fair ml software,” in Proceedings of the 28th ACM Joint Meeting on
European Software Engineering Conference and Symposium on the Foundations of
Software Engineering, 2020, pp. 654–665.

[44]
J. M. Zhang and M. Harman, “"ignorance and prejudice" in software fairness,”
in 43rd IEEE/ACM International Conference on Software Engineering,
ICSE 2021, Madrid, Spain, 22-30 May 2021.1em plus 0.5em minus
0.4emIEEE, 2021, pp. 1436–1447. [Online]. Available:

[45]
S. Tizpaz-Niari, A. Kumar, G. Tan, and A. Trivedi, “Fairness-aware
configuration of machine learning libraries,” in 44th IEEE/ACM 44th
International Conference on Software Engineering, ICSE 2022, Pittsburgh,
PA, USA, May 25-27, 2022.1em plus 0.5em minus 0.4emIEEE,
2022, pp. 909–920. [Online]. Available:

[46]
Z. Chen, J. M. Zhang, F. Sarro, and M. Harman, “Maat: A novel ensemble
approach to addressing fairness and performance bugs for machine learning
software,” ser. ESEC/FSE 2022, 2022, p. 1122–1134. [Online]. Available:

[47]
V. Nair, Z. Yu, T. Menzies, N. Siegmund, and S. Apel, “Finding faster
configurations using flash,” IEEE Transactions on Software
Engineering, vol. 46, no. 7, pp. 794–811, 2020.

[48]
Q. Guo, X. Xie, Y. Li, X. Zhang, Y. Liu, X. Li, and C. Shen, “Audee: Automated
testing for deep learning frameworks,” in 2020 35th IEEE/ACM
International Conference on Automated Software Engineering (ASE), 2020, pp.

[49]
B. H. Zhang, B. Lemoine, and M. Mitchell, “Mitigating unwanted biases with
adversarial learning,” in Proceedings of the 2018 AAAI/ACM Conference
on AI, Ethics, and Society, 2018, pp. 335–340.

[50]
F. Kamiran, A. Karim, and X. Zhang, “Decision theory for discrimination-aware
classification,” in 2012 IEEE 12th International Conference on Data
Mining, 2012, pp. 924–929.

[51]
A. Agarwal, A. Beygelzimer, M. Dudík, J. Langford, and H. Wallach, “A
reductions approach to fair classification,” in International
Conference on Machine Learning.1em plus 0.5em minus 0.4emPMLR,
2018, pp. 60–69.

[52]
A. Albarghouthi, L. D'Antoni, S. Drews, and A. V. Nori, “Fairsquare:
Probabilistic verification of program fairness,” no. OOPSLA, 2017. [Online].
Available: <https://doi.org/10.1145/3133904>

[53]
J. Wang, Y. Li, and C. Wang, “Synthesizing fair decision trees via iterative
constraint solving,” in Computer Aided Verification: 34th
International Conference, CAV 2022, Haifa, Israel, August 7–10, 2022,
Proceedings, Part II.1em plus 0.5em minus 0.4emSpringer, 2022,
pp. 364–385.

[54]
T. Piketty and E. Saez, “Income inequality in the united states, 1913–1998,”
The Quarterly journal of economics, vol. 118, no. 1, pp. 1–41, 2003.

[55]
E. Black, H. Elzayn, A. Chouldechova, J. Goldin, and D. Ho, “Algorithmic
fairness and vertical equity: Income fairness with irs tax audit models,” in
2022 ACM Conference on Fairness, Accountability, and Transparency,
ser. FAccT '22.1em plus 0.5em minus 0.4emNew York, NY, USA:
Association for Computing Machinery, 2022, p. 1479–1503. [Online].
Available: <https://doi.org/10.1145/3531146.3533204>

[56]
S. Tizpaz-Niari, V. Monjezi, M. Wagner, S. Darian, K. Reed, and A. Trivedi,
“Metamorphic testing and debugging of tax preparation software,” 2023.
[Online]. Available: <https://arxiv.org/abs/2205.04998>


11institutetext: Institute of Nuclear Physics, ul. Radzikowskiego 152, 31–342
Cracow, Poland

# TiC Lattice Dynamics from ab initio Calculations

P.T. Jochym  K. Parlinski  M. Sternik<EMAIL_ADDRESS>
###### Abstract

Ab initio calculations and a direct method have been applied to derive the
phonon dispersion curves and phonon density of states for the TiC crystal. The
results are compared and found to be in a good agreement with the experimental
neutron scattering data. The force constants have been determined from the
Hellmann-Feynman forces induced by atomic displacements in the
$2$$\times$$2$$\times$$2$ supercell. The calculated phonon density of states
suggests that vibrations of Ti atoms form acoustic branches, whereas the
motion of C atoms is confined to optic branches. The elastic constants have
been found using the deformation method and compared with the results obtained
from acoustic phonon slopes.

###### pacs:

63.20.–ePhonons and vibrations in crystal lattices and 71.15.NTotal energy
calculations

The transition-metal carbide compounds, of which TiC is a representative, are
of big scientific and technological interest because of their striking
mechanical properties, extreme hardness combined additionally with metallic
electrical and thermal conductivity. They have usually a rock-salt crystal
structure – like some ionic crystals – while having properties typical for
materials with covalent bonding. It suggests that the bonding have some mixed
nature, which makes these materials very interesting as an object of study.

Over a large range of fractional carbon content from TiC to TiC0.5 toth  the
TiC is stable in the NaCl structure, with space group Fm${\bar{3}}$m. In
practice, it is usually nonstoichiometric and contains carbon vacancy defects.
The phonon dispersion curves of nearly stoichiometric TiC0.95 have been
measured along main symmetry directions by Pintschovius et al pint  using the
inelastic neutron scattering technique. In the same paper, the influence of
the carbon content on the phonon dispersion curves has been additionally
checked by measuring a sample with lower carbon concentration TiC0.89. The
results shown that the optic modes could differ by 3.5%, while the acoustic
modes are less sensitive to carbon content. The elastic constants derived from
the measured acoustic dispersion curves agree very well with those obtained by
ultrasonic measurements for TiC0.91 crystal exp_dataA .

Very interesting properties and simple structure of TiC caused that efforts to
understand its properties using first-principle calculations have been
undertaken several times neckel ; zhukov ; price ; haglund . The electronic
structure, bulk modulus, and elastic constants have been found TiC_GGA  by
means of accurate first-principles total-energy calculations using the full
potential linear muffin-tin orbital method. The calculated values are
generally in very good agreement with experiment.

In this paper we intend to extend first-principle calculations to describe the
phonon dispersion curves and phonon density of TiC. The method which was is
based on the total energy calculation and Hellmann–Feynman (HF) forces.
Dispersion relations are calculated by direct method frank ; kresse ;
parlinski ; ackland ; sluiter , in which the force constants of the dynamical
matrix are calculated from HF forces. Elastic constants and bulk modulus have
been estimated by straight evaluation of energy derivatives with respect to
deformation.

The energies and forces of the TiC crystal were calculated by the method of
total energy minimization, using a norm–conserving pseudopotentials as an
approximation of the atomic core – valence electron interaction lin ;
goniakowskiA . For an excellent review of ab initio total–energy calculations
see Ref.payne  and further references given there. For the lattice dynamics
calculations the $2$$\times$$2$$\times$$2$ supercell with periodic boundary
conditions and 64 atoms has been used. For optimizing the structure and for
direct calculation of energy deformations $1$$\times$$1$$\times$$1$ supercell
has been utilized. The ab initio total energy calculations were done with the
CASTEP package castep  and standard pseudopotentials provided within this
package. These pseudopotentials were constructed within LDA approximation. The
Ti pseudopotential treats 3p electrons as belonging to the core. We have
employed non–local variant of those potentials parametrized in the reciprocal
space with 900 eV and 600 eV cut–off energy for $1$$\times$$1$$\times$$1$ and
$2$$\times$$2$$\times$$2$ supercell, respectively. Tests which had been made
with the cut-off energy of 600 eV applied to the $1$$\times$$1$$\times$$1$
supercell, showed that in this case the cohesive energy is only 0.2 eV higher,
and the lattice constant is longer by 0.0006Å(0.01%), comparing to the 900 eV
cut-off. Therefore, we have used 600 eV cut-off energy for a larger supercell
calculation.

A generalized gradient approximation (GGA) has been used for the exchange
energy term of the valence states of the Hamiltonian perdew . Although this
procedure is not quite consistent, there are indications goniakowskiB  that
LDA generated pseudopotentials used with GGA exchange term in crystals are
sufficiently accurate for energies and structures. The integration over the
Brillouin zone has been performed with weighted summation over wave vectors
generated by Monkhorst-Pack scheme pack  using 0.1 Å-1 and 0.06 Å-1 grid sizes
which leads to 4 and 32 wave vectors for $1$$\times$$1$$\times$$1$ supercell
and to 1 and 4 wave vectors for $2$$\times$$2$$\times$$2$ supercell,
respectively. No significant difference due to different number of wave
vectors, 1 or 4, has been found in calculated phonon frequencies. The results
of dispersion curve presented here are generated from the 0.06 Å-1 data set
with 4 wave vectors.

Quantity | Exper. | Non-metallic |  | Metallic |
---|---|---|---|---|---
# k |  | 4 | 4 | 14 | 32 | 63
$a_{0}$ | $4.327$ | $4.345$ | $4.35$ | $4.342$ | $4.345$ | $4.342$
$B$ | $2.4$ | $2.42$ | $2.41$ | $3.527$ | $2.418$ | $2.427$
$\omega_{\Gamma}$ | $16.3$ | $16.5$ | $19.0$ | $11.8$ | $15.6$ | $15.6$
$\omega_{X}^{1}$ | $8.5$ | $8.5$ | $10.2$ | $5.1$ | $9.2$ | $8.5$
$\omega_{X}^{2}$ | $10.8$ | $10.0$ | $11.1$ | $10.3$ | $10.7$ | $10.7$
$\omega_{X}^{3}$ | $16.5$ | $17.0$ | $18.2$ | $12.8$ | $15.7$ | $15.3$
$\omega_{X}^{4}$ | $18.5$ | $17.8$ | $19.3$ | $16.7$ | $18.1$ | $17.7$
Table 1: Comparison of results from metallic and non-metallic approach

Since TiC is a weakly metallic compound it could have been appropriate to use
a smearing of electron levels in the Brillouin zone integration. This
procedure is computationally more expensive. Thus we have performed a series
of tests to compare lattice constants, bulk moduli and phonons in $\Gamma$ and
X points obtained using gaussian smearing in the range from 4 to 0.1 eV and
denser k-space grids (4, 14, 32, 63 k-points) for one unit cell configuration.
The comparison of results shown in the Table 1 indicate that, for sufficiently
dense k-space grid used, these quantities are not very sensitive (specially
the lattice constant and bulk modulus) to the choice of method of Brillouin
zone integration. Thus we have decided to use much simpler and faster non-
metallic approach in the final $2$$\times$$2$$\times$$2$ supercell
calculation. Note also that results obtained in the metallic regime lead to
larger divergency of phonon frequencies from experimental data.

The minimization of the total energy with respect to the lattice constant was
performed within the CASTEP minimizer module with $1$$\times$$1$$\times$$1$
supercell and cross-checked with $2$$\times$$2$$\times$$2$ supercell. The
equilibrium lattice constant is $a=4.3448$Å which could be compared with the
experimental value $a=4.3269$Ånaray . Such a good agreament may be
coincidental, since inclusion of 3p electrons in the valence band as well as
gradient corrections in core states tends to change the equilibrium lattice
constant garcia ; juan ; seifert .

The HF forces are defined as

${\bf F}_{i}({\bf n},\nu)=-\partial E_{tot}/\partial{\bf R}_{i}({\bf n},\nu)$
(1)

where ${\bf n},\nu$ are the indices of the unit cell and atom within the unit
cell, respectively, and $i$ is the Cartesian component. At extremum all HF
forces vanish. Non-zeroth HF forces arise, when a single atom $({\bf m},\mu)$
is displaced by $u_{j}({\bf m},\mu)$ from its equilibrium position. The
arising forces are related with the cummulant force constants $\Phi_{ij}$ by a
relation parlinski ; ackland

${\bf F}_{i}({\bf n},\nu)=-\sum_{{\bf m},\mu,j}\Phi_{ij}({\bf n},\nu;{\bf
m},\mu)u_{j}({\bf m},\mu)$ (2)

Cummulant force constants appear as a result of periodic boundary conditions
imposed on the supercell. To calculate HF forces an atomic configuration with
a single displaced atom must be minimized with respect to the electronic part
only. Each of such runs provides $3n$ HF forces, where $n=8$ or $64$ is the
number of atoms in the $1$$\times$$1$$\times$$1$ or $2$$\times$$2$$\times$$2$
supercell, respectively. Two runs with displacements along $z$ direction, one
for Ti and second for C were performed. To minimize the systematic errors two
other runs with negative displacements were carried on. The displacement
amplitudes were set to $1\%$ of the lattice constant. This amplitude has been
selected after number of careful tests in which displacements ware varied from
$0.1\%$ to $1.6\%$. It has been proven that the anharmonic contributions are
negligable up to $1\%$ of displacement. At smaller displacements, close to
$0.1\%$ the HF forces became too small and the uncertainty of force constant
evaluation is not acceptable.

Figure 1: Absolute eigenvalues of the force constant matrices of the TiC
crystal as a function of interatomic distance.

The data of $4\times 3n$ HF forces ${\bf F}_{i}({\bf n},\nu)$ and four
displacements $u_{j}({\bf m},\mu)$ form an overdetermined set of equations,
Eq.(2), for the force constants. This system was solved by the singular value
decomposition algorithm parlinski ; numrec , which automatically provides the
least-squares solution. For $2$$\times$$2$$\times$$2$ supercell the 768
components of the HF forces lead to 33 non-zero independent parameters of the
force constants. Knowledge of these force constants allows one to test the
range of interaction potential. For that, similarly to procedure used in
sluiter , we have calculated the absolute values of the eigenvalues of the
$3\times 3$ $\Phi_{ij}({\bf n},\nu;{\bf m},\mu)$ matrices and plotted them
against a distance between atoms $({\bf n},\nu)$ and $({\bf m},\mu)$. The
result for $2$$\times$$2$$\times$$2$ supercell is depicted in Fig. 1. It shows
that the interaction drops down at least by two orders of magnitudes within
the size of the supercell. The scatter of points is too large to assign any
rule which could govern the decaying of the force constants parameters with
distance. From Fig. 1 it is clear that the $1$$\times$$1$$\times$$1$ supercell
is too small to be used for calculation of dispersion curves since the
considered force constants terminate at radius of $a\sqrt{3}/2=0.866a$.

The knowledge of the cummulant force constants allows one to define an
approximate dynamical matrix, in which the summation over all atoms is
confined to those atoms which reside within the volume of the supercell. The
approximate dynamical matrix becomes equal to the conventional one at discrete
wave vectors ${\bf k}_{L}$ given be the equation $\exp(2\pi i{\bf
k}_{L}\cdot{\bf L})=1$ parlinski , where ${\bf L}=({\bf L}_{1},{\bf
L}_{2},{\bf L}_{3})$ are the lattice vectors of the supercell. At the wave
vectors ${\bf k}_{L}$ the phonon frequencies $\omega^{2}({\bf k}_{L})$,
calculated by diagonalization of the approximate dynamical matrix, are the
same as those calculated from the exact dynamical matrix. For
$2$$\times$$2$$\times$$2$ supercell the exact wave vectors are at $\Gamma$, X,
L, W, and two other wave vector points, namely, the midpoint between $\Gamma$
and X along $\langle 1,0,0\rangle$ and $\langle 1,1,1\rangle$ directions. The
advantage of the above described direct method is that it does not impose any
limit to the range of interaction. When the supercell size is smaller than the
range of interaction, the direct method interpolates the dispersion curves
between the exact points.

Figure 2: Phonon dispersion relations of TiC crystal calculated from
$2$$\times$$2$$\times$$2$ supercell (full line). The same dispersion curves
divided by a factor $1.05$ (dashed line). Experimental points taken from Ref.
pint .

The dispersion curves calculated on the basis of $2$$\times$$2$$\times$$2$
supercell are shown in Fig. 2. They correspond to temperature T=0. At the same
Figure 2 we compare the calculated dispersion curves with the experimental
phonon frequencies measured by the inelastic neutron scattering pint .
Different symbols of experimental points correspond to different phonon
polarizations established experimentally.

The experimental points are about $\Delta\omega=0.05\omega$ lower then the
calculated values. One of the reasons for such situation could be that
experimental values are measured at room temperature and anharmonic effects
might deminish the phonon frequencies. Another reason may be related with a
non-stoichiometric concentration of carbon in measured samples, which diminish
in average phonon frequencies as well. The efect that experimental phonon
points are at lower frequencies then the calculated ones appears as well in
other ab-initio calculations frank ; parlinski .

The theory of lattice dynamics provides summation rules which follow from
translational and rotational invariances of the crystal. For TiC these
invariances lead to two equations, which can be used to set the values of on-
site force constant parameters. If these conditions are satisfied the acoustic
phonon branches at k=0 would point to $\omega=0$. Our dispersion curves,
without imposed translational-rotational invariances point at k$=0$ to
$\omega=0.31$THz. Since the invariance conditions should be satisfied within
the ab initio and direct methods, we have not corrected our dispersion curves
displayed in Fig. 2 according to this effect.

Figure 3: Phonon density spectrum $g(\omega)$ of the TiC crystal and partial
phonon density spectra $g_{i,\alpha}(\omega)$, where $\alpha=$ C or Ti and
$i=x,$ $y$ or $z$ are the displacement directions. Experimental points taken
from Ref. pint .

By sampling the dynamical matrix for many wave vectors, one can calculate the
phonon density function $g(\omega)$, and the partial phonon density functions
$g_{x,Ti}(\omega)$ and $g_{x,C}(\omega)$. The $g(\omega)$ describes the number
of phonon frequencies in an interval around $\omega$, while $g_{x,Ti}$ and
$g_{x,C}$ specify the number of phonon frequencies in the interval around
$\omega$ but modified by the population of the $x$ displacement component of
either Ti or C atom. The density of states functions are shown in Fig. 3. They
are conventionally normalized to $\int g(\omega)d\omega=1$ and $\int
g_{x,\alpha}(\omega)d\omega=1/6$, where $\alpha=$ Ti or C. The curves show
that motions within acoustic dispersion curves are almost entirely due to Ti
atoms. The sum of Titanium density functions
($g_{x,Ti}(\omega)+g_{y,Ti}(\omega)+g_{z,Ti}(\omega)$) fits to the total
density function $g(\omega)$ below 14 THz. Part of $g(\omega)$ above the gap
is mainly due to C atoms. Thus, TiC form a rather special crystal, in which
the Ti atoms, as heavier, form a frame for elastic motion, and the C atoms
vibrate within the optical modes. We add, however, that the magnitude of the
force constants between Ti - C and between Ti - Ti and C - C at the same
distances remains of the same order. Thus, TiC does not consist of two weakly
bounded subsystems. In Fig. 3 the experimentally determined phonon density of
states pint ; gompf  is also shown. The agreement within acoustic region is
quite qood. In optic part the experimental resolution of $\approx 2.5$THz and
the nonsoichiometry of carbon in the sample lead to broadening of the optic
band.

The elastic constants can be evaluated from the slopes of acoustic branches of
phonon dispersion curves. Here, we have used the invariance conditions to set
the acoustic modes to $\omega=0$ at k$=0$. For the fcc structure, the set of
seven equations for three unknown variables $c_{ij}$ has to be solved: $\rho
v^{2}_{\alpha,L}=c_{11}$; $\rho v^{2}_{\alpha,T}=c_{44}$; $\rho
v^{2}_{\beta,L}=c_{11}+4c_{44}+2c_{12}$; $\rho
v^{2}_{\beta,T}=c_{11}+c_{44}-c_{12}$; $\rho
v^{2}_{\gamma,L}=c_{11}+2c_{44}+c_{12}$; $\rho
v^{2}_{\gamma,T_{h}}=c_{11}-c_{12}$; $\rho v^{2}_{\gamma,T_{3}}=c_{44}$, where
the velocities $v_{i}$ are the slopes of appropriate acoustic branches. Here,
$\alpha,\beta$ and $\gamma$ indexes describe three main directions in the
reciprocal space $\langle 0,0,1\rangle$, $\langle 1,1,1\rangle$ and $\langle
0,1,1\rangle$, respectively. The density of TiC $\rho=4849.6$ kg/m3 was found
using the calculated lattice constant. The solution of this set of equations,
by means of the least-squares fit, results in $c_{11}=5.7$, $c_{12}=2.5$ and
$c_{44}=1.51$ MBar. Comparing these values with experimental data (Table 2)
one notices considerable discrepancies, especially for $c_{12}$. Closer look
to dispersion curves, Fig. 2, assures us that the transverse acoustic mode in
$\langle 1,1,0\rangle$ direction and the longitudinal in $\langle
1,1,1\rangle$ direction do not follow exactly the experimental points. We have
checked, that such discrepancy of slope is sufficient to produce the observed
deviations of the calculated elastic constants, specially of $c_{12}$. The
discrepancy in slopes is so small that it could be assigned to the missing
force constants describing mainly the long range Coulomb interaction. Indeed,
the direct method did not handle well interactions exceeding the size of the
supercell, except for the special points $\Gamma$, X, L, W mentioned earlier.

Result | $B$ | $c_{11}$ | $c_{12}$ | $c_{44}$
---|---|---|---|---
Exp. TiC0.91 exp_dataA  | $(2.422)$ | $5.145$ | $1.060$ | $1.788$
Exp. TiC exp_dataB  | $2.4$ | $5.00$ | $1.13$ | $1.75$
Exp. TiC0.95 pint  | $(2.5)$ | $5.4$ | $1.1$ | $1.8$
Exp. TiC0.89 pint  | $(2.5)$ | $5.2$ | $1.1$ | $1.8$
Calculations TiC_GGA  | $2.2$ | $4.7$ | $0.97$ | $1.67$
Ours, acoustic modes | $(3.57)$ | $5.7$ | $2.5$ | $1.51$
Ours, deform. energy | $2.42$ | $5.39$ | $0.94$ | $1.92$
Ours, stress–strain | $(2.46)$ | $5.21$ | $1.09$ | $1.93$
Table 2: Experimental and calculated bulk modulus $B$ and elastic constants
$c_{ij}$ of the TiC crystal, in units of $10^{11}{\rm Nm}^{-2}=1{\rm MBar}$.
Values in parenthesis are calculated from $B=\frac{1}{3}(c_{11}+2c_{12})$.

To rule out that the discrepancy for $c_{12}$ does not follow from ab-initio
calculations, we have derived the bulk modulus and the elastic constants by
two conventional methods TiC_GGA . The bulk modulus was calculated in the
$1$$\times$$1$$\times$$1$ supercell from the total energy as a function of the
lattice constant. This data have been fitted to a third order polynomial.
Hence, the derivative $\frac{\partial E}{\partial V}$, and the bulk modulus
$B$ were calculated from the relation: $B=-V\frac{\partial E}{\partial V}$,
where $V$ is the volume of TiC crystallographic unit cell. The $B$ values are
given in Table 2 and they agree quite well with experimental data.

A calculation of $c_{11}$, $c_{12}$, $c_{44}$ elastic constants involved a
similar procedure, in which the crystal was deformed by lengthening and sheer
deformations. Deformations from $0.5\%$ to $3\%$ in length and from $1$ to $5$
degrees in angle have been used. The results of small and large deformations
are consistent. The deformed lattices have space groups Fm$\bar{3}$m, I4/mmm
and I/mmm, for bulk expansion, elongation along z, and sheer modes. In all
these deformed lattices all atoms remain in the high-symmetry sites which
guaranties an extremum of potential energy at those positions. Thus, in the
case of cubic TiC, one is entitled not to carry on the relaxation of internal
degrees of freedom during supercell deformations. The calculated elastic
constants are compared to experimental data in Table 2. They agree with
calculated values given in Ref.TiC_GGA  and are in good agreement with
experiment.

We have also checked the above results using stresses calculated by CASTEP.
The elastic constants were found as first derivatives of the stress–strain
relationship. For that, we have used the same runs of electronic calculations
as for the elastic constant derivation given above using deformation
dependence of energy. The results are included in Table 2 and they show
slightly better agreement with the experimental data.

In summary, we have shown that the ab initio calculations of HF forces,
together with the direct method lead to satisfactory description of phonons in
the TiC crystal. As usual the frequency of the calculated dispersion curves
are slightly too high. The acoustic modes are not sufficiently accurate, thus
the elastic constants calculated from them differ from those found through the
supercell deformations. The acoustic part of density of states fits well to
neutron experimental data.

###### Acknowledgements.

The computational advices of J. Oleszkiewicz and the use of facilities of the
Computer and Network Center of the Jagellonian University, Cracow, where the
calculations have been done, are kindly acknowledged. This work was partially
supported by the Polish State Committee of Scientific Research (KBN), grant No
2 PO3B 145 10.

## References

  * (1) L.E. Toth, Transition Metal Carbides and Nitrides (Academic Press, New York 1971).
  * (2) L. Pintschovius, W. Reichardt and B. Scheerer, J. Phys. C: Solid State Phys. 11, 1557 (1978).
  * (3) R. Chang and L. J. Graham, J. Appl. Phys. 37, 3778 (1966).
  * (4) A. Neckel, R. Eibler, P.Weinberger and K. Schwartz, J. Phys C 9, 579 (1976).
  * (5) V.P. Zhukov, V.A. Gubanov, O. Jepsen, N.E. Christensen and O.K. Andersen, J. Phys. Chem. Solids 49, 841 (1988).
  * (6) D.L. Price and B.R. Cooper, Phys. Rev. B 39, 4945 (1989).
  * (7) J. Häglund, A.F. Guillerment, G. Grimvall and M. Körling, Phys. Rev. B 48, 11685 (1993).
  * (8) R. Ahuja, O. Eriksson, J. M. Wills, B. Johansson, Phys. Rev. B 53, 3073 (1996).
  * (9) W.Frank, C.Elsässer, and M,Fähnle, Phys.Rev.Lett. 74, 1791 (1995).
  * (10) G.Kresse, J.Furthmüller, and J.Hafner, Europhys.Lett. 32, 729 (1995).
  * (11) K. Parlinski, Z. Q. Li and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997).
  * (12) G. J. Ackland, M. C. Warren and S. J. Clark, J. Phys.: Condens. Matter 9, 7861–7872 (1997).
  * (13) M.H.F. Sluiter, M. Weinert, Y. Kawazoe, Europhys. Lett. 43, 183 (1998).
  * (14) J.S. Lin, A. Qteish, M.C. Payne and V. Heine, Phys. Rev. B 47, 4174 (1993).
  * (15) J. Goniakowski, J.M. Holender, L.N. Kantorovich, M.J. Gillan and J.A. White, Phys. Rev. B 53, 953 (1996).
  * (16) M.C. Payene, M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopoulos, Rev. of Mod. Phys. 64, 1045 (1992)
  * (17) CASTEP is a part of $Cerius^{2}$ system developed by Molecular Simulations.
  * (18) J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh, C. Fiolhais, Phys. Rev. B 46, 6671 (1992).
  * (19) J. Goniakowski, J. M. Holender, L. N. Kantorovich, M. J. Gillan, J. A. White, Phys. Rev. B 53, 957 (1996)
  * (20) H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188 (1976).
  * (21) I. Náray-Szabó, Inorganic Crystal Chemistry (Akadémiai kiadó, Budapest, 1969), p.301.
  * (22) A. García, C. Elsässer, J. Zhu, S.G. Louie, M.L. Cohen, Phys. Rev. B 46, 9829 (1992)
  * (23) Y.M. Juan and E. Kaxiras, Phys. Rev. B 48, 14 994 (1993)
  * (24) K. Seifert, J. Hafner, J. Furthmüller, G. Kresse, J. Phys. Condens. Matter 7, 3683 (1995)
  * (25) W. H. Press, B.P. Flamery, S.A. Tenkolsky, and W.T. Vetterling, Numerical Recipes (Cambridge Univ. Press, 1988), p.301
  * (26) F. Gompf and W. Reichardt, 1977, unpublished.
  * (27) M. M. Choy at. al., Numerical Data and Functional Relationships in Science and Technology (Springer, Berlin, 1979), vol.11.