Adding on to Bob's good advice, you can use 3dAllineate to compute the "in-between" space too. cat_matvec can compute this with from the affine transformation with the -S (square root) option. We have a script here that does something similar for a single subject across two sessions:

https://sscc.nimh.nih.gov/sscc/dglen/alignmentacross2sessions

That in-between volume is meant to remove bias from either of the source sessions. This part only concerns the affine part of the alignment. For the nonlinear part, as Bob noted, you will need to use 3dQwarp.

Another important consideration is to make the data intensity values similar if you are computing a mean "template" with 3dUnifize. Most of this kind of thing is done in make_template_dask.py (a successor for the similar scripts you mentioned), but that's not distributed in our standard AFNI distribution yet, and it's probably overkill for this similar but somewhat simpler situation.

If sizes are very similar, you might consider using a rigid_body or rigid_equiv version of the affine transformation. These are both options for align_epi_anat.py that includes calls to 3dAllineate. The rigid_body limits the alignment by using 3dAllineate's "-shift_rotate" option. The rigid equivalent also produces a transformation that does not scale or shear the dataset, but it's derived from the full affine transformation (compute with cat_matvec -P or use align_epi_anat.py). With data from different subjects, it's likely this would be preferable.

Twins, even "identical" twins, aren't really identical, so it's hard to say if their brains would be a very close match. Monochorionic twins, for example, can have very different anatomical characteristics.