#include <bits/stdc++.h>
const int mod=1e9+7;
using namespace std;
int n,head[200010],fa[200010],o,m;
int mx[200010],se[200010],th[200010],lim[200010],sz[200010],f[200010],ans;
struct edge{
	int to,link;
}e[400010];
void qmo(int &x){
	x+=(x>>31)&mod;
}
void add_edge(int u,int v){
	e[++o]={v,head[u]},head[u]=o;
	e[++o]={u,head[v]},head[v]=o;
}
void dfs(int u,int pre){
	sz[u]=1,fa[u]=pre;
	mx[u]=se[u]=th[u]=-2e9;
	for(int i=head[u],v;i;i=e[i].link) if((v=e[i].to)^pre){
		dfs(v,u),sz[u]+=sz[v];
		if(sz[v]>=mx[u]) th[u]=se[u],se[u]=mx[u],mx[u]=sz[v];
		else if(sz[v]>=se[u]) th[u]=se[u],se[u]=sz[v];
		else if(sz[v]>th[u]) th[u]=sz[v];
	}
	if(pre){
		if(n-sz[u]>=mx[u]) th[u]=se[u],se[u]=mx[u],mx[u]=n-sz[u];
		else if(n-sz[u]>=se[u]) th[u]=se[u],se[u]=n-sz[u];
		else if(n-sz[u]>th[u]) th[u]=n-sz[u];
	}
	m=min(m,mx[u]);
	if(se[u]>0) lim[se[u]]=max(lim[se[u]],th[u]);
} 
int get(int wl,int bl){//w<=wl,b<=bl,w>=b
	bl=min(bl,wl);
	if(bl==wl){
		return f[wl];
	}
	int res=bl*(bl+1ll)%mod*((wl-bl)*(wl+bl+1ll)/2%mod)%mod;
	qmo(res+=get(bl,bl)-mod);
	return res;
}
void sol(){
	scanf("%d",&n),o=0,m=2e9,ans=0;
	for(int i=1;i<=n;i++) head[i]=0;
	for(int i=2,w;i<=n;i++) scanf("%d",&w),add_edge(w,i);
	for(int i=1;i<=n;i++) f[i]=(f[i-1]+1ll*i*i%mod*i)%mod,lim[i]=0;
	dfs(1,0);
	for(int i=n;i>=1;i--) lim[i]=max(lim[i+1],lim[i]);
	if(m>=1){
		int w=m*(m+1ll)/2%mod;
		ans=2ll*w*w%mod;
		for(int w=1;w<=m;w++){
			int b=min(w,lim[w]);
			qmo(ans-=b*(b+1ll)/2%mod*w%mod);
			b=min(w-1,lim[w]);
			qmo(ans-=b*(b+1ll)/2%mod*w%mod);
		}
	}
	for(int u=1;u<=n;u++){
		int mxx=0,Up=-1,o;
		for(int i=head[u],v;i;i=e[i].link) if(mx[v=e[i].to]<=mx[u]&&mx[v]>mxx) Up=v,mxx=mx[v],o=(v==fa[u]?sz[u]:n-sz[v]);
		mxx=max(mxx,m);
		if(mxx<mx[u]&&Up>0) qmo(ans+=get(mx[u],o)-mod),qmo(ans-=get(mxx,o));
	}
	printf("%d\n",ans);
}
int main(){
	int _;
	scanf("%d",&_);
	while(_--) sol();

}