$\def\cX{\mathcal X}\def\cY{\mathcal Y}\def\Th{\mathrm{Th}}$UPDATE: The question was changed so that $\Th_2(\cX)$ no longer refers to the second-order theory of $\cX$, but the second-order diagram. This likely means the answer below could be considerbly simplified, as the answer to the secondary question basically already answers the main one. But I am not going to rewrite it.
The answer to the secondary question is an easy YES: for example, given an $L$-structure $\cX$ with domain $[n]$, let
$F(\cX)$ be the collection of pairs $(\cY,i)$ such that $\cY$ is an $L$-structure on $[n]$ isomorphic to $\cX$, $i\le n\log n$ is an integer, and the $i$th bit of the description of the lexicographically first isomorphism between $\cX$ and $\cY$ is $1$. Note that it is enough to have oracle access to one of $F(\cX)$ and $F(\cY)$; we do not need both.
If we want to make it depend less on details of the representation of the output, let e.g. $F(\cX)$ consist of pairs $(\cY,f)$ such that $\cY$ is an $L$-structure on $[n]$,
$f$ is a partial self-map of $[n]$, and there exists an isomorphism $F\colon\cX\simeq\cY$ such that $F\supseteq f$. (It suffices to consider $f$ being the identity on an initial segment of the domain.) Using this, we can construct an isomorphism between $\cX$ and $\cY$ by fixing it one element at a time, querying, for a given $a\notin\operatorname{dom}(f)$, the existence of an isomorphism extending $f\cup\{(a,b)\}$ for each $b$ until we are successful.
We can express this by the truth of an existential first-order formula: if $f=\{(a_i,b_i):i<k\}$, let $\delta(\vec b,\vec y)$ be the (conjunction of the) atomic diagram of $(\cY,\vec b)$; then $f$ extends to an isomorphism $\cX\simeq\cY$ iff $\mathcal X\models\exists\vec y\,\delta(\vec a,\vec y)$. Thus, you can also take for $F(\mathcal X)$ a truth predicate for existential first-order formulas with parameters from $\cX$.
The same kind of argument also works for the main question: given structures $\cX$ and $\cY$ with domain $[n]$, and a partial self-map $f$ on $[n]$, let $\delta_\cX(\vec x)$ be the atomic diagram of $\mathcal X$ in variables $\{x_i:i\in[n]\}$, and $\delta_\cY(\vec y)$ be the atomic diagram of $\cY$ in variables $\{y_i:i\in[n]\}$. Then $f$ extends to an isomorphism $\cX\simeq\cY$ iff the existential second-order sentence (with no extra-logical symbols)
$$\exists F\,\exists\vec R\,\exists\vec x\,\Bigl(\delta_\cX(\vec x)\land\delta_\cY(F(\vec x))\land\bigwedge_{f(i)=j}F(x_i)=x_j\Bigr)$$
is true in some/every structure with $n$ elements, where $\vec R$ quantifies over all relations and functions of $L$ as used in $\delta_\cX$ and $\delta_\cY$, and $F(\vec x)$ denotes the sequence of values $F(x_0),F(x_1),\dots$. Again, using this oracle, we can construct an isomorphism $\cX\simeq\cY$ one element at a time; we only use the oracle for one of $\Th_2(\cX)$ or $\Th_2(\cY)$ rather than their join. (Actually, the join is pointness anyway: since $\cX\simeq\cY$, we have $\Th_2(\cX)=\Th_2(\cY)$, and $\Th_2(\cX)\oplus\Th_2(\cY)$ is trivially equivalent to $\Th_2(\cX)$.)
The only issue is whether this can be done in linear time. This is quite a tall order. Since it may depend on all kinds of details, let me state for definiteness that we are using the standard multi-tape Turing machine model with a read-only input tape, write-only output tape, several work tapes, and an oracle query tape. The content of the oracle query tape is not erased or otherwise modified after making an oracle query. I’m assuming $\cX$ and $\cY$ are presented as input in the most obvious way as lists of tables of their relations and functions. Since $\Th_2(\cX)$ has to be presented as a set of strings over a finite alphabet, variables (both first- and second-order) in the formulas have to be identified using numerical indices; I will assue these are written in binary, so that if I use, say, $n$ variables, each takes $O(\log n)$ bits to write down.
Note that if we have only unary relations and constants, then the size of the input is $O(n)$, whereas the size of the output is $\Omega(n\log n)$. Thus, the task is impossible to do in linear time with whatever oracle, as a linear-time machine does not even have the time to write down the result.
I do not know how to implement the algorithm above in linear time, but as outlined below, it can be done with an extra $\log n$ factor (or a little bit more in the purely unary case).
The diagrams $\delta_\cX$ and $\delta_\cY$ are basically just copies of the tables of the relations and functions of $\cX$ and $\cY$, thus we can write down the $\exists F\,\exists\vec R\,\exists\vec x\,\bigl(\delta_\cX(\vec x)\land\delta_\cY(F(\vec x))\land\dots$ part of the formula on the query tape in time $O\bigl((|\cX|+|\cY|)\log n)$, where $|\cX|$ denotes the size of $\cX$ as an input string; the $\log n$ factor comes from variable indices of $x_i$ and $y_i$. (Actually, we also need to include the constraints that the $x_i$ are pairwise distinct, and likewise for $\vec y$. The obvious way, using a conjunction of $\binom n2$ inequalities, takes $O(n^2\log n)$ bits. We can do it more efficiently, but let’s not worry about that for now.)
Then we write down $\dots F(x_0)=x_0\bigr)$ and query the oracle; if the answer is negative, we modify it to $\dots F(x_0)=x_1\bigr)$, and so on, until we get a positive answer. Each such change takes amortized constant time by the standard analysis of a binary counter; that is, it takes time $O(n)$ to find the image of $0$. We copy it to the output tape, and go on to add $\dots\land F(x_1)=x_0$ on the query tape, and continue. In this way, we construct the isomorphism in time $O(n^2)$ on top of the $O\bigl((|\cX|+|\cY|)\log n)$ initialization part.
Thus, if $L$ contains at least one at least binary relation or function symbol, the algorithm as a whole works in time $O\bigl((|\cX|+|\cY|)\log n)$, as the size of the input, $|\cX|+|\cY|$, is $\Omega(n^2)$.
If all symbols in $L$ are unary, we need to work harder. As mentioned above, the best we can hope for in this case is to make it work in time $O(n\log n)$; that is, to have an algorithm that works in time linear in the combined size of the input and output. I’m not sure whether that is possible to arrange; the best I can do is $O(n\log n\log\log n)$.
We can reduce the number of oracle queries to $O(n\log n)$ if we use binary search to find the image of $i\in[n]$: i.e., instead of trying whether the image can be $0$, $1$, $2$, ..., we first determine the first bit of the image, then the second bit, and so on. To make sure that the individual queries do not take too much time to construct, we can introduce predicates for elements with a given bit set: the formula will now start with
$$\exists F\,\exists\vec R\,\exists\vec x\,\exists P_0,\dots,P_{\log n}\,\Bigl(\delta_\cX(\vec x)\land\delta_\cY(F(\vec x))\land\bigwedge_{\substack{i\in[n]\\j\le\log n}}(P_j(x_i))^{b(i,j)}\land\dots,$$
where $b(i,j)\in\{0,1\}$ is the $j$th bit of $i$, and $\phi^1=\phi$, $\phi^0=\neg\phi$. (NB: This also implies distinctness of the $x_i$s.) The indices of the predicate variables in the $n\log n$ new conjuncts take $O(n\log n\log\log n)$ bits. The indices of the $x_i$ variables would take $O(n(\log n)^2)$ bits as written, but we can reduce this to $O(n\log n)$ by judicious quantification: instead of $\bigwedge_j(P_j(x_i))^{b(i,j)}$, write $\exists z\,\bigl(z=x_i\land\bigwedge_j(P_j(z))^{b(i,j)}\bigr)$. One can then arrange the whole computation to take time $O(n\log n\log\log n)$.
