#include <bits/stdc++.h>
#define int long long
#define S(X, Y) (((X) + (Y)) * ((Y) - (X) + 1) / 2) % Mod
using namespace std;

const int N = 2e5 + 10, Mod = 1e9 + 7;
int n, fa[N], nxt[N], siz[N], ans;
pair<int, int> Max;
vector<int> g[N];

void cal(int u, int pre){
	int fir = 0, sec = 0, tmp;
	for(int v : g[u]) if(v != pre){
		tmp = siz[v];
		if(tmp >= fir) swap(fir, tmp);
		if(tmp >= sec) swap(sec, tmp);
	}
	nxt[fir] = max(nxt[fir], sec);
	nxt[sec] = max(nxt[sec], fir);
	for(int v : g[u]) if(v != pre)
		cal(v, u);
	return ;
};

void fun(int u, int pre){
	int tmp = S(max(n - siz[pre] + 1, Max.first + 1), n - siz[u]);
	tmp = 2 * tmp % Mod * S(1, siz[u]) % Mod; ans = (ans + tmp) % Mod;
	for(int v : g[u]) if(v != pre) fun(v, u);
	return ;
};

void dfs(int u, int pre){
	siz[u] = 1;
	for(int v : g[u]) if(v != pre)
		dfs(v, u), siz[u] += siz[v];
	return ;
};

void init(){
	Max = make_pair(0, 0);
	for(int i = 0; i <= n + 5; i++)
		siz[i] = fa[i] = nxt[i] = 0;
	ans = 0;
	for(int i = 0; i < N; i++) g[i].clear();
}

void solve(){
	init();

	cin >> n;
	for(int i = 2; i <= n; i++){
		cin >> fa[i];
		g[i].push_back(fa[i]);
		g[fa[i]].push_back(i);
	}
	for(int i = n; i >= 1; i--) siz[fa[i]] += ++siz[i];
	int G = 1, H = 0;

	find_G:;
	for(int v : g[G]) if(v != H)
		if(siz[v] > siz[1] / 2){
			H = G, G = v; goto find_G;}

	dfs(G, 0); siz[G]--;

	for(int v : g[G]) Max = max(Max, make_pair(siz[v], v));

	ans = S(1, Max.first) * S(1, Max.first) % Mod;

	cal(G, Max.second);
	cal(Max.second, G);
	for(int i = Max.first; i >= 1; i--){
		nxt[i] = max(nxt[i], nxt[i + 1]);
		if(nxt[i] + 1 <= Max.first)
			ans = (ans + i * S(nxt[i] + 1, Max.first)) % Mod;
	}

	for(int v : g[G]) fun(v, G);

	cout << ans << '\n';
	return ;
}

signed main(){
	int T; cin >> T;
	while(T--) solve();
	return 0;
}